lxsd / fvn

fvn_test/Makefile

	1	+
	2	+include $(BTREE)/Make.inc
	3	+
	4	+programs = test_fac test_matinv test_specfunc \
	5	+test_det test_matcon test_matev test_sparse test_inter1d \
	6	+test_inter2d test_inter3d test_akima test_lsp test_muller \
	7	+test_integ
	8	+
	9	+prog:$(programs)
	10	+
	11	+clean:
	12	+ rm -f {.o,.oo,.ipo,.a,*.mod}
	13	+ rm -f $(programs)
	14	+
	15	+%: %.o
	16	+ $(LINK) $(LINKFLAGS) $< init_random_seed.o $(LINKFVN) -o $@
	17	+
	18	+%.o: %.f90
	19	+ $(F95) $(F95FLAGS) -c $<
	20	+
	21	+$(programs): init_random_seed.o

fvn_test/init_random_seed.f90

Diff comments View file @ 27d3b84

	1	+! As the norm stated for random_seed "If no argument is present, the processor assigns a
	2	+! processor-dependent value to the seed."
	3	+! It seems that there are different possible interpretations for this :
	4	+! In gfortran a call to random_seed without arguments will initialize to a default value
	5	+! which will always be the same.
	6	+! In Portland fortran, a call to random_seed without arguments will initialize with something
	7	+! machine and time dependant
	8	+! Unfortunately as the underlying algorithm are (or at least can be) different in each
	9	+! compiler, we cannot generalize the use of this initialization provided in gfortran
	10	+! documentation
	11	+SUBROUTINE init_random_seed()
	12	+ INTEGER :: i, n, clock
	13	+ INTEGER, DIMENSION(:), ALLOCATABLE :: seed
	14	+ logical :: gfortran=.false.
	15	+ !
	16	+ ! Set gfortran to .true. if using gfortran to obtain the desired behaviour
	17	+ !
	18	+ if (gfortran) then
	19	+ CALL RANDOM_SEED(size = n)
	20	+ ALLOCATE(seed(n))
	21	+
	22	+ CALL SYSTEM_CLOCK(COUNT=clock)
	23	+
	24	+ seed = clock + 37 * (/ (i - 1, i = 1, n) /)
	25	+ CALL RANDOM_SEED(PUT = seed)
	26	+
	27	+ DEALLOCATE(seed)
	28	+ else
	29	+ call random_seed()
	30	+ end if
	31	+END SUBROUTINE

fvn_test/test_akima.f90

Diff comments View file @ 27d3b84

	1	+program akima
	2	+ use fvn
	3	+ implicit none
	4	+ integer :: nbpoints,nppoints,i
	5	+ real(8),dimension(:),allocatable :: x_d,y_d,breakpoints_d
	6	+ real(8),dimension(:,:),allocatable :: coeff_fvn_d
	7	+ real(8) :: xstep_d,xp_d,ty_d,fvn_y_d
	8	+ open(2,file='fvn_akima_double.dat')
	9	+ open(3,file='fvn_akima_breakpoints_double.dat')
	10	+ nbpoints=30
	11	+ allocate(x_d(nbpoints))
	12	+ allocate(y_d(nbpoints))
	13	+ allocate(breakpoints_d(nbpoints))
	14	+ allocate(coeff_fvn_d(4,nbpoints))
	15	+ xstep_d=20./dfloat(nbpoints)
	16	+ do i=1,nbpoints
	17	+ x_d(i)=-10.+dfloat(i)*xstep_d
	18	+ y_d(i)=dsin(x_d(i))
	19	+ write(3,44) x_d(i),y_d(i)
	20	+ end do
	21	+ close(3)
	22	+ call fvn_akima(nbpoints,x_d,y_d,breakpoints_d,coeff_fvn_d)
	23	+ nppoints=1000
	24	+ xstep_d=22./dfloat(nppoints)
	25	+ do i=1,nppoints
	26	+ xp_d=-11.+dfloat(i)*xstep_d
	27	+ ty_d=dsin(xp_d)
	28	+ fvn_y_d=fvn_spline_eval(xp_d,nbpoints-1,breakpoints_d,coeff_fvn_d)
	29	+ write(2,44) xp_d,ty_d,fvn_y_d
	30	+ end do
	31	+ close(2)
	32	+deallocate(coeff_fvn_d,breakpoints_d,y_d,x_d)
	33	+write(,) "All done, plot results with gnuplot using command :"
	34	+write(,) "pl 'fvn_akima_double.dat' u 1:2 w l,'fvn_akima_breakpoints_double.dat' w p"
	35	+
	36	+44 FORMAT(4(1X,1PE22.14))
	37	+end program

fvn_test/test_det.f90

Diff comments View file @ 27d3b84

	1	+program test_det
	2	+use fvn
	3	+implicit none
	4	+real(8),dimension(3,3) :: a
	5	+real(8) :: deta
	6	+integer :: status,i
	7	+
	8	+call init_random_seed()
	9	+call random_number(a)
	10	+a=a*100
	11	+deta=fvn_det(3,a,status)
	12	+do i=1,3
	13	+ write (*,'(3(f15.5))') a(i,:)
	14	+end do
	15	+write (,)
	16	+write (,) "Det = ",deta,"Status = ",status
	17	+end program

fvn_test/test_integ.f90

Diff comments View file @ 27d3b84

	1	+program integ
	2	+ use fvn
	3	+ implicit none
	4	+ real(8), external :: f1,f2,g,h
	5	+ real(8) :: a,b,epsabs,epsrel,abserr,res
	6	+ integer :: key,ier
	7	+ a=0.
	8	+ b=1.
	9	+ epsabs=1d-8
	10	+ epsrel=1d-8
	11	+ key=2
	12	+ call fvn_integ_1_gk(f1,a,b,epsabs,epsrel,key,res,abserr,ier)
	13	+ write(,) "Integration of x*x between 0 and 1 : "
	14	+ write(,) res
	15	+ call fvn_integ_2_gk(f2,a,b,g,h,epsabs,epsrel,key,res,abserr,ier)
	16	+ write(,) "Integration of x*y between 0 and 1 on both x and y : "
	17	+ write(,) res
	18	+
	19	+end program
	20	+function f1(x)
	21	+ implicit none
	22	+ real(8) :: x,f1
	23	+ f1=x*x
	24	+end function
	25	+function f2(x,y)
	26	+ implicit none
	27	+ real(8) :: x,y,f2
	28	+ f2=x*y
	29	+end function
	30	+function g(x)
	31	+ implicit none
	32	+ real(8) :: x,g
	33	+ g=0.
	34	+end function
	35	+function h(x)
	36	+ implicit none
	37	+ real(8) :: x,h
	38	+ h=1.
	39	+end function

fvn_test/test_inter1d.f90

Diff comments View file @ 27d3b84

	1	+program inter1d
	2	+use fvn
	3	+implicit none
	4	+integer(kind=4),parameter :: ndata=33
	5	+integer(kind=4) :: i,nout
	6	+real(kind=8) :: f,fdata(ndata),h,pi,q,sin,x,xdata(ndata)
	7	+real(kind=8) ::tv
	8	+intrinsic sin
	9	+f(x)=sin(x)
	10	+xdata(1)=0.
	11	+fdata(1)=f(xdata(1))
	12	+h=1./32.
	13	+do i=2,ndata
	14	+xdata(i)=xdata(i-1)+h
	15	+fdata(i)=f(xdata(i))
	16	+end do
	17	+call init_random_seed()
	18	+call random_number(x)
	19	+q=fvn_quad_interpol(x,ndata,xdata,fdata)
	20	+tv=f(x)
	21	+write(*,'("x y z ",1(f8.5))') x
	22	+write(*,'("Calculated (real) value :",f8.5)') tv
	23	+write(*,'("fvn interpolation : ",f8.5)') q
	24	+write(*,'("Relative fvn error :",e12.5)') abs((q-tv)/tv)
	25	+end program

fvn_test/test_inter2d.f90

Diff comments View file @ 27d3b84

	1	+program inter2d
	2	+use fvn
	3	+implicit none
	4	+integer(kind=4),parameter :: nx=21,ny=42
	5	+integer(kind=4) :: i,j
	6	+real(kind=8) :: f,fdata(nx,ny),dble,pi,q,sin,x,xdata(nx),y,ydata(ny)
	7	+real(kind=8) :: tv
	8	+intrinsic dble,sin
	9	+f(x,y)=sin(x+2.*y)
	10	+do i=1,nx
	11	+xdata(i)=dble(i-1)/dble(nx-1)
	12	+end do
	13	+do i=1,ny
	14	+ydata(i)=dble(i-1)/dble(ny-1)
	15	+end do
	16	+do i=1,nx
	17	+do j=1,ny
	18	+fdata(i,j)=f(xdata(i),ydata(j))
	19	+end do
	20	+end do
	21	+call init_random_seed()
	22	+call random_number(x)
	23	+call random_number(y)
	24	+q=fvn_quad_2d_interpol(x,y,nx,xdata,ny,ydata,fdata)
	25	+tv=f(x,y)
	26	+write(*,'("x y z ",2(f8.5))') x,y
	27	+write(*,'("Calculated (real) value :",f8.5)') tv
	28	+write(*,'("fvn interpolation : ",f8.5)') q
	29	+write(*,'("Relative fvn error :",e12.5)') abs((q-tv)/tv)
	30	+end program

fvn_test/test_inter3d.f90

Diff comments View file @ 27d3b84

	1	+program test_inter3d
	2	+use fvn
	3	+implicit none
	4	+integer(kind=4),parameter :: nx=21,ny=42,nz=18
	5	+integer(kind=4) :: i,j,k
	6	+real(kind=8) :: f,fdata(nx,ny,nz),dble,pi,q,sin,x,xdata(nx),y,ydata(ny),z,zdata(nz)
	7	+real(kind=8) :: tv
	8	+intrinsic dble,sin
	9	+f(x,y,z)=sin(x+2.y+3.z)
	10	+do i=1,nx
	11	+ xdata(i)=2.*(dble(i-1)/dble(nx-1))
	12	+end do
	13	+do i=1,ny
	14	+ ydata(i)=2.*(dble(i-1)/dble(ny-1))
	15	+end do
	16	+do i=1,nz
	17	+ zdata(i)=2.*(dble(i-1)/dble(nz-1))
	18	+end do
	19	+do i=1,nx
	20	+ do j=1,ny
	21	+ do k=1,nz
	22	+ fdata(i,j,k)=f(xdata(i),ydata(j),zdata(k))
	23	+ end do
	24	+ end do
	25	+end do
	26	+call init_random_seed()
	27	+call random_number(x)
	28	+call random_number(y)
	29	+call random_number(z)
	30	+q=fvn_quad_3d_interpol(x,y,z,nx,xdata,ny,ydata,nz,zdata,fdata)
	31	+tv=f(x,y,z)
	32	+write(*,'("x y z ",3(f8.5))') x,y,z
	33	+write(*,'("Calculated (real) value :",f8.5)') tv
	34	+write(*,'("fvn interpolation : ",f8.5)') q
	35	+write(*,'("Relative fvn error :",e12.5)') abs((q-tv)/tv)
	36	+end program

fvn_test/test_lsp.f90

Diff comments View file @ 27d3b84

	1	+program lsp
	2	+use fvn
	3	+implicit none
	4	+integer,parameter :: npoints=13,deg=3
	5	+ integer :: status,i
	6	+ real(kind=8) :: xm(npoints),ym(npoints),xstep,xc,yc
	7	+ real(kind=8) :: coeff(deg+1)
	8	+ xm = (/ -3.8,-2.7,-2.2,-1.9,-1.1,-0.7,0.5,1.7,2.,2.8,3.2,3.8,4. /)
	9	+ ym = (/ -3.1,-2.,-0.9,0.8,1.8,0.4,2.1,1.8,3.2,2.8,3.9,5.2,7.5 /)
	10	+ open(2,file='fvn_lsp_double_mesure.dat')
	11	+ open(3,file='fvn_lsp_double_poly.dat')
	12	+ do i=1,npoints
	13	+ write(2,44) xm(i),ym(i)
	14	+ end do
	15	+ close(2)
	16	+ call fvn_lspoly(npoints,xm,ym,deg,coeff,status)
	17	+ xstep=(xm(npoints)-xm(1))/1000.
	18	+ do i=1,1000
	19	+ xc=xm(1)+(i-1)*xstep
	20	+ yc=poly(xc,coeff)
	21	+ write(3,44) xc,yc
	22	+ end do
	23	+ close(3)
	24	+write(,) "All done, plot results with gnuplot using command :"
	25	+write(,) "pl 'fvn_lsp_double_mesure.dat' u 1:2 w p,'fvn_lsp_double_poly.dat' u 1:2 w l"
	26	+44 FORMAT(4(1X,1PE22.14))
	27	+contains
	28	+function poly(x,coeff)
	29	+ implicit none
	30	+ real(8) :: x
	31	+ real(8) :: coeff(deg+1)
	32	+ real(8) :: poly
	33	+ integer :: i
	34	+ poly=0.
	35	+ do i=1,deg+1
	36	+ poly=poly+coeff(i)x*(i-1)
	37	+ end do
	38	+end function
	39	+end program

fvn_test/test_matcon.f90

Diff comments View file @ 27d3b84

	1	+program test_matcon
	2	+use fvn
	3	+implicit none
	4	+real(8),dimension(3,3) :: a
	5	+real(8) :: rcond
	6	+integer :: status,i
	7	+call init_random_seed()
	8	+call random_number(a)
	9	+a=a*100
	10	+call fvn_matcon(3,a,rcond,status)
	11	+write(,) "Reasonnably conditionned matrix"
	12	+do i=1,3
	13	+ write (*,'(3(e12.5))') a(i,:)
	14	+end do
	15	+write (,)
	16	+write (,) "Cond = ",rcond
	17	+write (,)
	18	+write (,)
	19	+a(1,1)=a(1,1)*1d9
	20	+write(,) "Badly conditionned matrix"
	21	+do i=1,3
	22	+ write (*,'(3(e12.5))') a(i,:)
	23	+end do
	24	+call fvn_matcon(3,a,rcond,status)
	25	+write (,)
	26	+write (,) "Cond = ",rcond
	27	+
	28	+end program

fvn_test/test_matev.f90

Diff comments View file @ 27d3b84

	1	+program test_matev
	2	+use fvn
	3	+implicit none
	4	+real(8),dimension(3,3) :: a
	5	+complex(8),dimension(3) :: evala
	6	+complex(8),dimension(3,3) :: eveca
	7	+integer :: status,i,j
	8	+call init_random_seed()
	9	+call random_number(a)
	10	+a=a*100
	11	+call fvn_matev(3,a,evala,eveca,status)
	12	+
	13	+write(,) "The matrix :"
	14	+write (*,'(3(e12.5))') a
	15	+write (,)
	16	+do i=1,3
	17	+write(*,'("Eigenvalue ",I3," : (",e12.5,",",e12.5,") ")') i,evala(i)
	18	+write(,) "Associated Eigenvector :"
	19	+do j=1,3
	20	+write(*,'("(",e12.5,",",e12.5,") ")') eveca(j,i)
	21	+end do
	22	+write(,)
	23	+end do
	24	+end program

fvn_test/test_matinv.f90

Diff comments View file @ 27d3b84

	1	+program test_matinv
	2	+use fvn
	3	+implicit none
	4	+integer(4), parameter :: n=3
	5	+integer(4) :: status,i
	6	+complex(8),dimension(n,n) :: m,im,prod
	7	+real(8),dimension(n,n) :: rtmp,itmp
	8	+character(len=80) :: fmreal,fmcmplx
	9	+
	10	+fmcmplx='(3("(",f8.5,",",f8.5,") "))'
	11	+
	12	+! initialize pseudo random generator
	13	+call init_random_seed()
	14	+! fill real and imaginary part
	15	+call random_number(rtmp)
	16	+call random_number(itmp)
	17	+! create the complex matrix (fvn_i is defined in the fvn module)
	18	+m=rtmp+fvn_i*itmp
	19	+write(,) "Matrix M"
	20	+do i=1,n
	21	+ write(*,fmcmplx) m(i,:)
	22	+end do
	23	+
	24	+! Invertion
	25	+call fvn_matinv(n,m,im)
	26	+write(,) "Inverse of M"
	27	+do i=1,n
	28	+ write(*,fmcmplx) im(i,:)
	29	+end do
	30	+
	31	+! Result should be identity matrix
	32	+write(,) "Product of M and inverse of M :"
	33	+prod=matmul(m,im)
	34	+do i=1,n
	35	+ write (*,fmcmplx) prod(i,:)
	36	+end do
	37	+
	38	+end program

fvn_test/test_muller.f90

Diff comments View file @ 27d3b84

	1	+program muller
	2	+use fvn
	3	+implicit none
	4	+integer :: i,info
	5	+complex(8),dimension(10) :: roots
	6	+integer,dimension(10) :: infer
	7	+complex(8), external :: f
	8	+
	9	+real(8) :: eps1
	10	+eps1=1.d-10
	11	+call fvn_muller(f,1.d-12,1.d-10,0,0,10,roots,200,infer,info)
	12	+write(,) "Error code :",info
	13	+do i=1,10
	14	+ write(,) roots(i),infer(i)
	15	+enddo
	16	+end program
	17	+
	18	+function f(x)
	19	+ complex(8) :: x,f
	20	+ f=x**10-1
	21	+end function

fvn_test/test_sparse.f90

Diff comments View file @ 27d3b84

	1	+program test_sparse
	2	+use fvn
	3	+implicit none
	4	+integer(4), parameter :: nz=12
	5	+integer(4), parameter :: n=5
	6	+real(8),dimension(nz) :: A
	7	+real(8),dimension(n,n) :: As
	8	+integer(4),dimension(nz) :: Ti,Tj
	9	+real(8),dimension(n) :: B,x
	10	+integer(4) :: status,i
	11	+! Description of the matrix in triplet form
	12	+A = (/ 2.,3.,3.,-1.,4.,4.,-3.,1.,2.,2.,6.,1. /)
	13	+B = (/ 8., 45., -3., 3., 19./)
	14	+Ti = (/ 1,2,1,3,5,2,3,4,5,3,2,5 /)
	15	+Tj = (/ 1,1,2,2,2,3,3,3,3,4,5,5 /)
	16	+
	17	+! Reconstruction of the matrix in standard form
	18	+! just needed for printing the matrix here
	19	+As=0.
	20	+do i=1,nz
	21	+ As(Ti(i),Tj(i))=A(i)
	22	+end do
	23	+write(,) "Matrix in standard representation :"
	24	+do i=1,5
	25	+ write(*,'(5f8.4)') As(i,:)
	26	+end do
	27	+write(,)
	28	+write(*,'("Right hand side :",5f8.4)') B
	29	+
	30	+!specific routine that will be used here
	31	+!call fvn_di_sparse_solve(n,nz,A,Ti,Tj,B,x,status)
	32	+call fvn_sparse_solve(n,nz,A,Ti,Tj,B,x,status)
	33	+write(*,'("Solution :",5f8.4)') x
	34	+write(*,'("Product matrix Solution :",5f8.4)') matmul(As,x)
	35	+end program

fvn_test/test_specfunc.f90

Diff comments View file @ 27d3b84

	1	+program test_specfunc
	2	+use fvn_fnlib
	3	+implicit none
	4	+integer,parameter :: npoints=200
	5	+integer :: i
	6	+real(8), dimension(npoints) :: j0
	7	+real(8) :: xmin,xmax,xstep,x,y
	8	+
	9	+! bsj0
	10	+xmin=-20.
	11	+xmax=20.
	12	+xstep=(xmax-xmin)/dble(npoints)
	13	+open(2,file='bsj0.dat')
	14	+do i=1,npoints
	15	+x=xmin+i*xstep
	16	+y=bsj0(x)
	17	+write(2,'(2e22.14)') x,y
	18	+end do
	19	+close(2)
	20	+
	21	+! bsj1
	22	+xmin=-20.
	23	+xmax=20.
	24	+xstep=(xmax-xmin)/dble(npoints)
	25	+open(2,file='bsj1.dat')
	26	+do i=1,npoints
	27	+x=xmin+i*xstep
	28	+y=bsj1(x)
	29	+write(2,'(2e22.14)') x,y
	30	+end do
	31	+close(2)
	32	+
	33	+!bsi0
	34	+xmin=-4.
	35	+xmax=4.
	36	+xstep=(xmax-xmin)/dble(npoints)
	37	+open(2,file='bsi0.dat')
	38	+do i=1,npoints
	39	+x=xmin+i*xstep
	40	+y=bsi0(x)
	41	+write(2,'(2e22.14)') x,y
	42	+end do
	43	+close(2)
	44	+
	45	+!bsi1
	46	+xmin=-4.
	47	+xmax=4.
	48	+xstep=(xmax-xmin)/dble(npoints)
	49	+open(2,file='bsi1.dat')
	50	+do i=1,npoints
	51	+x=xmin+i*xstep
	52	+y=bsi1(x)
	53	+write(2,'(2e22.14)') x,y
	54	+end do
	55	+close(2)
	56	+
	57	+!bsy0
	58	+xmin=0.
	59	+xmax=20.
	60	+xstep=(xmax-xmin)/dble(npoints)
	61	+open(2,file='bsy0.dat')
	62	+do i=1,npoints
	63	+x=xmin+i*xstep
	64	+y=bsy0(x)
	65	+write(2,'(2e22.14)') x,y
	66	+end do
	67	+close(2)
	68	+
	69	+!bsy1
	70	+xmin=0.
	71	+xmax=20.
	72	+xstep=(xmax-xmin)/dble(npoints)
	73	+open(2,file='bsy1.dat')
	74	+do i=1,npoints
	75	+x=xmin+i*xstep
	76	+y=bsy1(x)
	77	+write(2,'(2e22.14)') x,y
	78	+end do
	79	+close(2)
	80	+
	81	+!bsk0
	82	+xmin=0.
	83	+xmax=4.
	84	+xstep=(xmax-xmin)/dble(npoints)
	85	+open(2,file='bsk0.dat')
	86	+do i=1,npoints
	87	+x=xmin+i*xstep
	88	+y=bsk0(x)
	89	+write(2,'(2e22.14)') x,y
	90	+end do
	91	+close(2)
	92	+
	93	+!bsk1
	94	+xmin=0.
	95	+xmax=4.
	96	+xstep=(xmax-xmin)/dble(npoints)
	97	+open(2,file='bsk1.dat')
	98	+do i=1,npoints
	99	+x=xmin+i*xstep
	100	+y=bsk1(x)
	101	+write(2,'(2e22.14)') x,y
	102	+end do
	103	+close(2)
	104	+
	105	+!erf
	106	+xmin=-4.
	107	+xmax=4.
	108	+xstep=(xmax-xmin)/dble(npoints)
	109	+open(2,file='erf.dat')
	110	+do i=1,npoints
	111	+x=xmin+i*xstep
	112	+y=erf(x)
	113	+write(2,'(2e22.14)') x,y
	114	+end do
	115	+close(2)
	116	+
	117	+! gamma
	118	+xmin=-3.
	119	+xmax=7.
	120	+xstep=(xmax-xmin)/2000.
	121	+open(2,file='gamma.dat')
	122	+do i=1,2000
	123	+x=xmin+i*xstep
	124	+if ((abs(x-nint(x)) >= 1d-6) .or. x>0. )then
	125	+ y=gamma(x)
	126	+ write(2,'(2e22.14)') x,y
	127	+end if
	128	+end do
	129	+close(2)
	130	+
	131	+! 1/gamma
	132	+xmin=-3.
	133	+xmax=7.
	134	+xstep=(xmax-xmin)/dble(npoints)
	135	+open(2,file='gamr.dat')
	136	+do i=1,npoints
	137	+x=xmin+i*xstep
	138	+y=gamr(x)
	139	+write(2,'(2e22.14)') x,y
	140	+end do
	141	+close(2)
	142	+
	143	+! ei
	144	+xmin=0.
	145	+xmax=1.
	146	+xstep=(xmax-xmin)/dble(npoints)
	147	+open(2,file='ei.dat')
	148	+do i=1,npoints
	149	+x=xmin+i*xstep
	150	+y=ei(x)
	151	+write(2,'(2e22.14)') x,y
	152	+end do
	153	+close(2)
	154	+
	155	+! e1
	156	+xmin=0.
	157	+xmax=1.
	158	+xstep=(xmax-xmin)/dble(npoints)
	159	+open(2,file='e1.dat')
	160	+do i=1,npoints
	161	+x=xmin+i*xstep
	162	+y=e1(x)
	163	+write(2,'(2e22.14)') x,y
	164	+end do
	165	+close(2)
	166	+
	167	+
	168	+
	169	+end program test_specfunc

GITLAB

git-svn-id: https://lxsd.femto-st.fr/svn/fvn@36 b657c933-2333-4658-acf2-d3c7c2708721