diff --git a/fvn_test/Makefile b/fvn_test/Makefile
new file mode 100644
index 0000000..f3cb1f5
--- /dev/null
+++ b/fvn_test/Makefile
@@ -0,0 +1,21 @@
+
+include $(BTREE)/Make.inc
+
+programs = test_fac test_matinv test_specfunc  \
+test_det test_matcon test_matev test_sparse test_inter1d \
+test_inter2d test_inter3d test_akima test_lsp test_muller \
+test_integ
+
+prog:$(programs)
+
+clean:
+		rm -f {*.o,*.oo,*.ipo,*.a,*.mod}
+		rm -f $(programs)
+
+%: %.o
+		$(LINK) $(LINKFLAGS) $< init_random_seed.o $(LINKFVN) -o $@
+
+%.o: %.f90
+		$(F95) $(F95FLAGS) -c $<
+
+$(programs): init_random_seed.o
diff --git a/fvn_test/init_random_seed.f90 b/fvn_test/init_random_seed.f90
new file mode 100644
index 0000000..4cc1d78
--- /dev/null
+++ b/fvn_test/init_random_seed.f90
@@ -0,0 +1,31 @@
+! As the norm stated for random_seed "If no argument is present, the processor assigns a 
+! processor-dependent value to the seed."
+! It seems that there are different possible interpretations for this :
+! In gfortran a call to random_seed without arguments will initialize to a default value
+! which will always be the same.
+! In Portland fortran, a call to random_seed without arguments will initialize with something 
+! machine and time dependant
+! Unfortunately as the underlying algorithm are (or at least can be) different in each
+! compiler, we cannot generalize the use of this initialization provided in gfortran 
+! documentation
+SUBROUTINE init_random_seed()
+    INTEGER :: i, n, clock
+    INTEGER, DIMENSION(:), ALLOCATABLE :: seed
+    logical :: gfortran=.false.
+    !
+    ! Set gfortran to .true. if using gfortran to obtain the desired behaviour
+    !
+    if (gfortran) then
+        CALL RANDOM_SEED(size = n)
+        ALLOCATE(seed(n))
+
+        CALL SYSTEM_CLOCK(COUNT=clock)
+
+        seed = clock + 37 * (/ (i - 1, i = 1, n) /)
+        CALL RANDOM_SEED(PUT = seed)
+
+        DEALLOCATE(seed)
+    else
+        call random_seed()
+    end if
+END SUBROUTINE
diff --git a/fvn_test/test_akima.f90 b/fvn_test/test_akima.f90
new file mode 100644
index 0000000..31088aa
--- /dev/null
+++ b/fvn_test/test_akima.f90
@@ -0,0 +1,37 @@
+program akima
+ use fvn
+ implicit none
+ integer :: nbpoints,nppoints,i
+ real(8),dimension(:),allocatable :: x_d,y_d,breakpoints_d
+ real(8),dimension(:,:),allocatable :: coeff_fvn_d
+ real(8) :: xstep_d,xp_d,ty_d,fvn_y_d
+ open(2,file='fvn_akima_double.dat')
+ open(3,file='fvn_akima_breakpoints_double.dat')
+ nbpoints=30
+ allocate(x_d(nbpoints))
+ allocate(y_d(nbpoints))
+ allocate(breakpoints_d(nbpoints))
+ allocate(coeff_fvn_d(4,nbpoints))
+ xstep_d=20./dfloat(nbpoints)
+ do i=1,nbpoints
+    x_d(i)=-10.+dfloat(i)*xstep_d
+    y_d(i)=dsin(x_d(i))
+    write(3,44) x_d(i),y_d(i)
+ end do
+ close(3)
+ call fvn_akima(nbpoints,x_d,y_d,breakpoints_d,coeff_fvn_d)
+ nppoints=1000
+ xstep_d=22./dfloat(nppoints)
+ do i=1,nppoints
+    xp_d=-11.+dfloat(i)*xstep_d
+    ty_d=dsin(xp_d)
+    fvn_y_d=fvn_spline_eval(xp_d,nbpoints-1,breakpoints_d,coeff_fvn_d)
+    write(2,44) xp_d,ty_d,fvn_y_d
+ end do
+ close(2)
+deallocate(coeff_fvn_d,breakpoints_d,y_d,x_d)
+write(*,*) "All done, plot results with gnuplot using command :"
+write(*,*) "pl 'fvn_akima_double.dat' u 1:2 w l,'fvn_akima_breakpoints_double.dat' w p"
+
+44       FORMAT(4(1X,1PE22.14))
+end program
diff --git a/fvn_test/test_det.f90 b/fvn_test/test_det.f90
new file mode 100644
index 0000000..376c520
--- /dev/null
+++ b/fvn_test/test_det.f90
@@ -0,0 +1,17 @@
+program test_det
+use fvn
+implicit none
+real(8),dimension(3,3) :: a
+real(8) :: deta
+integer :: status,i
+
+call init_random_seed()
+call random_number(a)
+a=a*100
+deta=fvn_det(3,a,status)
+do i=1,3
+    write (*,'(3(f15.5))') a(i,:)
+end do
+write (*,*)
+write (*,*) "Det = ",deta,"Status = ",status
+end program
diff --git a/fvn_test/test_integ.f90 b/fvn_test/test_integ.f90
new file mode 100644
index 0000000..383feda
--- /dev/null
+++ b/fvn_test/test_integ.f90
@@ -0,0 +1,39 @@
+program integ
+ use fvn
+ implicit none
+ real(8), external :: f1,f2,g,h
+ real(8) :: a,b,epsabs,epsrel,abserr,res
+ integer :: key,ier
+ a=0.
+ b=1.
+ epsabs=1d-8
+ epsrel=1d-8
+ key=2
+ call fvn_integ_1_gk(f1,a,b,epsabs,epsrel,key,res,abserr,ier)
+ write(*,*) "Integration of x*x between 0 and 1 : "
+ write(*,*) res
+ call fvn_integ_2_gk(f2,a,b,g,h,epsabs,epsrel,key,res,abserr,ier)
+ write(*,*) "Integration of x*y between 0 and 1 on both x and y : "
+ write(*,*) res
+
+end program
+function f1(x)
+ implicit none
+    real(8) :: x,f1
+    f1=x*x
+end function
+function f2(x,y)
+ implicit none
+    real(8) :: x,y,f2
+    f2=x*y
+end function
+function g(x)
+ implicit none
+    real(8) :: x,g
+    g=0.
+end function
+function h(x)
+ implicit none
+    real(8) :: x,h
+    h=1.
+end function
diff --git a/fvn_test/test_inter1d.f90 b/fvn_test/test_inter1d.f90
new file mode 100644
index 0000000..a1f0a08
--- /dev/null
+++ b/fvn_test/test_inter1d.f90
@@ -0,0 +1,25 @@
+program inter1d
+use fvn
+implicit none
+integer(kind=4),parameter :: ndata=33
+integer(kind=4) :: i,nout
+real(kind=8) :: f,fdata(ndata),h,pi,q,sin,x,xdata(ndata)
+real(kind=8) ::tv
+intrinsic sin
+f(x)=sin(x)
+xdata(1)=0.
+fdata(1)=f(xdata(1))
+h=1./32.
+do i=2,ndata
+xdata(i)=xdata(i-1)+h
+fdata(i)=f(xdata(i))
+end do
+call init_random_seed()
+call random_number(x)
+q=fvn_quad_interpol(x,ndata,xdata,fdata)
+tv=f(x)
+write(*,'("x y z ",1(f8.5))') x
+write(*,'("Calculated (real) value :",f8.5)') tv
+write(*,'("fvn interpolation : ",f8.5)') q
+write(*,'("Relative fvn error :",e12.5)') abs((q-tv)/tv)
+end program
\ No newline at end of file
diff --git a/fvn_test/test_inter2d.f90 b/fvn_test/test_inter2d.f90
new file mode 100644
index 0000000..a8bc359
--- /dev/null
+++ b/fvn_test/test_inter2d.f90
@@ -0,0 +1,30 @@
+program inter2d
+use fvn
+implicit none
+integer(kind=4),parameter :: nx=21,ny=42
+integer(kind=4) :: i,j
+real(kind=8) :: f,fdata(nx,ny),dble,pi,q,sin,x,xdata(nx),y,ydata(ny)
+real(kind=8) :: tv
+intrinsic dble,sin
+f(x,y)=sin(x+2.*y)
+do i=1,nx
+xdata(i)=dble(i-1)/dble(nx-1)
+end do
+do i=1,ny
+ydata(i)=dble(i-1)/dble(ny-1)
+end do
+do i=1,nx
+do j=1,ny
+fdata(i,j)=f(xdata(i),ydata(j))
+end do
+end do
+call init_random_seed()
+call random_number(x)
+call random_number(y)
+q=fvn_quad_2d_interpol(x,y,nx,xdata,ny,ydata,fdata)
+tv=f(x,y)
+write(*,'("x y z ",2(f8.5))') x,y
+write(*,'("Calculated (real) value :",f8.5)') tv
+write(*,'("fvn interpolation : ",f8.5)') q
+write(*,'("Relative fvn error :",e12.5)') abs((q-tv)/tv)
+end program
\ No newline at end of file
diff --git a/fvn_test/test_inter3d.f90 b/fvn_test/test_inter3d.f90
new file mode 100644
index 0000000..e00ade6
--- /dev/null
+++ b/fvn_test/test_inter3d.f90
@@ -0,0 +1,36 @@
+program test_inter3d
+use fvn
+implicit none
+integer(kind=4),parameter :: nx=21,ny=42,nz=18
+integer(kind=4) :: i,j,k
+real(kind=8) :: f,fdata(nx,ny,nz),dble,pi,q,sin,x,xdata(nx),y,ydata(ny),z,zdata(nz)
+real(kind=8) :: tv
+intrinsic dble,sin
+f(x,y,z)=sin(x+2.*y+3.*z)
+do i=1,nx
+      xdata(i)=2.*(dble(i-1)/dble(nx-1))
+end do
+do i=1,ny
+      ydata(i)=2.*(dble(i-1)/dble(ny-1))
+end do
+do i=1,nz
+      zdata(i)=2.*(dble(i-1)/dble(nz-1))
+end do
+do i=1,nx
+      do j=1,ny
+            do k=1,nz
+                  fdata(i,j,k)=f(xdata(i),ydata(j),zdata(k))
+            end do
+      end do
+end do
+call init_random_seed()
+call random_number(x)
+call random_number(y)
+call random_number(z)
+q=fvn_quad_3d_interpol(x,y,z,nx,xdata,ny,ydata,nz,zdata,fdata)
+tv=f(x,y,z)
+write(*,'("x y z ",3(f8.5))') x,y,z
+write(*,'("Calculated (real) value :",f8.5)')  tv
+write(*,'("fvn interpolation : ",f8.5)') q
+write(*,'("Relative fvn error :",e12.5)') abs((q-tv)/tv)
+end program
diff --git a/fvn_test/test_lsp.f90 b/fvn_test/test_lsp.f90
new file mode 100644
index 0000000..42922d7
--- /dev/null
+++ b/fvn_test/test_lsp.f90
@@ -0,0 +1,39 @@
+program lsp
+use fvn
+implicit none
+integer,parameter :: npoints=13,deg=3
+ integer :: status,i
+ real(kind=8) :: xm(npoints),ym(npoints),xstep,xc,yc
+ real(kind=8) :: coeff(deg+1)
+ 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. /)
+ 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 /)
+ open(2,file='fvn_lsp_double_mesure.dat')
+ open(3,file='fvn_lsp_double_poly.dat')
+ do i=1,npoints
+    write(2,44) xm(i),ym(i)
+ end do
+ close(2)
+ call fvn_lspoly(npoints,xm,ym,deg,coeff,status)
+ xstep=(xm(npoints)-xm(1))/1000.
+ do i=1,1000
+    xc=xm(1)+(i-1)*xstep
+    yc=poly(xc,coeff)
+    write(3,44) xc,yc
+ end do
+ close(3)
+write(*,*) "All done, plot results with gnuplot using command :"
+write(*,*) "pl 'fvn_lsp_double_mesure.dat' u 1:2 w p,'fvn_lsp_double_poly.dat' u 1:2 w l"
+44      FORMAT(4(1X,1PE22.14))
+contains
+function poly(x,coeff)
+    implicit none
+    real(8) :: x
+    real(8) :: coeff(deg+1)
+    real(8) :: poly
+    integer :: i
+    poly=0.
+    do i=1,deg+1
+         poly=poly+coeff(i)*x**(i-1)
+    end do
+end function
+end program
diff --git a/fvn_test/test_matcon.f90 b/fvn_test/test_matcon.f90
new file mode 100644
index 0000000..cc956fb
--- /dev/null
+++ b/fvn_test/test_matcon.f90
@@ -0,0 +1,28 @@
+program test_matcon
+use fvn
+implicit none
+real(8),dimension(3,3) :: a
+real(8) :: rcond
+integer :: status,i
+call init_random_seed()
+call random_number(a)
+a=a*100
+call fvn_matcon(3,a,rcond,status)
+write(*,*) "Reasonnably conditionned matrix"
+do i=1,3
+    write (*,'(3(e12.5))') a(i,:)
+end do
+write (*,*)
+write (*,*) "Cond = ",rcond
+write (*,*)
+write (*,*)
+a(1,1)=a(1,1)*1d9
+write(*,*) "Badly conditionned matrix"
+do i=1,3
+    write (*,'(3(e12.5))') a(i,:)
+end do
+call fvn_matcon(3,a,rcond,status)
+write (*,*)
+write (*,*) "Cond = ",rcond
+
+end program
\ No newline at end of file
diff --git a/fvn_test/test_matev.f90 b/fvn_test/test_matev.f90
new file mode 100644
index 0000000..982a1e4
--- /dev/null
+++ b/fvn_test/test_matev.f90
@@ -0,0 +1,24 @@
+program test_matev
+use fvn
+implicit none
+real(8),dimension(3,3) :: a
+complex(8),dimension(3) :: evala
+complex(8),dimension(3,3) :: eveca
+integer :: status,i,j
+call init_random_seed()
+call random_number(a)
+a=a*100
+call fvn_matev(3,a,evala,eveca,status)
+
+write(*,*) "The matrix :"
+write (*,'(3(e12.5))') a
+write (*,*)
+do i=1,3
+write(*,'("Eigenvalue ",I3," : (",e12.5,",",e12.5,")  ")') i,evala(i)
+write(*,*) "Associated Eigenvector :"
+do j=1,3
+write(*,'("(",e12.5,",",e12.5,")  ")') eveca(j,i)
+end do
+write(*,*)
+end do
+end program
\ No newline at end of file
diff --git a/fvn_test/test_matinv.f90 b/fvn_test/test_matinv.f90
new file mode 100644
index 0000000..e493a53
--- /dev/null
+++ b/fvn_test/test_matinv.f90
@@ -0,0 +1,38 @@
+program test_matinv
+use fvn
+implicit none
+integer(4), parameter :: n=3
+integer(4) :: status,i
+complex(8),dimension(n,n) :: m,im,prod
+real(8),dimension(n,n) :: rtmp,itmp
+character(len=80) :: fmreal,fmcmplx
+
+fmcmplx='(3("(",f8.5,",",f8.5,")  "))'
+
+! initialize pseudo random generator
+call init_random_seed()
+! fill real and imaginary part
+call random_number(rtmp)
+call random_number(itmp)
+! create the complex matrix (fvn_i is defined in the fvn module)
+m=rtmp+fvn_i*itmp
+write(*,*) "Matrix M"
+do i=1,n
+    write(*,fmcmplx) m(i,:)
+end do
+
+! Invertion
+call fvn_matinv(n,m,im)
+write(*,*) "Inverse of M"
+do i=1,n
+    write(*,fmcmplx) im(i,:)
+end do
+
+! Result should be identity matrix
+write(*,*) "Product of M and inverse of M :"
+prod=matmul(m,im)
+do i=1,n
+    write (*,fmcmplx) prod(i,:)
+end do
+
+end program
diff --git a/fvn_test/test_muller.f90 b/fvn_test/test_muller.f90
new file mode 100644
index 0000000..f0bfcc4
--- /dev/null
+++ b/fvn_test/test_muller.f90
@@ -0,0 +1,21 @@
+program muller
+use fvn
+implicit none
+integer :: i,info
+complex(8),dimension(10) :: roots
+integer,dimension(10) :: infer
+complex(8), external :: f
+
+real(8) :: eps1
+eps1=1.d-10
+call fvn_muller(f,1.d-12,1.d-10,0,0,10,roots,200,infer,info)
+write(*,*) "Error code :",info
+do i=1,10
+   write(*,*) roots(i),infer(i)
+enddo
+end program
+
+function f(x)
+   complex(8) :: x,f
+   f=x**10-1
+end function
diff --git a/fvn_test/test_sparse.f90 b/fvn_test/test_sparse.f90
new file mode 100644
index 0000000..7eb54b3
--- /dev/null
+++ b/fvn_test/test_sparse.f90
@@ -0,0 +1,35 @@
+program test_sparse
+use fvn
+implicit none
+integer(4), parameter :: nz=12
+integer(4), parameter :: n=5
+real(8),dimension(nz) :: A
+real(8),dimension(n,n) :: As
+integer(4),dimension(nz) :: Ti,Tj
+real(8),dimension(n) :: B,x
+integer(4) :: status,i
+! Description of the matrix in triplet form
+A = (/ 2.,3.,3.,-1.,4.,4.,-3.,1.,2.,2.,6.,1. /)
+B = (/ 8., 45., -3., 3., 19./)
+Ti = (/ 1,2,1,3,5,2,3,4,5,3,2,5 /)
+Tj = (/ 1,1,2,2,2,3,3,3,3,4,5,5 /)
+
+! Reconstruction of the matrix in standard form
+! just needed for printing the matrix here
+As=0.
+do i=1,nz
+    As(Ti(i),Tj(i))=A(i)
+end do
+write(*,*) "Matrix in standard representation :"
+do i=1,5
+    write(*,'(5f8.4)') As(i,:)
+end do
+write(*,*)
+write(*,'("Right hand side :",5f8.4)') B
+
+!specific routine that will be used here
+!call fvn_di_sparse_solve(n,nz,A,Ti,Tj,B,x,status)
+call fvn_sparse_solve(n,nz,A,Ti,Tj,B,x,status)
+write(*,'("Solution :",5f8.4)') x
+write(*,'("Product matrix Solution :",5f8.4)') matmul(As,x)
+end program
\ No newline at end of file
diff --git a/fvn_test/test_specfunc.f90 b/fvn_test/test_specfunc.f90
new file mode 100644
index 0000000..cb9fc10
--- /dev/null
+++ b/fvn_test/test_specfunc.f90
@@ -0,0 +1,170 @@
+program test_specfunc
+use fvn_fnlib
+implicit none
+integer,parameter :: npoints=200
+integer :: i
+real(8), dimension(npoints) :: j0
+real(8) :: xmin,xmax,xstep,x,y
+
+! bsj0
+xmin=-20.
+xmax=20.
+xstep=(xmax-xmin)/dble(npoints)
+open(2,file='bsj0.dat')
+do i=1,npoints
+x=xmin+i*xstep
+y=bsj0(x)
+write(2,'(2e22.14)') x,y
+end do
+close(2)
+
+! bsj1
+xmin=-20.
+xmax=20.
+xstep=(xmax-xmin)/dble(npoints)
+open(2,file='bsj1.dat')
+do i=1,npoints
+x=xmin+i*xstep
+y=bsj1(x)
+write(2,'(2e22.14)') x,y
+end do
+close(2)
+
+!bsi0
+xmin=-4.
+xmax=4.
+xstep=(xmax-xmin)/dble(npoints)
+open(2,file='bsi0.dat')
+do i=1,npoints
+x=xmin+i*xstep
+y=bsi0(x)
+write(2,'(2e22.14)') x,y
+end do
+close(2)
+
+!bsi1
+xmin=-4.
+xmax=4.
+xstep=(xmax-xmin)/dble(npoints)
+open(2,file='bsi1.dat')
+do i=1,npoints
+x=xmin+i*xstep
+y=bsi1(x)
+write(2,'(2e22.14)') x,y
+end do
+close(2)
+
+!bsy0
+xmin=0.
+xmax=20.
+xstep=(xmax-xmin)/dble(npoints)
+open(2,file='bsy0.dat')
+do i=1,npoints
+x=xmin+i*xstep
+y=bsy0(x)
+write(2,'(2e22.14)') x,y
+end do
+close(2)
+
+!bsy1
+xmin=0.
+xmax=20.
+xstep=(xmax-xmin)/dble(npoints)
+open(2,file='bsy1.dat')
+do i=1,npoints
+x=xmin+i*xstep
+y=bsy1(x)
+write(2,'(2e22.14)') x,y
+end do
+close(2)
+
+!bsk0
+xmin=0.
+xmax=4.
+xstep=(xmax-xmin)/dble(npoints)
+open(2,file='bsk0.dat')
+do i=1,npoints
+x=xmin+i*xstep
+y=bsk0(x)
+write(2,'(2e22.14)') x,y
+end do
+close(2)
+
+!bsk1
+xmin=0.
+xmax=4.
+xstep=(xmax-xmin)/dble(npoints)
+open(2,file='bsk1.dat')
+do i=1,npoints
+x=xmin+i*xstep
+y=bsk1(x)
+write(2,'(2e22.14)') x,y
+end do
+close(2)
+
+!erf
+xmin=-4.
+xmax=4.
+xstep=(xmax-xmin)/dble(npoints)
+open(2,file='erf.dat')
+do i=1,npoints
+x=xmin+i*xstep
+y=erf(x)
+write(2,'(2e22.14)') x,y
+end do
+close(2)
+
+! gamma
+xmin=-3.
+xmax=7.
+xstep=(xmax-xmin)/2000.
+open(2,file='gamma.dat')
+do i=1,2000
+x=xmin+i*xstep
+if ((abs(x-nint(x)) >= 1d-6) .or. x>0. )then
+    y=gamma(x)
+    write(2,'(2e22.14)') x,y
+end if
+end do
+close(2)
+
+! 1/gamma
+xmin=-3.
+xmax=7.
+xstep=(xmax-xmin)/dble(npoints)
+open(2,file='gamr.dat')
+do i=1,npoints
+x=xmin+i*xstep
+y=gamr(x)
+write(2,'(2e22.14)') x,y
+end do
+close(2)
+
+! ei
+xmin=0.
+xmax=1.
+xstep=(xmax-xmin)/dble(npoints)
+open(2,file='ei.dat')
+do i=1,npoints
+x=xmin+i*xstep
+y=ei(x)
+write(2,'(2e22.14)') x,y
+end do
+close(2)
+
+! e1
+xmin=0.
+xmax=1.
+xstep=(xmax-xmin)/dble(npoints)
+open(2,file='e1.dat')
+do i=1,npoints
+x=xmin+i*xstep
+y=e1(x)
+write(2,'(2e22.14)') x,y
+end do
+close(2)
+
+
+
+end program test_specfunc
+