program test_calculate

use Compact

implicit none

real*8, parameter :: pi=3.14159265358979323846d0
real*8, parameter :: me=1.00d-20

integer :: nxp,nyp,nzp

real*8 :: hxp,hyp,hzp

real*8 :: l_0

real*8, allocatable, dimension(:,:,:) :: aaa,bbb,ccc,ddd

real*8, allocatable, dimension(:,:) :: aa,bb,cc,dd,ee

integer :: ierr

real*8 :: xx,yy,zz,fxyz

integer :: i,j,k

real*4, dimension(2) :: startt
real*4, dimension(2) :: endt
real*4 :: result

nxp = 512
nyp = 256
nzp = 256

l_0 = 2.0d0
hyp=l_0*pi/REAL(nyp)
hxp=hyp
hzp=hyp

allocate(aa(1, nxp), stat=ierr); aa=0.0
allocate(bb(1, nxp), stat=ierr); bb=0.0
allocate(cc(1, nxp), stat=ierr); cc=0.0
allocate(dd(1, nxp), stat=ierr); dd=0.0
allocate(ee(1, nxp), stat=ierr); ee=0.0


allocate(aaa(nxp, nyp, nzp), stat=ierr); aaa=0.0
allocate(bbb(nxp, nyp, nzp), stat=ierr); bbb=0.0
allocate(ccc(nxp, nyp, nzp), stat=ierr); ccc=0.0
allocate(ddd(nxp, nyp, nzp), stat=ierr); ddd=0.0


call ludcmp(nxp, nyp, nzp, 1, 0, 0)

do i = 1,nxp
xx = dfloat(i-1) * hxp
aa(1,i) = sin(1.1d0 * xx)
dd(1,i) = 1.1d0 * cos(1.1d0 * xx)
end do

call rhs1np(nxp, hxp, aa, bb, 1)

call dfnonp(nxp, hxp, aa, cc, 1, 1)


print *, lxf

print *, lxs

print *, aa

print *, bb

print *, cc

print *, dd

ee = ((cc - dd)/dd)


write(*,*) "relerr_min_max", minval(ee(1,nxp/4:3*nxp/4)), maxval(ee(1,nxp/4:3*nxp/4))
write(*,*) "relerr_min_max", minval(ee), maxval(ee)
write(*,*) "  exact_min_max", minval(dd), maxval(dd)
write(*,*) "compact_min_max", minval(cc), maxval(cc)


write(*,*) "test ddx"

do k = 1,nzp
do j = 1,nyp
do i = 1,nxp

xx = i * hxp
yy = j * hyp
zz = k * hzp

aaa(i,j,k) = sin(1.1 * xx) * sin(3.0 * yy) * sin(2.0 * zz)

end do
end do
end do





do k = 1,nzp
do j = 1,nyp

call dfnonp(nxp, hxp, aaa(:,j,k), bbb(:,j,k), 1, 1)

end do
end do

do k = 1,nzp
do j = 1,nyp
do i = 1,nxp ! nxp/4,3*nxp/4

xx = i * hxp
yy = j * hyp
zz = k * hzp
fxyz = (1.1 * cos(1.1 * xx) * sin(3.0 * yy) * sin(2.0 * zz))
ccc(i,j,k) = (bbb(i,j,k) - fxyz) / (fxyz)

if ((abs(ccc(i,j,k)) > 0.01 ) .and. (abs(bbb(i,j,k)) > 1.0e-14 )) write(*,*) bbb(i,j,k), fxyz

end do
end do
end do

write(*,*) "relerr_min_max", minval(ccc), maxval(ccc)







deallocate(aa)
deallocate(bb)
deallocate(cc)
deallocate(dd)
deallocate(ee)

deallocate(aaa)
deallocate(bbb)
deallocate(ccc)
deallocate(ddd)

end program