module fvn
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! fvn : a f95 module replacement for some imsl routines
! it uses lapack for linear algebra
! it uses modified quadpack for integration
!
! William Daniau 2007->today
! william.daniau@femto-st.fr
!
! Routines naming scheme :
!
!           fvn_x_name
!           where x can be  s : real 
!                           d : real double precision
!                           c : complex
!                           z : double complex
!
!
! This piece of code is totally free! Do whatever you want with it. However
! if you find it usefull it would be kind to give credits ;-) 
!
! svn version
! February 2008 : added fnlib to repository so there's no use to have
!                 special functions and trigonometry here. Some functions
!                 are then removed
! January 2008 : added quadratic interpolation, gamma/factorial function,
!                 a function which return identity matrix,
!                 evaluation of nterm chebyshev series
! September 2007 : added sparse system solving by interfacing umfpack
! June 2007 : added some complex trigonometric functions
!
! TO DO LIST :
! + Order eigenvalues and vectors in decreasing eigenvalue's modulus order -> atm 
! eigenvalues are given with no particular order.
! + Generic interface for fvn_x_name family  -> fvn_name
! + Make some parameters optional, status for example
! + use f95 kinds "double complex" -> complex(kind=8)
! + unify quadpack routines
! + ...
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

implicit none
! We define pi and i for the module
real(kind=8),parameter :: fvn_pi = 3.141592653589793_8
complex(kind=8),parameter :: fvn_i = (0._8,1._8)


! All quadpack routines are private to the module
private ::  d1mach,dqag,dqag_2d_inner,dqag_2d_outer,dqage,dqage_2d_inner, &
            dqage_2d_outer,dqk15,dqk15_2d_inner,dqk15_2d_outer,dqk21,dqk21_2d_inner,dqk21_2d_outer, &
            dqk31,dqk31_2d_inner,dqk31_2d_outer,dqk41,dqk41_2d_inner,dqk41_2d_outer, &
            dqk51,dqk51_2d_inner,dqk51_2d_outer,dqk61,dqk61_2d_inner,dqk61_2d_outer,dqpsrt

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Generic interface Definition
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

! Identity Matrix
interface fvn_ident
      module procedure fvn_s_ident,fvn_d_ident,fvn_c_ident,fvn_z_ident
end interface fvn_ident

! Matrix inversion
interface fvn_matinv
      module procedure fvn_s_matinv,fvn_d_matinv,fvn_c_matinv,fvn_z_matinv
end interface fvn_matinv

! Determinant
interface fvn_det
    module procedure fvn_s_det,fvn_d_det,fvn_c_det,fvn_z_det
end interface fvn_det

! Condition
interface fvn_matcon
    module procedure fvn_s_matcon,fvn_d_matcon,fvn_c_matcon,fvn_z_matcon
end interface fvn_matcon

! Eigen
interface fvn_matev
    module procedure fvn_s_matev,fvn_d_matev,fvn_c_matev,fvn_z_matev
end interface fvn_matev

! Utility procedure find interval
interface fvn_find_interval
    module procedure fvn_s_find_interval,fvn_d_find_interval
end interface fvn_find_interval

! Quadratic 1D interpolation
interface fvn_quad_interpol
    module procedure fvn_s_quad_interpol,fvn_d_quad_interpol
end interface fvn_quad_interpol

! Quadratic 2D interpolation
interface fvn_quad_2d_interpol
    module procedure fvn_s_quad_2d_interpol,fvn_d_quad_2d_interpol
end interface fvn_quad_2d_interpol

! Quadratic 3D interpolation
interface fvn_quad_3d_interpol
    module procedure fvn_s_quad_3d_interpol,fvn_d_quad_3d_interpol
end interface fvn_quad_3d_interpol

! Akima interpolation
interface fvn_akima
    module procedure fvn_s_akima,fvn_d_akima
end interface fvn_akima

! Akima evaluation
interface fvn_spline_eval
    module procedure fvn_s_spline_eval,fvn_d_spline_eval
end interface fvn_spline_eval

! Least square polynomial
interface fvn_lspoly
    module procedure fvn_s_lspoly,fvn_d_lspoly
end interface fvn_lspoly

! Muller
interface fvn_muller
    module procedure fvn_z_muller
end interface fvn_muller

! Gauss legendre
interface fvn_gauss_legendre
    module procedure fvn_d_gauss_legendre
end interface fvn_gauss_legendre

! Simple Gauss Legendre integration
interface fvn_gl_integ
    module procedure fvn_d_gl_integ
end interface fvn_gl_integ

! Adaptative Gauss Kronrod integration f(x)
interface fvn_integ_1_gk
    module procedure fvn_d_integ_1_gk
end interface fvn_integ_1_gk

! Adaptative Gauss Kronrod integration f(x,y)
interface fvn_integ_2_gk
    module procedure fvn_d_integ_2_gk
end interface fvn_integ_2_gk

! Sparse solving
interface fvn_sparse_solve
    module procedure fvn_zl_sparse_solve,fvn_zi_sparse_solve,fvn_dl_sparse_solve,fvn_di_sparse_solve
end interface fvn_sparse_solve


contains 

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! Identity Matrix
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
function fvn_d_ident(n)
      implicit none
      integer(kind=4) :: n
      real(kind=8), dimension(n,n) :: fvn_d_ident

      real(kind=8),dimension(n*n) :: vect
      integer(kind=4) :: i

      vect=0._8
      vect(1:n*n:n+1) = 1._8
      fvn_d_ident=reshape(vect, shape = (/ n,n /))
end function

function fvn_s_ident(n)
      implicit none
      integer(kind=4) :: n
      real(kind=4), dimension(n,n) :: fvn_s_ident

      real(kind=4),dimension(n*n) :: vect
      integer(kind=4) :: i

     vect=0._4
      vect(1:n*n:n+1) = 1._4
      fvn_s_ident=reshape(vect, shape = (/ n,n /))
end function

function fvn_c_ident(n)
      implicit none
      integer(kind=4) :: n
      complex(kind=4), dimension(n,n) :: fvn_c_ident

      complex(kind=4),dimension(n*n) :: vect
      integer(kind=4) :: i

      vect=(0._4,0._4)
      vect(1:n*n:n+1) = (1._4,0._4)
      fvn_c_ident=reshape(vect, shape = (/ n,n /))
end function

function fvn_z_ident(n)
      implicit none
      integer(kind=4) :: n
      complex(kind=8), dimension(n,n) :: fvn_z_ident

      complex(kind=8),dimension(n*n) :: vect
      integer(kind=4) :: i

      vect=(0._8,0._8)
      vect(1:n*n:n+1) = (1._8,0._8)
      fvn_z_ident=reshape(vect, shape = (/ n,n /))
end function



!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! Matrix inversion subroutines
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

subroutine fvn_s_matinv(d,a,inva,status)
    !
    ! Matrix inversion of a real matrix using BLAS and LAPACK
    !
    ! d (in) : matrix rank
    ! a (in) : input matrix
    ! inva (out) : inversed matrix
    ! status (ou) : =0 if something failed
    !
    implicit none
    integer, intent(in) :: d
    real, intent(in) :: a(d,d)
    real, intent(out) :: inva(d,d)
    integer, intent(out),optional :: status

    integer, allocatable :: ipiv(:)
    real, allocatable :: work(:)
    real twork(1)
    integer :: info
    integer :: lwork

    if (present(status)) status=1

    allocate(ipiv(d))
    ! copy a into inva using BLAS
    !call scopy(d*d,a,1,inva,1)
    inva(:,:)=a(:,:)
    ! LU factorization using LAPACK
    call sgetrf(d,d,inva,d,ipiv,info)
    ! if info is not equal to 0, something went wrong we exit setting status to 0
    if (info /= 0) then
        if (present(status)) status=0
        deallocate(ipiv)
        return
    end if
    ! we use the query fonction of xxxtri to obtain the optimal workspace size
    call sgetri(d,inva,d,ipiv,twork,-1,info)
    lwork=int(twork(1))
    allocate(work(lwork))
    ! Matrix inversion using LAPACK
    call sgetri(d,inva,d,ipiv,work,lwork,info)
    ! again if info is not equal to 0, we exit setting status to 0
    if (info /= 0) then
        if (present(status)) status=0
    end if
    deallocate(work)
    deallocate(ipiv)
end subroutine

subroutine fvn_d_matinv(d,a,inva,status)
    !
    ! Matrix inversion of a double precision matrix using BLAS and LAPACK
    !
    ! d (in) : matrix rank
    ! a (in) : input matrix
    ! inva (out) : inversed matrix
    ! status (ou) : =0 if something failed
    !
    implicit none
    integer, intent(in), optional :: d
    double precision, intent(in) :: a(d,d)
    double precision, intent(out) :: inva(d,d)
    integer, intent(out),optional :: status

    integer, allocatable :: ipiv(:)
    double precision, allocatable :: work(:)
    double precision :: twork(1)
    integer :: info
    integer :: lwork

    if (present(status)) status=1

    allocate(ipiv(d))
    ! copy a into inva using BLAS
    !call dcopy(d*d,a,1,inva,1)
    inva(:,:)=a(:,:)
    ! LU factorization using LAPACK
    call dgetrf(d,d,inva,d,ipiv,info)
    ! if info is not equal to 0, something went wrong we exit setting status to 0
    if (info /= 0) then
        if (present(status)) status=0
        deallocate(ipiv)
        return
    end if
    ! we use the query fonction of xxxtri to obtain the optimal workspace size
    call dgetri(d,inva,d,ipiv,twork,-1,info)
    lwork=int(twork(1))
    allocate(work(lwork))
    ! Matrix inversion using LAPACK
    call dgetri(d,inva,d,ipiv,work,lwork,info)
    ! again if info is not equal to 0, we exit setting status to 0
    if (info /= 0) then
        if (present(status)) status=0
    end if
    deallocate(work)
    deallocate(ipiv)
end subroutine

subroutine fvn_c_matinv(d,a,inva,status)
    !
    ! Matrix inversion of a complex matrix using BLAS and LAPACK
    !
    ! d (in) : matrix rank
    ! a (in) : input matrix
    ! inva (out) : inversed matrix
    ! status (ou) : =0 if something failed
    !
    implicit none
    integer, intent(in) :: d
    complex, intent(in) :: a(d,d)
    complex, intent(out) :: inva(d,d)
    integer, intent(out),optional :: status

    integer, allocatable :: ipiv(:)
    complex, allocatable :: work(:)
    complex :: twork(1)
    integer :: info
    integer :: lwork

     if (present(status)) status=1

    allocate(ipiv(d))
    ! copy a into inva using BLAS
    !call ccopy(d*d,a,1,inva,1)
    inva(:,:)=a(:,:)
    
    ! LU factorization using LAPACK
    call cgetrf(d,d,inva,d,ipiv,info)
    ! if info is not equal to 0, something went wrong we exit setting status to 0
    if (info /= 0) then
        if (present(status)) status=0
        deallocate(ipiv)
        return
    end if
    ! we use the query fonction of xxxtri to obtain the optimal workspace size
    call cgetri(d,inva,d,ipiv,twork,-1,info)
    lwork=int(twork(1))
    allocate(work(lwork))
    ! Matrix inversion using LAPACK
    call cgetri(d,inva,d,ipiv,work,lwork,info)
    ! again if info is not equal to 0, we exit setting status to 0
    if (info /= 0) then
        if (present(status)) status=0
    end if
    deallocate(work)
    deallocate(ipiv)
end subroutine

subroutine fvn_z_matinv(d,a,inva,status)
    !
    ! Matrix inversion of a double complex matrix using BLAS and LAPACK
    !
    ! d (in) : matrix rank
    ! a (in) : input matrix
    ! inva (out) : inversed matrix
    ! status (ou) : =0 if something failed
    !
    implicit none
    integer, intent(in) :: d
    double complex, intent(in) :: a(d,d)
    double complex, intent(out) :: inva(d,d)
    integer, intent(out),optional :: status

    integer, allocatable :: ipiv(:)
    double complex, allocatable :: work(:)
    double complex :: twork(1)
    integer :: info
    integer :: lwork

     if (present(status)) status=1

    allocate(ipiv(d))
    ! copy a into inva using BLAS
    !call zcopy(d*d,a,1,inva,1)
    inva(:,:)=a(:,:)
    
    ! LU factorization using LAPACK
    call zgetrf(d,d,inva,d,ipiv,info)
    ! if info is not equal to 0, something went wrong we exit setting status to 0
    if (info /= 0) then
        if (present(status)) status=0
        deallocate(ipiv)
        return
    end if
    ! we use the query fonction of xxxtri to obtain the optimal workspace size
    call zgetri(d,inva,d,ipiv,twork,-1,info)
    lwork=int(twork(1))
    allocate(work(lwork))
    ! Matrix inversion using LAPACK
    call zgetri(d,inva,d,ipiv,work,lwork,info)
    ! again if info is not equal to 0, we exit setting status to 0
    if (info /= 0) then
        if (present(status)) status=0
    end if
    deallocate(work)
    deallocate(ipiv)
end subroutine


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! Determinants
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
function fvn_s_det(d,a,status)
    !
    ! Evaluate the determinant of a square matrix using lapack LU factorization
    !
    ! d (in) : matrix rank
    ! a (in) : The Matrix
    ! status (out) : =0 if LU factorization failed
    !
    implicit none
    integer, intent(in) :: d
    real, intent(in) :: a(d,d)
    integer, intent(out), optional :: status
    real :: fvn_s_det
    
    real, allocatable :: wc_a(:,:)
    integer, allocatable :: ipiv(:)
    integer :: info,i
    
     if (present(status)) status=1
    allocate(wc_a(d,d))
    allocate(ipiv(d))
    wc_a(:,:)=a(:,:)
    call sgetrf(d,d,wc_a,d,ipiv,info)
    if (info/= 0) then
         if (present(status)) status=0
        fvn_s_det=0.e0
        deallocate(ipiv)
        deallocate(wc_a)
        return
    end if
    fvn_s_det=1.e0
    do i=1,d
        if (ipiv(i)==i) then
            fvn_s_det=fvn_s_det*wc_a(i,i)
        else
            fvn_s_det=-fvn_s_det*wc_a(i,i)
        end if
    end do
    deallocate(ipiv)
    deallocate(wc_a)

end function

function fvn_d_det(d,a,status)
    !
    ! Evaluate the determinant of a square matrix using lapack LU factorization
    !
    ! d (in) : matrix rank
    ! a (in) : The Matrix
    ! status (out) : =0 if LU factorization failed
    !
    implicit none
    integer, intent(in) :: d
    double precision, intent(in) :: a(d,d)
    integer, intent(out), optional :: status
    double precision :: fvn_d_det
    
    double precision, allocatable :: wc_a(:,:)
    integer, allocatable :: ipiv(:)
    integer :: info,i
    
     if (present(status)) status=1
    allocate(wc_a(d,d))
    allocate(ipiv(d))
    wc_a(:,:)=a(:,:)
    call dgetrf(d,d,wc_a,d,ipiv,info)
    if (info/= 0) then
         if (present(status)) status=0
        fvn_d_det=0.d0
        deallocate(ipiv)
        deallocate(wc_a)
        return
    end if
    fvn_d_det=1.d0
    do i=1,d
        if (ipiv(i)==i) then
            fvn_d_det=fvn_d_det*wc_a(i,i)
        else
            fvn_d_det=-fvn_d_det*wc_a(i,i)
        end if
    end do
    deallocate(ipiv)
    deallocate(wc_a)

end function

function fvn_c_det(d,a,status)    !
    ! Evaluate the determinant of a square matrix using lapack LU factorization
    !
    ! d (in) : matrix rank
    ! a (in) : The Matrix
    ! status (out) : =0 if LU factorization failed
    !
    implicit none
    integer, intent(in) :: d
    complex, intent(in) :: a(d,d)
    integer, intent(out), optional :: status
    complex :: fvn_c_det
    
    complex, allocatable :: wc_a(:,:)
    integer, allocatable :: ipiv(:)
    integer :: info,i
    
     if (present(status)) status=1
    allocate(wc_a(d,d))
    allocate(ipiv(d))
    wc_a(:,:)=a(:,:)
    call cgetrf(d,d,wc_a,d,ipiv,info)
    if (info/= 0) then
         if (present(status)) status=0
        fvn_c_det=(0.e0,0.e0)
        deallocate(ipiv)
        deallocate(wc_a)
        return
    end if
    fvn_c_det=(1.e0,0.e0)
    do i=1,d
        if (ipiv(i)==i) then
            fvn_c_det=fvn_c_det*wc_a(i,i)
        else
            fvn_c_det=-fvn_c_det*wc_a(i,i)
        end if
    end do
    deallocate(ipiv)
    deallocate(wc_a)

end function

function fvn_z_det(d,a,status)
    !
    ! Evaluate the determinant of a square matrix using lapack LU factorization
    !
    ! d (in) : matrix rank
    ! a (in) : The Matrix
    ! det (out) : determinant
    ! status (out) : =0 if LU factorization failed
    !
    implicit none
    integer, intent(in) :: d
    double complex, intent(in) :: a(d,d)
    integer, intent(out), optional :: status
    double complex :: fvn_z_det
    
    double complex, allocatable :: wc_a(:,:)
    integer, allocatable :: ipiv(:)
    integer :: info,i
    
     if (present(status)) status=1
    allocate(wc_a(d,d))
    allocate(ipiv(d))
    wc_a(:,:)=a(:,:)
    call zgetrf(d,d,wc_a,d,ipiv,info)
    if (info/= 0) then
         if (present(status)) status=0
        fvn_z_det=(0.d0,0.d0)
        deallocate(ipiv)
        deallocate(wc_a)
        return
    end if
    fvn_z_det=(1.d0,0.d0)
    do i=1,d
        if (ipiv(i)==i) then
            fvn_z_det=fvn_z_det*wc_a(i,i)
        else
            fvn_z_det=-fvn_z_det*wc_a(i,i)
        end if
    end do
    deallocate(ipiv)
    deallocate(wc_a)

end function

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! Condition test
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! 1-norm 
! fonction lapack slange,dlange,clange,zlange pour obtenir la 1-norm
! fonction lapack sgecon,dgecon,cgecon,zgecon pour calculer la rcond
!
subroutine fvn_s_matcon(d,a,rcond,status)
    ! Matrix condition (reciprocal of condition number)
    ! 
    ! d (in) : matrix rank
    ! a (in) : The Matrix
    ! rcond (out) : guess what
    ! status (out) : =0 if something went wrong
    !
    implicit none
    integer, intent(in) :: d
    real, intent(in) :: a(d,d)
    real, intent(out) :: rcond
    integer, intent(out), optional :: status
    
    real, allocatable :: work(:)
    integer, allocatable :: iwork(:)
    real :: anorm
    real, allocatable :: wc_a(:,:) ! working copy of a
    integer :: info
    integer, allocatable :: ipiv(:)
    
    real, external :: slange
    
    
     if (present(status)) status=1
    
    anorm=slange('1',d,d,a,d,work) ! work is unallocated as it is only used when computing infinity norm
    
    allocate(wc_a(d,d))
    !call scopy(d*d,a,1,wc_a,1)
    wc_a(:,:)=a(:,:)
    
    allocate(ipiv(d))
    call sgetrf(d,d,wc_a,d,ipiv,info)
    if (info /= 0) then
         if (present(status)) status=0
        deallocate(ipiv)
        deallocate(wc_a)
        return
    end if
    allocate(work(4*d))
    allocate(iwork(d))
    call sgecon('1',d,wc_a,d,anorm,rcond,work,iwork,info)
    if (info /= 0) then
         if (present(status)) status=0
    end if
    deallocate(iwork)
    deallocate(work)
    deallocate(ipiv)
    deallocate(wc_a)

end subroutine

subroutine fvn_d_matcon(d,a,rcond,status)
    ! Matrix condition (reciprocal of condition number)
    ! 
    ! d (in) : matrix rank
    ! a (in) : The Matrix
    ! rcond (out) : guess what
    ! status (out) : =0 if something went wrong
    !
    implicit none
    integer, intent(in) :: d
    double precision, intent(in) :: a(d,d)
    double precision, intent(out) :: rcond
    integer, intent(out), optional :: status
    
    double precision, allocatable :: work(:)
    integer, allocatable :: iwork(:)
    double precision :: anorm
    double precision, allocatable :: wc_a(:,:) ! working copy of a
    integer :: info
    integer, allocatable :: ipiv(:)
    
    double precision, external :: dlange
    
    
     if (present(status)) status=1
    
    anorm=dlange('1',d,d,a,d,work) ! work is unallocated as it is only used when computing infinity norm
    
    allocate(wc_a(d,d))
    !call dcopy(d*d,a,1,wc_a,1)
    wc_a(:,:)=a(:,:)
    
    allocate(ipiv(d))
    call dgetrf(d,d,wc_a,d,ipiv,info)
    if (info /= 0) then
         if (present(status)) status=0
        deallocate(ipiv)
        deallocate(wc_a)
        return
    end if

    allocate(work(4*d))
    allocate(iwork(d))
    call dgecon('1',d,wc_a,d,anorm,rcond,work,iwork,info)
    if (info /= 0) then
         if (present(status)) status=0
    end if
    deallocate(iwork)
    deallocate(work)
    deallocate(ipiv)
    deallocate(wc_a)

end subroutine

subroutine fvn_c_matcon(d,a,rcond,status)
    ! Matrix condition (reciprocal of condition number)
    ! 
    ! d (in) : matrix rank
    ! a (in) : The Matrix
    ! rcond (out) : guess what
    ! status (out) : =0 if something went wrong
    !
    implicit none
    integer, intent(in) :: d
    complex, intent(in) :: a(d,d)
    real, intent(out) :: rcond
    integer, intent(out), optional :: status
    
    real, allocatable :: rwork(:)
    complex, allocatable :: work(:)
    integer, allocatable :: iwork(:)
    real :: anorm
    complex, allocatable :: wc_a(:,:) ! working copy of a
    integer :: info
    integer, allocatable :: ipiv(:)
    
    real, external :: clange
    
    
     if (present(status)) status=1
    
    anorm=clange('1',d,d,a,d,rwork) ! rwork is unallocated as it is only used when computing infinity norm
    
    allocate(wc_a(d,d))
    !call ccopy(d*d,a,1,wc_a,1)
    wc_a(:,:)=a(:,:)
    
    allocate(ipiv(d))
    call cgetrf(d,d,wc_a,d,ipiv,info)
    if (info /= 0) then
         if (present(status)) status=0
        deallocate(ipiv)
        deallocate(wc_a)
        return
    end if
    allocate(work(2*d))
    allocate(rwork(2*d))
    call cgecon('1',d,wc_a,d,anorm,rcond,work,rwork,info)
    if (info /= 0) then
         if (present(status)) status=0
    end if
    deallocate(rwork)
    deallocate(work)
    deallocate(ipiv)
    deallocate(wc_a)
end subroutine

subroutine fvn_z_matcon(d,a,rcond,status)
    ! Matrix condition (reciprocal of condition number)
    ! 
    ! d (in) : matrix rank
    ! a (in) : The Matrix
    ! rcond (out) : guess what
    ! status (out) : =0 if something went wrong
    !
    implicit none
    integer, intent(in) :: d
    double complex, intent(in) :: a(d,d)
    double precision, intent(out) :: rcond
    integer, intent(out), optional :: status
    
    double complex, allocatable :: work(:)
    double precision, allocatable :: rwork(:)
    double precision :: anorm
    double complex, allocatable :: wc_a(:,:) ! working copy of a
    integer :: info
    integer, allocatable :: ipiv(:)
    
    double precision, external :: zlange
    
    
     if (present(status)) status=1
    
    anorm=zlange('1',d,d,a,d,rwork) ! rwork is unallocated as it is only used when computing infinity norm
    
    allocate(wc_a(d,d))
    !call zcopy(d*d,a,1,wc_a,1)
    wc_a(:,:)=a(:,:)
    
    allocate(ipiv(d))
    call zgetrf(d,d,wc_a,d,ipiv,info)
    if (info /= 0) then
         if (present(status)) status=0
        deallocate(ipiv)
        deallocate(wc_a)
        return
    end if

    allocate(work(2*d))
    allocate(rwork(2*d))
    call zgecon('1',d,wc_a,d,anorm,rcond,work,rwork,info)
    if (info /= 0) then
         if (present(status)) status=0
    end if
    deallocate(rwork)
    deallocate(work)
    deallocate(ipiv)
    deallocate(wc_a)
end subroutine

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! Valeurs propres/ Vecteurs propre
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

subroutine fvn_s_matev(d,a,evala,eveca,status)
    ! 
    ! integer d (in) : matrice rank
    ! real a(d,d) (in) : The Matrix
    ! complex evala(d) (out) : eigenvalues
    ! complex eveca(d,d) (out) : eveca(:,j) = jth eigenvector
    ! integer (out) : status =0 if something went wrong
    !
    ! interfacing Lapack routine SGEEV
    implicit none
    integer, intent(in) :: d
    real, intent(in) :: a(d,d)
    complex, intent(out) :: evala(d)
    complex, intent(out) :: eveca(d,d)
    integer, intent(out), optional :: status
    
    real, allocatable :: wc_a(:,:)  ! a working copy of a
    integer :: info
    integer :: lwork
    real, allocatable :: wr(:),wi(:)
    real :: vl      ! unused but necessary for the call
    real, allocatable :: vr(:,:)
    real, allocatable :: work(:)
    real :: twork(1)
    integer i
    integer j
    
     if (present(status)) status=1

    ! making a working copy of a
    allocate(wc_a(d,d))
    !call scopy(d*d,a,1,wc_a,1)
    wc_a(:,:)=a(:,:)
    
    allocate(wr(d))
    allocate(wi(d))
    allocate(vr(d,d))
    ! query optimal work size
    call sgeev('N','V',d,wc_a,d,wr,wi,vl,1,vr,d,twork,-1,info)
    lwork=int(twork(1))
    allocate(work(lwork))
    call sgeev('N','V',d,wc_a,d,wr,wi,vl,1,vr,d,work,lwork,info)

    if (info /= 0) then
         if (present(status)) status=0
        deallocate(work)
        deallocate(vr)
        deallocate(wi)
        deallocate(wr)
        deallocate(wc_a)
        return
    end if

    ! now fill in the results
    i=1
    do while(i<=d)
        evala(i)=cmplx(wr(i),wi(i))
        if (wi(i) == 0.) then ! eigenvalue is real
            eveca(:,i)=cmplx(vr(:,i),0.)
        else ! eigenvalue is complex
            evala(i+1)=cmplx(wr(i+1),wi(i+1))
            eveca(:,i)=cmplx(vr(:,i),vr(:,i+1))
            eveca(:,i+1)=cmplx(vr(:,i),-vr(:,i+1))
            i=i+1
        end if
        i=i+1
    enddo
    deallocate(work)
    deallocate(vr)
    deallocate(wi)
    deallocate(wr)
    deallocate(wc_a)

end subroutine

subroutine fvn_d_matev(d,a,evala,eveca,status)
    ! 
    ! integer d (in) : matrice rank
    ! double precision a(d,d) (in) : The Matrix
    ! double complex evala(d) (out) : eigenvalues
    ! double complex eveca(d,d) (out) : eveca(:,j) = jth eigenvector
    ! integer (out) : status =0 if something went wrong
    !
    ! interfacing Lapack routine DGEEV
    implicit none
    integer, intent(in) :: d
    double precision, intent(in) :: a(d,d)
    double complex, intent(out) :: evala(d)
    double complex, intent(out) :: eveca(d,d)
    integer, intent(out), optional :: status
    
    double precision, allocatable :: wc_a(:,:)  ! a working copy of a
    integer :: info
    integer :: lwork
    double precision, allocatable :: wr(:),wi(:)
    double precision :: vl      ! unused but necessary for the call
    double precision, allocatable :: vr(:,:)
    double precision, allocatable :: work(:)
    double precision :: twork(1)
    integer i
    integer j
    
    if (present(status)) status=1

    ! making a working copy of a
    allocate(wc_a(d,d))
    !call dcopy(d*d,a,1,wc_a,1)
    wc_a(:,:)=a(:,:)
    
    allocate(wr(d))
    allocate(wi(d))
    allocate(vr(d,d))
    ! query optimal work size
    call dgeev('N','V',d,wc_a,d,wr,wi,vl,1,vr,d,twork,-1,info)
    lwork=int(twork(1))
    allocate(work(lwork))
    call dgeev('N','V',d,wc_a,d,wr,wi,vl,1,vr,d,work,lwork,info)

    if (info /= 0) then
         if (present(status)) status=0
        deallocate(work)
        deallocate(vr)
        deallocate(wi)
        deallocate(wr)
        deallocate(wc_a)
        return
    end if

    ! now fill in the results
    i=1
    do while(i<=d)
        evala(i)=dcmplx(wr(i),wi(i))
        if (wi(i) == 0.) then ! eigenvalue is real
            eveca(:,i)=dcmplx(vr(:,i),0.)
        else ! eigenvalue is complex
            evala(i+1)=dcmplx(wr(i+1),wi(i+1))
            eveca(:,i)=dcmplx(vr(:,i),vr(:,i+1))
            eveca(:,i+1)=dcmplx(vr(:,i),-vr(:,i+1))
            i=i+1
        end if
        i=i+1
    enddo

    deallocate(work)
    deallocate(vr)
    deallocate(wi)
    deallocate(wr)
    deallocate(wc_a)

end subroutine

subroutine fvn_c_matev(d,a,evala,eveca,status)
    ! 
    ! integer d (in) : matrice rank
    ! complex a(d,d) (in) : The Matrix
    ! complex evala(d) (out) : eigenvalues
    ! complex eveca(d,d) (out) : eveca(:,j) = jth eigenvector
    ! integer (out) : status =0 if something went wrong
    !
    ! interfacing Lapack routine CGEEV
    implicit none
    integer, intent(in) :: d
    complex, intent(in) :: a(d,d)
    complex, intent(out) :: evala(d)
    complex, intent(out) :: eveca(d,d)
    integer, intent(out), optional :: status

    complex, allocatable :: wc_a(:,:) ! a working copy of a
    integer :: info
    integer :: lwork
    complex, allocatable :: work(:)
    complex :: twork(1)
    real, allocatable :: rwork(:)
    complex :: vl   ! unused but necessary for the call

     if (present(status)) status=1

    ! making a working copy of a
    allocate(wc_a(d,d))
    !call ccopy(d*d,a,1,wc_a,1)
    wc_a(:,:)=a(:,:)
    
    
    ! query optimal work size
    call cgeev('N','V',d,wc_a,d,evala,vl,1,eveca,d,twork,-1,rwork,info)
    lwork=int(twork(1))
    allocate(work(lwork))
    allocate(rwork(2*d))
    call cgeev('N','V',d,wc_a,d,evala,vl,1,eveca,d,work,lwork,rwork,info)
    
    if (info /= 0) then
         if (present(status)) status=0
    end if
    deallocate(rwork)
    deallocate(work)
    deallocate(wc_a)

end subroutine

subroutine fvn_z_matev(d,a,evala,eveca,status)
    ! 
    ! integer d (in) : matrice rank
    ! double complex a(d,d) (in) : The Matrix
    ! double complex evala(d) (out) : eigenvalues
    ! double complex eveca(d,d) (out) : eveca(:,j) = jth eigenvector
    ! integer (out) : status =0 if something went wrong
    !
    ! interfacing Lapack routine ZGEEV
    implicit none
    integer, intent(in) :: d
    double complex, intent(in) :: a(d,d)
    double complex, intent(out) :: evala(d)
    double complex, intent(out) :: eveca(d,d)
    integer, intent(out), optional :: status

    double complex, allocatable :: wc_a(:,:) ! a working copy of a
    integer :: info
    integer :: lwork
    double complex, allocatable :: work(:)
    double complex :: twork(1)
    double precision, allocatable :: rwork(:)
    double complex :: vl   ! unused but necessary for the call
    
     if (present(status)) status=1

    ! making a working copy of a
    allocate(wc_a(d,d))
    !call zcopy(d*d,a,1,wc_a,1)
    wc_a(:,:)=a(:,:)
    
    ! query optimal work size
    call zgeev('N','V',d,wc_a,d,evala,vl,1,eveca,d,twork,-1,rwork,info)
    lwork=int(twork(1))
    allocate(work(lwork))
    allocate(rwork(2*d))
    call zgeev('N','V',d,wc_a,d,evala,vl,1,eveca,d,work,lwork,rwork,info)
    
    if (info /= 0) then
         if (present(status)) status=0
    end if
    deallocate(rwork)
    deallocate(work)
    deallocate(wc_a)

end subroutine


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! Quadratic interpolation of tabulated function of 1,2 or 3 variables
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

subroutine fvn_s_find_interval(x,i,xdata,n)
      implicit none
      ! This routine find the indice i where xdata(i) <= x < xdata(i+1)
      ! xdata(n) must contains a set of increasingly ordered values 
      ! if x < xdata(1) i=0 is returned
      ! if x > xdata(n) i=n is returned
      ! special case is where x=xdata(n) then n-1 is returned so 
      ! we will not exclude the upper limit
      ! a simple dichotomy method is used

      real(kind=4), intent(in) :: x
      real(kind=4), intent(in), dimension(n) :: xdata
      integer(kind=4), intent(in) :: n
      integer(kind=4), intent(out) :: i

      integer(kind=4) :: imin,imax,imoyen

      ! special case is where x=xdata(n) then n-1 is returned so 
      ! we will not exclude the upper limit
      if (x == xdata(n)) then
            i=n-1
            return
      end if

      ! if x < xdata(1) i=0 is returned
      if (x < xdata(1)) then
            i=0
            return
      end if

      ! if x > xdata(n) i=n is returned
      if (x > xdata(n)) then
            i=n
            return
      end if

      ! here xdata(1) <= x <= xdata(n)
      imin=0
      imax=n+1

      do while((imax-imin) > 1)
            imoyen=(imax+imin)/2
            if (x >= xdata(imoyen)) then
                  imin=imoyen
            else
                  imax=imoyen
            end if
      end do

      i=imin

end subroutine


subroutine fvn_d_find_interval(x,i,xdata,n)
      implicit none
      ! This routine find the indice i where xdata(i) <= x < xdata(i+1)
      ! xdata(n) must contains a set of increasingly ordered values 
      ! if x < xdata(1) i=0 is returned
      ! if x > xdata(n) i=n is returned
      ! special case is where x=xdata(n) then n-1 is returned so 
      ! we will not exclude the upper limit
      ! a simple dichotomy method is used

      real(kind=8), intent(in) :: x
      real(kind=8), intent(in), dimension(n) :: xdata
      integer(kind=4), intent(in) :: n
      integer(kind=4), intent(out) :: i

      integer(kind=4) :: imin,imax,imoyen

      ! special case is where x=xdata(n) then n-1 is returned so 
      ! we will not exclude the upper limit
      if (x == xdata(n)) then
            i=n-1
            return
      end if

      ! if x < xdata(1) i=0 is returned
      if (x < xdata(1)) then
            i=0
            return
      end if

      ! if x > xdata(n) i=n is returned
      if (x > xdata(n)) then
            i=n
            return
      end if

      ! here xdata(1) <= x <= xdata(n)
      imin=0
      imax=n+1

      do while((imax-imin) > 1)
            imoyen=(imax+imin)/2
            if (x >= xdata(imoyen)) then
                  imin=imoyen
            else
                  imax=imoyen
            end if
      end do

      i=imin

end subroutine


function fvn_s_quad_interpol(x,n,xdata,ydata)
      implicit none
      ! This function evaluate the value of a function defined by a set of points
      ! and values, using a quadratic interpolation
      ! xdata must be increasingly ordered
      ! x must be within xdata(1) and xdata(n) to actually do interpolation
      ! otherwise extrapolation is done
      integer(kind=4), intent(in) :: n
      real(kind=4), intent(in), dimension(n) :: xdata,ydata
      real(kind=4), intent(in) :: x
      real(kind=4) ::  fvn_s_quad_interpol

      integer(kind=4) :: iinf,base,i,j
      real(kind=4) :: p

      call fvn_s_find_interval(x,iinf,xdata,n)

      ! Settings for extrapolation
      if (iinf==0) then
            ! TODO -> Lower bound extrapolation warning
            iinf=1
      end if

      if (iinf==n) then
            ! TODO -> Higher bound extrapolation warning
            iinf=n-1
      end if

      ! The three points we will use are iinf-1,iinf and iinf+1 with the
      ! exception of the first interval, where iinf=1 we will use 1,2 and 3
      if (iinf==1) then
            base=0
      else
            base=iinf-2
      end if

      ! The three points we will use are :
      ! xdata/ydata(base+1),xdata/ydata(base+2),xdata/ydata(base+3)

      ! Straight forward Lagrange polynomial
      fvn_s_quad_interpol=0.
      do i=1,3
            !  polynome i
            p=ydata(base+i)
            do j=1,3
                  if (j /= i) then
                        p=p*(x-xdata(base+j))/(xdata(base+i)-xdata(base+j))
                  end if
            end do
            fvn_s_quad_interpol=fvn_s_quad_interpol+p
      end do

end function


function fvn_d_quad_interpol(x,n,xdata,ydata)
      implicit none
      ! This function evaluate the value of a function defined by a set of points
      ! and values, using a quadratic interpolation
      ! xdata must be increasingly ordered
      ! x must be within xdata(1) and xdata(n) to actually do interpolation
      ! otherwise extrapolation is done
      integer(kind=4), intent(in) :: n
      real(kind=8), intent(in), dimension(n) :: xdata,ydata
      real(kind=8), intent(in) :: x
      real(kind=8) ::  fvn_d_quad_interpol

      integer(kind=4) :: iinf,base,i,j
      real(kind=8) :: p

      call fvn_d_find_interval(x,iinf,xdata,n)

      ! Settings for extrapolation
      if (iinf==0) then
            ! TODO -> Lower bound extrapolation warning
            iinf=1
      end if

      if (iinf==n) then
            ! TODO Higher bound extrapolation warning
            iinf=n-1
      end if

      ! The three points we will use are iinf-1,iinf and iinf+1 with the
      ! exception of the first interval, where iinf=1 we will use 1,2 and 3
      if (iinf==1) then
            base=0
      else
            base=iinf-2
      end if

      ! The three points we will use are :
      ! xdata/ydata(base+1),xdata/ydata(base+2),xdata/ydata(base+3)

      ! Straight forward Lagrange polynomial
      fvn_d_quad_interpol=0.
      do i=1,3
            !  polynome i
            p=ydata(base+i)
            do j=1,3
                  if (j /= i) then
                        p=p*(x-xdata(base+j))/(xdata(base+i)-xdata(base+j))
                  end if
            end do
            fvn_d_quad_interpol=fvn_d_quad_interpol+p
      end do

end function


function fvn_s_quad_2d_interpol(x,y,nx,xdata,ny,ydata,zdata)
      implicit none
      ! This function evaluate the value of a two variable function defined by a
      ! set of points and values, using a quadratic interpolation
      ! xdata and ydata must be increasingly ordered
      ! the couple (x,y) must be as x within xdata(1) and xdata(nx) and
      ! y within ydata(1) and ydata(ny) to actually do interpolation
      ! otherwise extrapolation is done
      integer(kind=4), intent(in) :: nx,ny
      real(kind=4), intent(in) :: x,y
      real(kind=4), intent(in), dimension(nx) :: xdata
      real(kind=4), intent(in), dimension(ny) :: ydata
      real(kind=4), intent(in), dimension(nx,ny) :: zdata
      real(kind=4) :: fvn_s_quad_2d_interpol

      integer(kind=4) :: ixinf,iyinf,basex,basey,i
      real(kind=4),dimension(3) :: ztmp
      !real(kind=4), external :: fvn_s_quad_interpol

      call fvn_s_find_interval(x,ixinf,xdata,nx)
      call fvn_s_find_interval(y,iyinf,ydata,ny)

      ! Settings for extrapolation
      if (ixinf==0) then
            ! TODO -> Lower x  bound extrapolation warning
            ixinf=1
      end if

      if (ixinf==nx) then
            ! TODO -> Higher x bound extrapolation warning
            ixinf=nx-1
      end if

      if (iyinf==0) then
            ! TODO -> Lower y  bound extrapolation warning
            iyinf=1
      end if

      if (iyinf==ny) then
            ! TODO -> Higher y bound extrapolation warning
            iyinf=ny-1
      end if

      ! The three points we will use are iinf-1,iinf and iinf+1 with the
      ! exception of the first interval, where iinf=1 we will use 1,2 and 3
      if (ixinf==1) then
            basex=0
      else
            basex=ixinf-2
      end if

      if (iyinf==1) then
            basey=0
      else
            basey=iyinf-2
      end if

      ! First we make 3 interpolations for x at y(base+1),y(base+2),y(base+3)
      ! stored in ztmp(1:3)
      do i=1,3
            ztmp(i)=fvn_s_quad_interpol(x,nx,xdata,zdata(:,basey+i))
      end do

      ! Then we make an interpolation for y using previous interpolations
      fvn_s_quad_2d_interpol=fvn_s_quad_interpol(y,3,ydata(basey+1:basey+3),ztmp)
end function


function fvn_d_quad_2d_interpol(x,y,nx,xdata,ny,ydata,zdata)
      implicit none
      ! This function evaluate the value of a two variable function defined by a
      ! set of points and values, using a quadratic interpolation
      ! xdata and ydata must be increasingly ordered
      ! the couple (x,y) must be as x within xdata(1) and xdata(nx) and
      ! y within ydata(1) and ydata(ny) to actually do interpolation
      ! otherwise extrapolation is done
      integer(kind=4), intent(in) :: nx,ny
      real(kind=8), intent(in) :: x,y
      real(kind=8), intent(in), dimension(nx) :: xdata
      real(kind=8), intent(in), dimension(ny) :: ydata
      real(kind=8), intent(in), dimension(nx,ny) :: zdata
      real(kind=8) :: fvn_d_quad_2d_interpol

      integer(kind=4) :: ixinf,iyinf,basex,basey,i
      real(kind=8),dimension(3) :: ztmp
      !real(kind=8), external :: fvn_d_quad_interpol

      call fvn_d_find_interval(x,ixinf,xdata,nx)
      call fvn_d_find_interval(y,iyinf,ydata,ny)

      ! Settings for extrapolation
      if (ixinf==0) then
            ! TODO -> Lower x  bound extrapolation warning
            ixinf=1
      end if

      if (ixinf==nx) then
            ! TODO -> Higher x bound extrapolation warning
            ixinf=nx-1
      end if

      if (iyinf==0) then
            ! TODO -> Lower y  bound extrapolation warning
            iyinf=1
      end if

      if (iyinf==ny) then
            ! TODO -> Higher y bound extrapolation warning
            iyinf=ny-1
      end if

      ! The three points we will use are iinf-1,iinf and iinf+1 with the
      ! exception of the first interval, where iinf=1 we will use 1,2 and 3
      if (ixinf==1) then
            basex=0
      else
            basex=ixinf-2
      end if

      if (iyinf==1) then
            basey=0
      else
            basey=iyinf-2
      end if

      ! First we make 3 interpolations for x at y(base+1),y(base+2),y(base+3)
      ! stored in ztmp(1:3)
      do i=1,3
            ztmp(i)=fvn_d_quad_interpol(x,nx,xdata,zdata(:,basey+i))
      end do

      ! Then we make an interpolation for y using previous interpolations
      fvn_d_quad_2d_interpol=fvn_d_quad_interpol(y,3,ydata(basey+1:basey+3),ztmp)
end function


function fvn_s_quad_3d_interpol(x,y,z,nx,xdata,ny,ydata,nz,zdata,tdata)
      implicit none
      ! This function evaluate the value of a 3 variables function defined by a
      ! set of points and values, using a quadratic interpolation
      ! xdata, ydata and zdata must be increasingly ordered
      ! The triplet (x,y,z) must be within xdata,ydata and zdata to actually 
      ! perform an interpolation, otherwise extrapolation is done
      integer(kind=4), intent(in) :: nx,ny,nz
      real(kind=4), intent(in) :: x,y,z
      real(kind=4), intent(in), dimension(nx) :: xdata
      real(kind=4), intent(in), dimension(ny) :: ydata
      real(kind=4), intent(in), dimension(nz) :: zdata
      real(kind=4), intent(in), dimension(nx,ny,nz) :: tdata
      real(kind=4) :: fvn_s_quad_3d_interpol

      integer(kind=4) :: ixinf,iyinf,izinf,basex,basey,basez,i,j
      !real(kind=4), external :: fvn_s_quad_interpol,fvn_s_quad_2d_interpol
      real(kind=4),dimension(3,3) :: ttmp

      call fvn_s_find_interval(x,ixinf,xdata,nx)
      call fvn_s_find_interval(y,iyinf,ydata,ny)
      call fvn_s_find_interval(z,izinf,zdata,nz)

      ! Settings for extrapolation
      if (ixinf==0) then
            ! TODO -> Lower x  bound extrapolation warning
            ixinf=1
      end if

      if (ixinf==nx) then
            ! TODO -> Higher x bound extrapolation warning
            ixinf=nx-1
      end if

      if (iyinf==0) then
            ! TODO -> Lower y  bound extrapolation warning
            iyinf=1
      end if

      if (iyinf==ny) then
            ! TODO -> Higher y bound extrapolation warning
            iyinf=ny-1
      end if

      if (izinf==0) then
            ! TODO -> Lower z bound extrapolation warning
            izinf=1
      end if

      if (izinf==nz) then
            ! TODO -> Higher z bound extrapolation warning
            izinf=nz-1
      end if

      ! The three points we will use are iinf-1,iinf and iinf+1 with the
      ! exception of the first interval, where iinf=1 we will use 1,2 and 3
      if (ixinf==1) then
            basex=0
      else
            basex=ixinf-2
      end if

      if (iyinf==1) then
            basey=0
      else
            basey=iyinf-2
      end if

      if (izinf==1) then
            basez=0
      else
            basez=izinf-2
      end if

      ! We first make 9 one dimensional interpolation on variable x.
      ! results are stored in ttmp
      do i=1,3
            do j=1,3
                  ttmp(i,j)=fvn_s_quad_interpol(x,nx,xdata,tdata(:,basey+i,basez+j))
            end do
      end do

      ! We then make a 2 dimensionnal interpolation on variables y and z
      fvn_s_quad_3d_interpol=fvn_s_quad_2d_interpol(y,z, &
            3,ydata(basey+1:basey+3),3,zdata(basez+1:basez+3),ttmp)
end function


function fvn_d_quad_3d_interpol(x,y,z,nx,xdata,ny,ydata,nz,zdata,tdata)
      implicit none
      ! This function evaluate the value of a 3 variables function defined by a
      ! set of points and values, using a quadratic interpolation
      ! xdata, ydata and zdata must be increasingly ordered
      ! The triplet (x,y,z) must be within xdata,ydata and zdata to actually 
      ! perform an interpolation, otherwise extrapolation is done
      integer(kind=4), intent(in) :: nx,ny,nz
      real(kind=8), intent(in) :: x,y,z
      real(kind=8), intent(in), dimension(nx) :: xdata
      real(kind=8), intent(in), dimension(ny) :: ydata
      real(kind=8), intent(in), dimension(nz) :: zdata
      real(kind=8), intent(in), dimension(nx,ny,nz) :: tdata
      real(kind=8) :: fvn_d_quad_3d_interpol

      integer(kind=4) :: ixinf,iyinf,izinf,basex,basey,basez,i,j
      !real(kind=8), external :: fvn_d_quad_interpol,fvn_d_quad_2d_interpol
      real(kind=8),dimension(3,3) :: ttmp

      call fvn_d_find_interval(x,ixinf,xdata,nx)
      call fvn_d_find_interval(y,iyinf,ydata,ny)
      call fvn_d_find_interval(z,izinf,zdata,nz)

      ! Settings for extrapolation
      if (ixinf==0) then
            ! TODO -> Lower x  bound extrapolation warning
            ixinf=1
      end if

      if (ixinf==nx) then
            ! TODO -> Higher x bound extrapolation warning
            ixinf=nx-1
      end if

      if (iyinf==0) then
            ! TODO -> Lower y  bound extrapolation warning
            iyinf=1
      end if

      if (iyinf==ny) then
            ! TODO -> Higher y bound extrapolation warning
            iyinf=ny-1
      end if

      if (izinf==0) then
            ! TODO -> Lower z bound extrapolation warning
            izinf=1
      end if

      if (izinf==nz) then
            ! TODO -> Higher z bound extrapolation warning
            izinf=nz-1
      end if

      ! The three points we will use are iinf-1,iinf and iinf+1 with the
      ! exception of the first interval, where iinf=1 we will use 1,2 and 3
      if (ixinf==1) then
            basex=0
      else
            basex=ixinf-2
      end if

      if (iyinf==1) then
            basey=0
      else
            basey=iyinf-2
      end if

      if (izinf==1) then
            basez=0
      else
            basez=izinf-2
      end if

      ! We first make 9 one dimensional interpolation on variable x.
      ! results are stored in ttmp
      do i=1,3
            do j=1,3
                  ttmp(i,j)=fvn_d_quad_interpol(x,nx,xdata,tdata(:,basey+i,basez+j))
            end do
      end do

      ! We then make a 2 dimensionnal interpolation on variables y and z
      fvn_d_quad_3d_interpol=fvn_d_quad_2d_interpol(y,z, &
            3,ydata(basey+1:basey+3),3,zdata(basez+1:basez+3),ttmp)
end function




!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! Akima spline interpolation and spline evaluation
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

! Single precision
subroutine fvn_s_akima(n,x,y,br,co)
    implicit none
    integer, intent(in)  :: n
    real, intent(in) :: x(n)
    real, intent(in) :: y(n)
    real, intent(out) :: br(n)
    real, intent(out) :: co(4,n)
    
    real, allocatable :: var(:),z(:)
    real :: wi_1,wi
    integer :: i
    real :: dx,a,b

    ! br is just a copy of x
    br(:)=x(:)
    
    allocate(var(n))
    allocate(z(n))
    ! evaluate the variations
    do i=1, n-1
        var(i+2)=(y(i+1)-y(i))/(x(i+1)-x(i))
    end do
    var(n+2)=2.e0*var(n+1)-var(n)
    var(n+3)=2.e0*var(n+2)-var(n+1)
    var(2)=2.e0*var(3)-var(4)
    var(1)=2.e0*var(2)-var(3)
  
    do i = 1, n
    wi_1=abs(var(i+3)-var(i+2))
    wi=abs(var(i+1)-var(i))
    if ((wi_1+wi).eq.0.e0) then
        z(i)=(var(i+2)+var(i+1))/2.e0
    else
        z(i)=(wi_1*var(i+1)+wi*var(i+2))/(wi_1+wi)
    end if
    end do
    
    do i=1, n-1
        dx=x(i+1)-x(i)
        a=(z(i+1)-z(i))*dx                      ! coeff intermediaires pour calcul wd
        b=y(i+1)-y(i)-z(i)*dx                   ! coeff intermediaires pour calcul wd
        co(1,i)=y(i)
        co(2,i)=z(i)
        !co(3,i)=-(a-3.*b)/dx**2                ! méthode wd
        !co(4,i)=(a-2.*b)/dx**3                 ! méthode wd
        co(3,i)=(3.e0*var(i+2)-2.e0*z(i)-z(i+1))/dx   ! méthode JP Moreau
        co(4,i)=(z(i)+z(i+1)-2.e0*var(i+2))/dx**2  !
        ! les coefficients donnés par imsl sont co(3,i)*2 et co(4,i)*6
        ! etrangement la fonction csval corrige et donne la bonne valeur ...
    end do
    co(1,n)=y(n)
    co(2,n)=z(n)
    co(3,n)=0.e0
    co(4,n)=0.e0

    deallocate(z)
    deallocate(var)

end subroutine

! Double precision
subroutine fvn_d_akima(n,x,y,br,co)

    implicit none
    integer, intent(in)  :: n
    double precision, intent(in) :: x(n)
    double precision, intent(in) :: y(n)
    double precision, intent(out) :: br(n)
    double precision, intent(out) :: co(4,n)
    
    double precision, allocatable :: var(:),z(:)
    double precision :: wi_1,wi
    integer :: i
    double precision :: dx,a,b
    
    ! br is just a copy of x
    br(:)=x(:)

    allocate(var(n))
    allocate(z(n))
    ! evaluate the variations
    do i=1, n-1
        var(i+2)=(y(i+1)-y(i))/(x(i+1)-x(i))
    end do
    var(n+2)=2.d0*var(n+1)-var(n)
    var(n+3)=2.d0*var(n+2)-var(n+1)
    var(2)=2.d0*var(3)-var(4)
    var(1)=2.d0*var(2)-var(3)
  
    do i = 1, n
    wi_1=dabs(var(i+3)-var(i+2))
    wi=dabs(var(i+1)-var(i))
    if ((wi_1+wi).eq.0.d0) then
        z(i)=(var(i+2)+var(i+1))/2.d0
    else
        z(i)=(wi_1*var(i+1)+wi*var(i+2))/(wi_1+wi)
    end if
    end do
    
    do i=1, n-1
        dx=x(i+1)-x(i)
        a=(z(i+1)-z(i))*dx                      ! coeff intermediaires pour calcul wd
        b=y(i+1)-y(i)-z(i)*dx                   ! coeff intermediaires pour calcul wd
        co(1,i)=y(i)
        co(2,i)=z(i)
        !co(3,i)=-(a-3.*b)/dx**2                ! méthode wd
        !co(4,i)=(a-2.*b)/dx**3                 ! méthode wd
        co(3,i)=(3.d0*var(i+2)-2.d0*z(i)-z(i+1))/dx   ! méthode JP Moreau
        co(4,i)=(z(i)+z(i+1)-2.d0*var(i+2))/dx**2  !
        ! les coefficients donnés par imsl sont co(3,i)*2 et co(4,i)*6
        ! etrangement la fonction csval corrige et donne la bonne valeur ...
    end do
    co(1,n)=y(n)
    co(2,n)=z(n)
    co(3,n)=0.d0
    co(4,n)=0.d0

    deallocate(z)
    deallocate(var)

end subroutine

!
! Single precision spline evaluation
!
function fvn_s_spline_eval(x,n,br,co)
    implicit none
    real, intent(in) :: x           ! x must be br(1)<= x <= br(n+1) otherwise value is extrapolated
    integer, intent(in) :: n        ! number of intervals
    real, intent(in) :: br(n+1)     ! breakpoints
    real, intent(in) :: co(4,n+1)   ! spline coeeficients
    real :: fvn_s_spline_eval
    
    integer :: i
    real :: dx
    
    if (x<=br(1)) then
        i=1
    else if (x>=br(n+1)) then
        i=n
    else
    i=1
    do while(x>=br(i))
        i=i+1
    end do
    i=i-1
    end if
    dx=x-br(i)
    fvn_s_spline_eval=co(1,i)+co(2,i)*dx+co(3,i)*dx**2+co(4,i)*dx**3

end function

! Double precision spline evaluation
function fvn_d_spline_eval(x,n,br,co)
    implicit none
    double precision, intent(in) :: x           ! x must be br(1)<= x <= br(n+1) otherwise value is extrapolated
    integer, intent(in) :: n        ! number of intervals
    double precision, intent(in) :: br(n+1)     ! breakpoints
    double precision, intent(in) :: co(4,n+1)   ! spline coeeficients
    double precision :: fvn_d_spline_eval
    
    integer :: i
    double precision :: dx
    
    
    if (x<=br(1)) then
        i=1
    else if (x>=br(n+1)) then
        i=n
    else
    i=1
    do while(x>=br(i))
        i=i+1
    end do
    i=i-1
    end if
    
    dx=x-br(i)
    fvn_d_spline_eval=co(1,i)+co(2,i)*dx+co(3,i)*dx**2+co(4,i)*dx**3

end function


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! 
! Least square problem
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
!




subroutine fvn_d_lspoly(np,x,y,deg,coeff,status)
!
!   Least square polynomial fitting
!
!   Find the coefficients of the least square polynomial of a given degree
!   for a set of coordinates.
!
!   The degree must be lower than the number of points
!
!   np (in)             : number of points
!   x(np) (in)          : x data
!   y(np) (in)          : y data
!   deg (in)            : polynomial's degree
!   coeff(deg+1) (out)  : polynomial coefficients
!   status (out)        : =0 if a problem occurs
!
implicit none

integer, intent(in) :: np,deg
real(kind=8), intent(in), dimension(np) :: x,y
real(kind=8), intent(out), dimension(deg+1) :: coeff
integer, intent(out), optional :: status

real(kind=8), allocatable, dimension(:,:) :: mat,bmat
real(kind=8),dimension(:),allocatable :: work
real(kind=8),dimension(1) :: twork
integer :: lwork,info

integer :: i,j

 if (present(status)) status=1
allocate(mat(np,deg+1),bmat(np,1))

! Design matrix valorisation
mat=reshape( (/ ((x(i)**(j-1),i=1,np),j=1,deg+1) /),shape=(/ np,deg+1 /) )

! second member valorisation
bmat=reshape ( (/  (y(i),i=1,np)   /)  ,shape = (/ np,1 /))

! query workspace size
call dgels('N',np,deg+1,1,mat,np,bmat,np,twork,-1,info)
lwork=twork(1)
allocate(work(int(lwork)))
! real call
call dgels('N',np,deg+1,1,mat,np,bmat,np,work,lwork,info)

if (info /= 0) then
     if (present(status)) status=0
end if

 coeff = (/ (bmat(i,1),i=1,deg+1) /)

deallocate(work)
deallocate(mat,bmat)
end subroutine

subroutine fvn_s_lspoly(np,x,y,deg,coeff,status)
!
!   Least square polynomial fitting
!
!   Find the coefficients of the least square polynomial of a given degree
!   for a set of coordinates.
!
!   The degree must be lower than the number of points
!
!   np (in)             : number of points
!   x(np) (in)          : x data
!   y(np) (in)          : y data
!   deg (in)            : polynomial's degree
!   coeff(deg+1) (out)  : polynomial coefficients
!   status (out)        : =0 if a problem occurs
!
implicit none

integer, intent(in) :: np,deg
real(kind=4), intent(in), dimension(np) :: x,y
real(kind=4), intent(out), dimension(deg+1) :: coeff
integer, intent(out), optional :: status

real(kind=4), allocatable, dimension(:,:) :: mat,bmat
real(kind=4),dimension(:),allocatable :: work
real(kind=4),dimension(1) :: twork
integer :: lwork,info

integer :: i,j

 if (present(status)) status=1
allocate(mat(np,deg+1),bmat(np,1))

! Design matrix valorisation
mat=reshape( (/ ((x(i)**(j-1),i=1,np),j=1,deg+1) /),shape=(/ np,deg+1 /) )

! second member valorisation
bmat=reshape ( (/  (y(i),i=1,np)   /)  ,shape = (/ np,1 /))

! query workspace size
call sgels('N',np,deg+1,1,mat,np,bmat,np,twork,-1,info)
lwork=twork(1)
allocate(work(int(lwork)))
! real call
call sgels('N',np,deg+1,1,mat,np,bmat,np,work,lwork,info)

if (info /= 0) then
     if (present(status)) status=0
end if

 coeff = (/ (bmat(i,1),i=1,deg+1) /)

deallocate(work)
deallocate(mat,bmat)
end subroutine








subroutine fvn_d_lspoly_svd(np,x,y,deg,coeff,status)
!
!   Least square polynomial fitting using singular value decomposition
!
!   Find the coefficients of the least square polynomial of a given degree
!   for a set of coordinates.
!
!   The degree must be lower than the number of points
!
!   np (in)             : number of points
!   x(np) (in)          : x data
!   y(np) (in)          : y data
!   deg (in)            : polynomial's degree
!   coeff(deg+1) (out)  : polynomial coefficients
!   status (out)        : =0 if a problem occurs
!
implicit none

integer, intent(in) :: np,deg
real(kind=8), intent(in), dimension(np) :: x,y
real(kind=8), intent(out), dimension(deg+1) :: coeff
integer, intent(out), optional :: status

real(kind=8), allocatable, dimension(:,:) :: mat,bmat
real(kind=8),dimension(:),allocatable :: work,singval
real(kind=8),dimension(1) :: twork
integer :: lwork,info,rank

integer :: i,j

 if (present(status)) status=1
allocate(mat(np,deg+1),bmat(np,1),singval(deg+1))

! Design matrix valorisation
mat=reshape( (/ ((x(i)**(j-1),i=1,np),j=1,deg+1) /),shape=(/ np,deg+1 /) )

! second member valorisation
bmat=reshape ( (/  (y(i),i=1,np)   /)  ,shape = (/ np,1 /))

! query workspace size
call dgelss(np,deg+1,1,mat,np,bmat,np,singval,-1.,rank,twork,-1,info)
lwork=twork(1)
allocate(work(int(lwork)))
! real call
call dgelss(np,deg+1,1,mat,np,bmat,np,singval,-1.,rank,work,lwork,info)

if (info /= 0) then
     if (present(status)) status=0
end if

 coeff = (/ (bmat(i,1),i=1,deg+1) /)

deallocate(work)
deallocate(mat,bmat,singval)
end subroutine

subroutine fvn_s_lspoly_svd(np,x,y,deg,coeff,status)
!
!   Least square polynomial fitting using singular value decomposition
!
!   Find the coefficients of the least square polynomial of a given degree
!   for a set of coordinates.
!
!   The degree must be lower than the number of points
!
!   np (in)             : number of points
!   x(np) (in)          : x data
!   y(np) (in)          : y data
!   deg (in)            : polynomial's degree
!   coeff(deg+1) (out)  : polynomial coefficients
!   status (out)        : =0 if a problem occurs
!
implicit none

integer, intent(in) :: np,deg
real(kind=4), intent(in), dimension(np) :: x,y
real(kind=4), intent(out), dimension(deg+1) :: coeff
integer, intent(out), optional :: status

real(kind=4), allocatable, dimension(:,:) :: mat,bmat
real(kind=4),dimension(:),allocatable :: work,singval
real(kind=4),dimension(1) :: twork
integer :: lwork,info,rank

integer :: i,j

 if (present(status)) status=1
allocate(mat(np,deg+1),bmat(np,1),singval(deg+1))

! Design matrix valorisation
mat=reshape( (/ ((x(i)**(j-1),i=1,np),j=1,deg+1) /),shape=(/ np,deg+1 /) )

! second member valorisation
bmat=reshape ( (/  (y(i),i=1,np)   /)  ,shape = (/ np,1 /))

! query workspace size
call sgelss(np,deg+1,1,mat,np,bmat,np,singval,-1.,rank,twork,-1,info)
lwork=twork(1)
allocate(work(int(lwork)))
! real call
call sgelss(np,deg+1,1,mat,np,bmat,np,singval,-1.,rank,work,lwork,info)

if (info /= 0) then
     if (present(status)) status=0
end if

 coeff = (/ (bmat(i,1),i=1,deg+1) /)

deallocate(work)
deallocate(mat,bmat,singval)
end subroutine


!
! Muller
!
!
!
! William Daniau 2007
!
! This routine is a fortran 90 port of Hans D. Mittelmann's routine muller.f
! http://plato.asu.edu/ftp/other_software/muller.f
!
! it can be used as a replacement for imsl routine dzanly with minor changes
!
!-----------------------------------------------------------------------
!
!   purpose             - zeros of an analytic complex function
!                           using the muller method with deflation
!
!   usage               - call fvn_z_muller (f,eps,eps1,kn,n,nguess,x,itmax,
!                           infer,ier)
!
!   arguments    f      - a complex function subprogram, f(z), written
!                           by the user specifying the equation whose
!                           roots are to be found.  f must appear in
!                           an external statement in the calling pro-
!                           gram.
!                eps    - 1st stopping criterion.  let fp(z)=f(z)/p
!                           where p = (z-z(1))*(z-z(2))*,,,*(z-z(k-1))
!                           and z(1),...,z(k-1) are previously found
!                           roots.  if ((cdabs(f(z)).le.eps) .and.
!                           (cdabs(fp(z)).le.eps)), then z is accepted
!                           as a root. (input)
!                eps1   - 2nd stopping criterion.  a root is accepted
!                           if two successive approximations to a given
!                           root agree within eps1. (input)
!                             note. if either or both of the stopping
!                             criteria are fulfilled, the root is
!                             accepted.
!                kn     - the number of known roots which must be stored
!                           in x(1),...,x(kn), prior to entry to muller
!                nguess - the number of initial guesses provided. these
!                           guesses must be stored in x(kn+1),...,
!                           x(kn+nguess).  nguess must be set equal
!                           to zero if no guesses are provided. (input)
!                n      - the number of new roots to be found by
!                           muller (input)
!                x      - a complex vector of length kn+n.  x(1),...,
!                           x(kn) on input must contain any known
!                           roots.  x(kn+1),..., x(kn+n) on input may,
!                           on user option, contain initial guesses for
!                           the n new roots which are to be computed.
!                           if the user does not provide an initial
!                           guess, zero is used.
!                           on output, x(kn+1),...,x(kn+n) contain the
!                           approximate roots found by muller.
!                itmax  - the maximum allowable number of iterations
!                           per root (input)
!                infer  - an integer vector of length kn+n.  on
!                           output infer(j) contains the number of
!                           iterations used in finding the j-th root
!                           when convergence was achieved.  if
!                           convergence was not obtained in itmax
!                           iterations, infer(j) will be greater than
!                           itmax (output).
!                ier    - error parameter (output)
!                         warning error
!                           ier = 33 indicates failure to converge with-
!                             in itmax iterations for at least one of
!                             the (n) new roots.
!
!
!   remarks      muller always returns the last approximation for root j
!                in x(j). if the convergence criterion is satisfied,
!                then infer(j) is less than or equal to itmax. if the
!                convergence criterion is not satisified, then infer(j)
!                is set to either itmax+1 or itmax+k, with k greater
!                than 1. infer(j) = itmax+1 indicates that muller did
!                not obtain convergence in the allowed number of iter-
!                ations. in this case, the user may wish to set itmax
!                to a larger value. infer(j) = itmax+k means that con-
!                vergence was obtained (on iteration k) for the defla-
!                ted function
!                              fp(z) = f(z)/((z-z(1)...(z-z(j-1)))
!
!                but failed for f(z). in this case, better initial
!                guesses might help or, it might be necessary to relax
!                the convergence criterion.
!
!-----------------------------------------------------------------------
!
subroutine fvn_z_muller (f,eps,eps1,kn,nguess,n,x,itmax,infer,ier)
     implicit none
      double precision :: rzero,rten,rhun,rp01,ax,eps1,qz,eps,tpq
      double complex ::   d,dd,den,fprt,frt,h,rt,t1,t2,t3, &
                          tem,z0,z1,z2,bi,xx,xl,y0,y1,y2,x0, &
                          zero,p1,one,four,p5
      
      double complex, external :: f
      integer :: ickmax,kn,nguess,n,itmax,ier,knp1,knpn,i,l,ic, &
                    knpng,jk,ick,nn,lm1,errcode
      double complex :: x(kn+n)
      integer :: infer(kn+n)
      
      
      data                zero/(0.0,0.0)/,p1/(0.1,0.0)/, &
                          one/(1.0,0.0)/,four/(4.0,0.0)/, &
                          p5/(0.5,0.0)/, &
                          rzero/0.0/,rten/10.0/,rhun/100.0/, &
                          ax/0.1/,ickmax/3/,rp01/0.01/

            ier = 0
            if (n .lt. 1) then ! What the hell are doing here then ...
                return
            end if
            !eps1 = rten **(-nsig)
            eps1 = min(eps1,rp01)
            
            knp1 = kn+1
            knpn = kn+n
            knpng = kn+nguess
            do i=1,knpn
                infer(i) = 0
                if (i .gt. knpng) x(i) = zero
            end do
            l= knp1

            ic=0
rloop:      do while (l<=knpn)   ! Main loop over new roots
                jk = 0
                ick = 0
                xl = x(l)
icloop:         do
                    ic = 0
                    h = ax
                    h = p1*h
                    if (cdabs(xl) .gt. ax) h = p1*xl
!                                  first three points are
!                                    xl+h,  xl-h,  xl
                    rt = xl+h
                    call deflated_work(errcode)
                    if (errcode == 1) then
                        exit icloop
                    end if

                    z0 = fprt
                    y0 = frt
                    x0 = rt
                    rt = xl-h
                    call deflated_work(errcode)
                    if (errcode == 1) then
                        exit icloop
                    end if

                    z1 = fprt
                    y1 = frt
                    h = xl-rt
                    d = h/(rt-x0)
                    rt = xl

                    call deflated_work(errcode)
                    if (errcode == 1) then
                        exit icloop
                    end if

   
                    z2 = fprt
                    y2 = frt
!                                  begin main algorithm
 iloop:             do
                        dd = one + d
                        t1 = z0*d*d
                        t2 = z1*dd*dd
                        xx = z2*dd
                        t3 = z2*d
                        bi = t1-t2+xx+t3
                        den = bi*bi-four*(xx*t1-t3*(t2-xx))
!                                  use denominator of maximum amplitude 
                        t1 = cdsqrt(den)
                        qz = rhun*max(cdabs(bi),cdabs(t1))
                        t2 = bi + t1
                        tpq = cdabs(t2)+qz
                        if (tpq .eq. qz) t2 = zero
                        t3 = bi - t1
                        tpq = cdabs(t3) + qz
                        if (tpq .eq. qz) t3 = zero
                        den = t2
                        qz = cdabs(t3)-cdabs(t2)
                        if (qz .gt. rzero) den = t3
!                                  test for zero denominator            
                        if (cdabs(den) .eq. rzero) then
                            call trans_rt()
                            call deflated_work(errcode)
                            if (errcode == 1) then
                                exit icloop
                            end if
                            z2 = fprt
                            y2 = frt
                            cycle iloop
                        end if


                        d = -xx/den
                        d = d+d
                        h = d*h
                        rt = rt + h
!                                  check convergence of the first kind  
                        if (cdabs(h) .le. eps1*max(cdabs(rt),ax)) then
                            if (ic .ne. 0) then
                                exit icloop
                            end if
                            ic = 1
                            z0 = y1
                            z1 = y2
                            z2 = f(rt)
                            xl = rt
                            ick = ick+1
                            if (ick .le. ickmax) then
                                cycle iloop 
                            end if
!                                  warning error, itmax = maximum
                            jk = itmax + jk
                            ier = 33
                        end if
                        if (ic .ne. 0) then
                            cycle icloop
                        end if
                        call deflated_work(errcode)
                        if (errcode == 1) then
                            exit icloop
                        end if

                        do while ( (cdabs(fprt)-cdabs(z2)*rten) .ge. rzero)
                            !   take remedial action to induce
                            !   convergence
                            d = d*p5
                            h = h*p5
                            rt = rt-h
                            call deflated_work(errcode)
                            if (errcode == 1) then
                                exit icloop
                            end if
                        end do
                        z0 = z1
                        z1 = z2
                        z2 = fprt
                        y0 = y1
                        y1 = y2
                        y2 = frt
                    end do iloop
                end do icloop
        x(l) = rt
        infer(l) = jk
        l = l+1
      end do rloop
      
      contains
        subroutine trans_rt()
           tem = rten*eps1
           if (cdabs(rt) .gt. ax) tem = tem*rt
           rt = rt+tem
           d = (h+tem)*d/h
           h = h+tem
        end subroutine trans_rt
        
        subroutine deflated_work(errcode)
            ! errcode=0 => no errors
            ! errcode=1 => jk>itmax or convergence of second kind achieved
            integer :: errcode,flag
            
            flag=1
    loop1:  do while(flag==1)
                errcode=0
                jk = jk+1
                if (jk .gt. itmax) then
                    ier=33
                    errcode=1
                    return
                end if
                frt = f(rt)
                fprt = frt
                if (l /= 1) then
                    lm1 = l-1
                    do i=1,lm1
                        tem = rt - x(i)
                        if (cdabs(tem) .eq. rzero) then
                        !if (ic .ne. 0) go to 15 !! ?? possible?
                            call trans_rt()
                            cycle loop1
                        end if
                        fprt = fprt/tem
                    end do
                end if
                flag=0
            end do loop1
 
            if (cdabs(fprt) .le. eps .and. cdabs(frt) .le. eps) then
                errcode=1
                return
            end if
            
        end subroutine deflated_work
      
      end subroutine


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
!   Integration
!
!   Only double precision coded atm
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!


subroutine fvn_d_gauss_legendre(n,qx,qw)
!
! This routine compute the n Gauss Legendre abscissas and weights
! Adapted from Numerical Recipes routine gauleg
!
! n (in) : number of points
! qx(out) : abscissas
! qw(out) : weights
!
implicit none
double precision,parameter :: pi=3.141592653589793d0
integer, intent(in) :: n
double precision, intent(out) :: qx(n),qw(n)

integer :: m,i,j
double precision :: z,z1,p1,p2,p3,pp
m=(n+1)/2
do i=1,m
    z=cos(pi*(dble(i)-0.25d0)/(dble(n)+0.5d0))
iloop:  do 
            p1=1.d0
            p2=0.d0
            do j=1,n
                p3=p2
                p2=p1
                p1=((2.d0*dble(j)-1.d0)*z*p2-(dble(j)-1.d0)*p3)/dble(j)
            end do
            pp=dble(n)*(z*p1-p2)/(z*z-1.d0)
            z1=z
            z=z1-p1/pp
            if (dabs(z-z1)<=epsilon(z)) then
                exit iloop
            end if
        end do iloop
    qx(i)=-z
    qx(n+1-i)=z
    qw(i)=2.d0/((1.d0-z*z)*pp*pp)
    qw(n+1-i)=qw(i)
end do
end subroutine



subroutine fvn_d_gl_integ(f,a,b,n,res)
!
! This is a simple non adaptative integration routine 
! using n gauss legendre abscissas and weights
!
!   f(in)   : the function to integrate
!   a(in)   : lower bound
!   b(in)   : higher bound
!   n(in)   : number of gauss legendre pairs
!   res(out): the evaluation of the integral
!
double precision,external :: f
double precision, intent(in) :: a,b
integer, intent(in):: n
double precision, intent(out) :: res

double precision, allocatable :: qx(:),qw(:)
double precision :: xm,xr
integer :: i

! First compute n gauss legendre abs and weight
allocate(qx(n))
allocate(qw(n))
call fvn_d_gauss_legendre(n,qx,qw)

xm=0.5d0*(b+a)
xr=0.5d0*(b-a)

res=0.d0

do i=1,n
    res=res+qw(i)*f(xm+xr*qx(i))
end do

res=xr*res

deallocate(qw)
deallocate(qx)

end subroutine

!!!!!!!!!!!!!!!!!!!!!!!!
!
! Simple and double adaptative Gauss Kronrod integration based on
! a modified version of quadpack ( http://www.netlib.org/quadpack
!
! Common parameters :
!
!       key (in)
!       epsabs
!       epsrel
!
!
!!!!!!!!!!!!!!!!!!!!!!!!

subroutine fvn_d_integ_1_gk(f,a,b,epsabs,epsrel,key,res,abserr,ier,limit)
!
! Evaluate the integral of function f(x) between a and b
!
! f(in) : the function
! a(in) : lower bound
! b(in) : higher bound
! epsabs(in) : desired absolute error
! epsrel(in) : desired relative error
! key(in) : gauss kronrod rule
!                     1:   7 - 15 points
!                     2:  10 - 21 points
!                     3:  15 - 31 points
!                     4:  20 - 41 points
!                     5:  25 - 51 points
!                     6:  30 - 61 points
!
! limit(in) : maximum number of subintervals in the partition of the 
!               given integration interval (a,b). A value of 500 will give the same
!               behaviour as the imsl routine dqdag
!
! res(out) : estimated integral value
! abserr(out) : estimated absolute error
! ier(out) : error flag from quadpack routines
!               0 : no error
!               1 : maximum number of subdivisions allowed
!                   has been achieved. one can allow more
!                   subdivisions by increasing the value of
!                   limit (and taking the according dimension
!                   adjustments into account). however, if
!                   this yield no improvement it is advised
!                   to analyze the integrand in order to
!                   determine the integration difficulaties.
!                   if the position of a local difficulty can
!                   be determined (i.e.singularity,
!                   discontinuity within the interval) one
!                   will probably gain from splitting up the
!                   interval at this point and calling the
!                   integrator on the subranges. if possible,
!                   an appropriate special-purpose integrator
!                   should be used which is designed for
!                   handling the type of difficulty involved.
!               2 : the occurrence of roundoff error is
!                   detected, which prevents the requested
!                   tolerance from being achieved.
!               3 : extremely bad integrand behaviour occurs
!                   at some points of the integration
!                   interval.
!               6 : the input is invalid, because
!                   (epsabs.le.0 and
!                   epsrel.lt.max(50*rel.mach.acc.,0.5d-28))
!                   or limit.lt.1 or lenw.lt.limit*4.
!                   result, abserr, neval, last are set
!                   to zero.
!                   except when lenw is invalid, iwork(1),
!                   work(limit*2+1) and work(limit*3+1) are
!                   set to zero, work(1) is set to a and
!                   work(limit+1) to b.

implicit none
double precision, external :: f
double precision, intent(in) :: a,b,epsabs,epsrel
integer, intent(in) :: key
integer, intent(in) :: limit
double precision, intent(out) :: res,abserr
integer, intent(out) :: ier

double precision, allocatable :: work(:)
integer, allocatable :: iwork(:)
integer :: lenw,neval,last

! imsl value for limit is 500
lenw=limit*4

allocate(iwork(limit))
allocate(work(lenw))

call dqag(f,a,b,epsabs,epsrel,key,res,abserr,neval,ier,limit,lenw,last,iwork,work)

deallocate(work)
deallocate(iwork)

end subroutine



subroutine fvn_d_integ_2_gk(f,a,b,g,h,epsabs,epsrel,key,res,abserr,ier,limit)
!
! Evaluate the double integral of function f(x,y) for x between a and b and y between g(x) and h(x)
!
! f(in) : the function
! a(in) : lower bound
! b(in) : higher bound
! g(in) : external function describing lower bound for y
! h(in) : external function describing higher bound for y
! epsabs(in) : desired absolute error
! epsrel(in) : desired relative error
! key(in) : gauss kronrod rule
!                     1:   7 - 15 points
!                     2:  10 - 21 points
!                     3:  15 - 31 points
!                     4:  20 - 41 points
!                     5:  25 - 51 points
!                     6:  30 - 61 points
!
! limit(in) : maximum number of subintervals in the partition of the 
!               given integration interval (a,b). A value of 500 will give the same
!               behaviour as the imsl routine dqdag
!
! res(out) : estimated integral value
! abserr(out) : estimated absolute error
! ier(out) : error flag from quadpack routines
!               0 : no error
!               1 : maximum number of subdivisions allowed
!                   has been achieved. one can allow more
!                   subdivisions by increasing the value of
!                   limit (and taking the according dimension
!                   adjustments into account). however, if
!                   this yield no improvement it is advised
!                   to analyze the integrand in order to
!                   determine the integration difficulaties.
!                   if the position of a local difficulty can
!                   be determined (i.e.singularity,
!                   discontinuity within the interval) one
!                   will probably gain from splitting up the
!                   interval at this point and calling the
!                   integrator on the subranges. if possible,
!                   an appropriate special-purpose integrator
!                   should be used which is designed for
!                   handling the type of difficulty involved.
!               2 : the occurrence of roundoff error is
!                   detected, which prevents the requested
!                   tolerance from being achieved.
!               3 : extremely bad integrand behaviour occurs
!                   at some points of the integration
!                   interval.
!               6 : the input is invalid, because
!                   (epsabs.le.0 and
!                   epsrel.lt.max(50*rel.mach.acc.,0.5d-28))
!                   or limit.lt.1 or lenw.lt.limit*4.
!                   result, abserr, neval, last are set
!                   to zero.
!                   except when lenw is invalid, iwork(1),
!                   work(limit*2+1) and work(limit*3+1) are
!                   set to zero, work(1) is set to a and
!                   work(limit+1) to b.

implicit none
double precision, external:: f,g,h
double precision, intent(in) :: a,b,epsabs,epsrel
integer, intent(in) :: key,limit
integer, intent(out) :: ier
double precision, intent(out) :: res,abserr


double precision, allocatable :: work(:)
integer, allocatable :: iwork(:)
integer :: lenw,neval,last

! imsl value for limit is 500
lenw=limit*4
allocate(work(lenw))
allocate(iwork(limit))

call dqag_2d_outer(f,a,b,g,h,epsabs,epsrel,key,res,abserr,neval,ier,limit,lenw,last,iwork,work)

deallocate(iwork)
deallocate(work)
end subroutine



subroutine fvn_d_integ_2_inner_gk(f,x,a,b,epsabs,epsrel,key,res,abserr,ier,limit)
!
! Evaluate the single integral of function f(x,y) for y between a and b with a
! given x value
!
! This function is used for the evaluation of the double integral fvn_d_integ_2_gk
!
! f(in) : the function
! x(in) : x
! a(in) : lower bound
! b(in) : higher bound
! epsabs(in) : desired absolute error
! epsrel(in) : desired relative error
! key(in) : gauss kronrod rule
!                     1:   7 - 15 points
!                     2:  10 - 21 points
!                     3:  15 - 31 points
!                     4:  20 - 41 points
!                     5:  25 - 51 points
!                     6:  30 - 61 points
!
! limit(in) : maximum number of subintervals in the partition of the 
!               given integration interval (a,b). A value of 500 will give the same
!               behaviour as the imsl routine dqdag
!
! res(out) : estimated integral value
! abserr(out) : estimated absolute error
! ier(out) : error flag from quadpack routines
!               0 : no error
!               1 : maximum number of subdivisions allowed
!                   has been achieved. one can allow more
!                   subdivisions by increasing the value of
!                   limit (and taking the according dimension
!                   adjustments into account). however, if
!                   this yield no improvement it is advised
!                   to analyze the integrand in order to
!                   determine the integration difficulaties.
!                   if the position of a local difficulty can
!                   be determined (i.e.singularity,
!                   discontinuity within the interval) one
!                   will probably gain from splitting up the
!                   interval at this point and calling the
!                   integrator on the subranges. if possible,
!                   an appropriate special-purpose integrator
!                   should be used which is designed for
!                   handling the type of difficulty involved.
!               2 : the occurrence of roundoff error is
!                   detected, which prevents the requested
!                   tolerance from being achieved.
!               3 : extremely bad integrand behaviour occurs
!                   at some points of the integration
!                   interval.
!               6 : the input is invalid, because
!                   (epsabs.le.0 and
!                   epsrel.lt.max(50*rel.mach.acc.,0.5d-28))
!                   or limit.lt.1 or lenw.lt.limit*4.
!                   result, abserr, neval, last are set
!                   to zero.
!                   except when lenw is invalid, iwork(1),
!                   work(limit*2+1) and work(limit*3+1) are
!                   set to zero, work(1) is set to a and
!                   work(limit+1) to b.

implicit none
double precision, external:: f
double precision, intent(in) :: x,a,b,epsabs,epsrel
integer, intent(in) :: key,limit
integer, intent(out) :: ier
double precision, intent(out) :: res,abserr


double precision, allocatable :: work(:)
integer, allocatable :: iwork(:)
integer :: lenw,neval,last

! imsl value for limit is 500
lenw=limit*4
allocate(work(lenw))
allocate(iwork(limit))

call dqag_2d_inner(f,x,a,b,epsabs,epsrel,key,res,abserr,neval,ier,limit,lenw,last,iwork,work)

deallocate(iwork)
deallocate(work)
end subroutine

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Include the modified quadpack files
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
include "fvn_quadpack/dqag_2d_inner.f"
include "fvn_quadpack/dqk15_2d_inner.f"
include "fvn_quadpack/dqk31_2d_outer.f"
include "fvn_quadpack/d1mach.f"
include "fvn_quadpack/dqk31_2d_inner.f"
include "fvn_quadpack/dqage.f"
include "fvn_quadpack/dqk15.f"
include "fvn_quadpack/dqk21.f"
include "fvn_quadpack/dqk31.f"
include "fvn_quadpack/dqk41.f"
include "fvn_quadpack/dqk51.f"
include "fvn_quadpack/dqk61.f"
include "fvn_quadpack/dqk41_2d_outer.f"
include "fvn_quadpack/dqk41_2d_inner.f"
include "fvn_quadpack/dqag_2d_outer.f"
include "fvn_quadpack/dqpsrt.f"
include "fvn_quadpack/dqag.f"
include "fvn_quadpack/dqage_2d_outer.f"
include "fvn_quadpack/dqage_2d_inner.f"
include "fvn_quadpack/dqk51_2d_outer.f"
include "fvn_quadpack/dqk51_2d_inner.f"
include "fvn_quadpack/dqk61_2d_outer.f"
include "fvn_quadpack/dqk21_2d_outer.f"
include "fvn_quadpack/dqk61_2d_inner.f"
include "fvn_quadpack/dqk21_2d_inner.f"
include "fvn_quadpack/dqk15_2d_outer.f"

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! Trigonometric functions
!
! fvn_z_acos, fvn_z_asin : complex arc cosine and sine
! fvn_d_acosh : arc cosinus hyperbolic
! fvn_d_asinh : arc sinus hyperbolic
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! February 2008 
! All Trigonometric functions removed due to implementation of fnlib

function fvn_z_acos(z)
    ! double complex arccos function adapted from 
    ! the c gsl library
    ! http://www.gnu.org/software/gsl/
    implicit none
    complex(kind=8) :: fvn_z_acos
    complex(kind=8) :: z
    real(kind=8) :: rz,iz,x,y,a,b,y2,r,s,d,apx,am1
    real(kind=8),parameter :: a_crossover=1.5_8,b_crossover = 0.6417_8
    complex(kind=8),parameter :: i=(0._8,1._8)
    real(kind=8) :: r_res,i_res

    rz=dreal(z)
    iz=dimag(z)
    if ( iz == 0._8 ) then
        fvn_z_acos=fvn_z_acos_real(rz)
        return
    end if

    x=dabs(rz)
    y=dabs(iz)
    r=fvn_d_hypot(x+1.,y)
    s=fvn_d_hypot(x-1.,y)
    a=0.5*(r + s)
    b=x/a
    y2=y*y

    if (b <= b_crossover) then
        r_res=dacos(b)
    else
        if (x <= 1.) then
            d=0.5*(a+x)*(y2/(r+x+1)+(s + (1 - x)))
            r_res=datan(dsqrt(d)/x)
        else
            apx=a+x
            d=0.5*(apx/(r+x+1)+apx/(s + (x - 1)))
            r_res=datan((y*dsqrt(d))/x);
        end if
    end if

    if (a <= a_crossover) then
        if (x < 1.) then
            am1=0.5*(y2 / (r + (x + 1)) + y2 / (s + (1 - x)))
        else
            am1=0.5*(y2 / (r + (x + 1)) + (s + (x - 1)))
        end if
        i_res = dlog(1.+(am1 + sqrt (am1 * (a + 1))));
    else
        i_res = dlog (a + dsqrt (a*a - 1.));
    end if
    if (rz <0.) then
        r_res=fvn_pi-r_res
    end if
    i_res=-sign(1._8,iz)*i_res
    fvn_z_acos=dcmplx(r_res)+fvn_i*dcmplx(i_res)

end function fvn_z_acos

function fvn_z_acos_real(r)
    ! return the double complex arc cosinus for a 
    ! double precision argument
    implicit none
    real(kind=8) :: r
    complex(kind=8) :: fvn_z_acos_real

    if (dabs(r)<=1._8) then
        fvn_z_acos_real=dcmplx(dacos(r))
        return
    end if
    if (r < 0._8) then
        fvn_z_acos_real=dcmplx(fvn_pi)-fvn_i*dcmplx(fvn_d_acosh(-r))
    else
        fvn_z_acos_real=fvn_i*dcmplx(fvn_d_acosh(r))
    end if
end function


function fvn_z_asin(z)
    ! double complex arcsin function derived from 
    ! the c gsl library
    ! http://www.gnu.org/software/gsl/
    implicit none
    complex(kind=8) :: fvn_z_asin
    complex(kind=8) :: z
    real(kind=8) :: rz,iz,x,y,a,b,y2,r,s,d,apx,am1
    real(kind=8),parameter :: a_crossover=1.5_8,b_crossover = 0.6417_8
    real(kind=8) :: r_res,i_res

    rz=dreal(z)
    iz=dimag(z)
    if ( iz == 0._8 ) then
        ! z is real
        fvn_z_asin=fvn_z_asin_real(rz)
        return
    end if
    
    x=dabs(rz)
    y=dabs(iz)
    r=fvn_d_hypot(x+1.,y)
    s=fvn_d_hypot(x-1.,y)
    a=0.5*(r + s)
    b=x/a
    y2=y*y

    if (b <= b_crossover) then
        r_res=dasin(b)
    else
        if (x <= 1.) then
            d=0.5*(a+x)*(y2/(r+x+1)+(s + (1 - x)))
            r_res=datan(x/dsqrt(d))
        else
            apx=a+x
            d=0.5*(apx/(r+x+1)+apx/(s + (x - 1)))
            r_res=datan(x/(y*dsqrt(d)));
        end if
    end if

    if (a <= a_crossover) then
        if (x < 1.) then
            am1=0.5*(y2 / (r + (x + 1)) + y2 / (s + (1 - x)))
        else
            am1=0.5*(y2 / (r + (x + 1)) + (s + (x - 1)))
        end if
        i_res = dlog(1.+(am1 + sqrt (am1 * (a + 1))));
    else
        i_res = dlog (a + dsqrt (a*a - 1.));
    end if
    r_res=sign(1._8,rz)*r_res
    i_res=sign(1._8,iz)*i_res
    fvn_z_asin=dcmplx(r_res)+fvn_i*dcmplx(i_res)

end function fvn_z_asin

function fvn_z_asin_real(r)
    ! return the double complex arc sinus for a 
    ! double precision argument
    implicit none
    real(kind=8) :: r
    complex(kind=8) :: fvn_z_asin_real

    if (dabs(r)<=1._8) then
        fvn_z_asin_real=dcmplx(dasin(r))
        return
    end if
    if (r < 0._8) then
        fvn_z_asin_real=dcmplx(-fvn_pi/2._8)+fvn_i*dcmplx(fvn_d_acosh(-r))
    else
        fvn_z_asin_real=dcmplx(fvn_pi/2._8)-fvn_i*dcmplx(fvn_d_acosh(r))
    end if
end function fvn_z_asin_real

function fvn_d_acosh(r)
    ! return the arc hyperbolic cosine
    implicit none
    real(kind=8) :: r
    real(kind=8) :: fvn_d_acosh
    if (r >=1) then
        fvn_d_acosh=dlog(r+dsqrt(r*r-1))
    else
        !! TODO : Better error handling!!!!!!
        stop "Argument to fvn_d_acosh lesser than 1"
    end if
end function fvn_d_acosh

function fvn_d_asinh(r)
    ! return the arc hyperbolic sine
    implicit none
    real(kind=8) :: r
    real(kind=8) :: fvn_d_asinh
    fvn_d_asinh=dlog(r+dsqrt(r*r+1))
end function fvn_d_asinh

function fvn_d_hypot(a,b)
    implicit none
    ! return the euclidian norm of vector(a,b)
    real(kind=8) :: a,b
    real(kind=8) :: fvn_d_hypot
    fvn_d_hypot=dsqrt(a*a+b*b)
end function

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! SPARSE RESOLUTION
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! 
! Sparse resolution is done by interfaçing Tim Davi's UMFPACK
! http://www.cise.ufl.edu/research/sparse/SuiteSparse/
! Used packages from SuiteSparse : AMD,UMFPACK,UFconfig
!
! Solve Ax=B using UMFPACK
!
! Where A is a sparse matrix given in its triplet form 
! T -> non zero elements
! Ti,Tj -> row and column index (1-based) of the given elt
! n : rank of matrix A
! nz : number of non zero elts
!
! fvn_*_sparse_solve
! * = zl : double complex + integer(8)
! * = zi : double complex + integer(4)
!
subroutine fvn_zl_sparse_solve(n,nz,T,Ti,Tj,B,x,status)
implicit none
integer(8), intent(in) :: n,nz
complex(8),dimension(nz),intent(in) :: T
integer(8),dimension(nz),intent(in)  :: Ti,Tj
complex(8),dimension(n),intent(in) :: B
complex(8),dimension(n),intent(out) :: x
integer(8), intent(out) :: status

integer(8),dimension(:),allocatable :: wTi,wTj
real(8),dimension(:),allocatable :: Tx,Tz
real(8),dimension(:),allocatable :: Ax,Az
integer(8),dimension(:),allocatable :: Ap,Ai
integer(8) :: symbolic,numeric
real(8),dimension(:),allocatable :: xx,xz,bx,bz
real(8),dimension(90) :: info
real(8),dimension(20) :: control
integer(8) :: sys


status=0

! we use a working copy of Ti and Tj to perform 1-based to 0-based translation
! Tx and Tz are the real and imaginary parts of T
allocate(wTi(nz),wTj(nz))
allocate(Tx(nz),Tz(nz))
Tx=dble(T)
Tz=aimag(T)
wTi=Ti-1
wTj=Tj-1
allocate(Ax(nz),Az(nz))
allocate(Ap(n+1),Ai(nz))

! perform the triplet to compressed column form -> Ap,Ai,Ax,Az
call umfpack_zl_triplet_to_col(n,n,nz,wTi,wTj,Tx,Tz,Ap,Ai,Ax,Az,status)
! if status is not zero a problem has occured
if (status /= 0) then
    write(*,*) "Problem during umfpack_zl_triplet_to_col"
endif

! Define defaults control values
call umfpack_zl_defaults(control)

! Symbolic analysis
call umfpack_zl_symbolic(n,n,Ap,Ai,Ax,Az,symbolic, control, info)
! info(1) should be zero
if (info(1) /= 0) then
    write(*,*) "Problem during symbolic analysis"
    status=info(1)
endif

! Numerical factorization
call umfpack_zl_numeric (Ap, Ai, Ax, Az, symbolic, numeric, control, info)
! info(1) should be zero
if (info(1) /= 0) then
    write(*,*) "Problem during numerical factorization"
    status=info(1)
endif

! free the C symbolic pointer
call umfpack_zl_free_symbolic (symbolic)

allocate(bx(n),bz(n),xx(n),xz(n))
bx=dble(B)
bz=aimag(B)
sys=0
! sys may be used to define type of solving -> see umfpack.h

! Solving
call umfpack_zl_solve (sys, Ap, Ai, Ax,Az, xx,xz, bx,bz, numeric, control, info)
! info(1) should be zero
if (info(1) /= 0) then
    write(*,*) "Problem during solving"
    status=info(1)
endif


! free the C numeric pointer
call umfpack_zl_free_numeric (numeric)

x=dcmplx(xx,xz)

deallocate(bx,bz,xx,xz)
deallocate(Ax,Az)
deallocate(Tx,Tz)
deallocate(wTi,wTj)
end subroutine





subroutine fvn_zi_sparse_solve(n,nz,T,Ti,Tj,B,x,status)
implicit none
integer(4), intent(in) :: n,nz
complex(8),dimension(nz),intent(in) :: T
integer(4),dimension(nz),intent(in)  :: Ti,Tj
complex(8),dimension(n),intent(in) :: B
complex(8),dimension(n),intent(out) :: x
integer(4), intent(out) :: status

integer(4),dimension(:),allocatable :: wTi,wTj
real(8),dimension(:),allocatable :: Tx,Tz
real(8),dimension(:),allocatable :: Ax,Az
integer(4),dimension(:),allocatable :: Ap,Ai
!integer(8) :: symbolic,numeric
integer(4),dimension(2) :: symbolic,numeric
! As symbolic and numeric are used to store a C pointer, it is necessary to
! still use an integer(8) for 64bits machines
! An other possibility : integer(4),dimension(2) :: symbolic,numeric
real(8),dimension(:),allocatable :: xx,xz,bx,bz
real(8),dimension(90) :: info
real(8),dimension(20) :: control
integer(4) :: sys

status=0
! we use a working copy of Ti and Tj to perform 1-based to 0-based translation
! Tx and Tz are the real and imaginary parts of T
allocate(wTi(nz),wTj(nz))
allocate(Tx(nz),Tz(nz))
Tx=dble(T)
Tz=aimag(T)
wTi=Ti-1
wTj=Tj-1
allocate(Ax(nz),Az(nz))
allocate(Ap(n+1),Ai(nz))

! perform the triplet to compressed column form -> Ap,Ai,Ax,Az
call umfpack_zi_triplet_to_col(n,n,nz,wTi,wTj,Tx,Tz,Ap,Ai,Ax,Az,status)
! if status is not zero a problem has occured
if (status /= 0) then
    write(*,*) "Problem during umfpack_zl_triplet_to_col"
endif

! Define defaults control values
call umfpack_zi_defaults(control)

! Symbolic analysis
call umfpack_zi_symbolic(n,n,Ap,Ai,Ax,Az,symbolic, control, info)
! info(1) should be zero
if (info(1) /= 0) then
    write(*,*) "Problem during symbolic analysis"
    status=info(1)
endif

! Numerical factorization
call umfpack_zi_numeric (Ap, Ai, Ax, Az, symbolic, numeric, control, info)
! info(1) should be zero
if (info(1) /= 0) then
    write(*,*) "Problem during numerical factorization"
    status=info(1)
endif

! free the C symbolic pointer
call umfpack_zi_free_symbolic (symbolic)

allocate(bx(n),bz(n),xx(n),xz(n))
bx=dble(B)
bz=aimag(B)
sys=0
! sys may be used to define type of solving -> see umfpack.h

! Solving
call umfpack_zi_solve (sys, Ap, Ai, Ax,Az, xx,xz, bx,bz, numeric, control, info)
! info(1) should be zero
if (info(1) /= 0) then
    write(*,*) "Problem during solving"
    status=info(1)
endif

! free the C numeric pointer
call umfpack_zi_free_numeric (numeric)

x=dcmplx(xx,xz)

deallocate(bx,bz,xx,xz)
deallocate(Ax,Az)
deallocate(Tx,Tz)
deallocate(wTi,wTj)
end subroutine






subroutine fvn_dl_sparse_solve(n,nz,T,Ti,Tj,B,x,status)
implicit none
integer(8), intent(in) :: n,nz
real(8),dimension(nz),intent(in) :: T
integer(8),dimension(nz),intent(in)  :: Ti,Tj
real(8),dimension(n),intent(in) :: B
real(8),dimension(n),intent(out) :: x
integer(8), intent(out) :: status

integer(8),dimension(:),allocatable :: wTi,wTj
real(8),dimension(:),allocatable :: A
integer(8),dimension(:),allocatable :: Ap,Ai
!integer(8) :: symbolic,numeric
integer(8) :: symbolic,numeric
real(8),dimension(90) :: info
real(8),dimension(20) :: control
integer(8) :: sys

status=0
! we use a working copy of Ti and Tj to perform 1-based to 0-based translation
allocate(wTi(nz),wTj(nz))
wTi=Ti-1
wTj=Tj-1
allocate(A(nz))
allocate(Ap(n+1),Ai(nz))

! perform the triplet to compressed column form -> Ap,Ai,Ax,Az
call umfpack_dl_triplet_to_col(n,n,nz,wTi,wTj,T,Ap,Ai,A,status)
! if status is not zero a problem has occured
if (status /= 0) then
    write(*,*) "Problem during umfpack_dl_triplet_to_col"
endif

! Define defaults control values
call umfpack_dl_defaults(control)

! Symbolic analysis
call umfpack_dl_symbolic(n,n,Ap,Ai,A,symbolic, control, info)
! info(1) should be zero
if (info(1) /= 0) then
    write(*,*) "Problem during symbolic analysis"
    status=info(1)
endif

! Numerical factorization
call umfpack_dl_numeric (Ap, Ai, A, symbolic, numeric, control, info)
! info(1) should be zero
if (info(1) /= 0) then
    write(*,*) "Problem during numerical factorization"
    status=info(1)
endif

! free the C symbolic pointer
call umfpack_dl_free_symbolic (symbolic)

sys=0
! sys may be used to define type of solving -> see umfpack.h

! Solving
call umfpack_dl_solve (sys, Ap, Ai, A, x, B, numeric, control, info)
! info(1) should be zero
if (info(1) /= 0) then
    write(*,*) "Problem during solving"
    status=info(1)
endif

! free the C numeric pointer
call umfpack_dl_free_numeric (numeric)

deallocate(A)
deallocate(wTi,wTj)
end subroutine






subroutine fvn_di_sparse_solve(n,nz,T,Ti,Tj,B,x,status)
implicit none
integer(4), intent(in) :: n,nz
real(8),dimension(nz),intent(in) :: T
integer(4),dimension(nz),intent(in)  :: Ti,Tj
real(8),dimension(n),intent(in) :: B
real(8),dimension(n),intent(out) :: x
integer(4), intent(out) :: status

integer(4),dimension(:),allocatable :: wTi,wTj
real(8),dimension(:),allocatable :: A
integer(4),dimension(:),allocatable :: Ap,Ai
!integer(8) :: symbolic,numeric
integer(4),dimension(2) :: symbolic,numeric
! As symbolic and numeric are used to store a C pointer, it is necessary to
! still use an integer(8) for 64bits machines
! An other possibility : integer(4),dimension(2) :: symbolic,numeric
real(8),dimension(90) :: info
real(8),dimension(20) :: control
integer(4) :: sys

status=0
! we use a working copy of Ti and Tj to perform 1-based to 0-based translation
allocate(wTi(nz),wTj(nz))
wTi=Ti-1
wTj=Tj-1
allocate(A(nz))
allocate(Ap(n+1),Ai(nz))

! perform the triplet to compressed column form -> Ap,Ai,Ax,Az
call umfpack_di_triplet_to_col(n,n,nz,wTi,wTj,T,Ap,Ai,A,status)
! if status is not zero a problem has occured
if (status /= 0) then
    write(*,*) "Problem during umfpack_di_triplet_to_col"
endif

! Define defaults control values
call umfpack_di_defaults(control)

! Symbolic analysis
call umfpack_di_symbolic(n,n,Ap,Ai,A,symbolic, control, info)
! info(1) should be zero
if (info(1) /= 0) then
    write(*,*) "Problem during symbolic analysis"
    status=info(1)
endif

! Numerical factorization
call umfpack_di_numeric (Ap, Ai, A, symbolic, numeric, control, info)
! info(1) should be zero
if (info(1) /= 0) then
    write(*,*) "Problem during numerical factorization"
    status=info(1)
endif

! free the C symbolic pointer
call umfpack_di_free_symbolic (symbolic)

sys=0
! sys may be used to define type of solving -> see umfpack.h

! Solving
call umfpack_di_solve (sys, Ap, Ai, A, x, B, numeric, control, info)
! info(1) should be zero
if (info(1) /= 0) then
    write(*,*) "Problem during solving"
    status=info(1)
endif

! free the C numeric pointer
call umfpack_di_free_numeric (numeric)

deallocate(A)
deallocate(wTi,wTj)
end subroutine

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! Special Functions
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! February 2008 
! All Special functions removed due to implementation of fnlib
!
function fvn_d_lngamma(x)
      ! This function returns ln(gamma(x))
      ! adapted from Numerical Recipes
      implicit none
      real(kind=8) :: x
      real(kind=8) :: fvn_d_lngamma

      real(kind=8) :: ser,stp,tmp,y,cof(6)
      integer(kind=4) :: i

      cof = (/ 76.18009172947146d0,-86.50532032941677d0, &
                  24.01409824083091d0,-1.231739572450155d0, &
                  .1208650973866179d-2,-.5395239384953d-5 /)
      stp = 2.5066282746310005d0

      tmp=x+5.5d0
      tmp=(x+0.5d0)*log(tmp)-tmp

      ser=1.000000000190015d0

      y=x

      do i=1,6
            y=y+1.d0
            ser=ser+cof(i)/y
      end do
      fvn_d_lngamma=tmp+log(stp*ser/x)
end function

function fvn_d_factorial(n)
      ! This function returns factorial(n) as a real(8)
      ! adapted from Numerical Recipes
      ! real value is calculated for integers lower than 32
      implicit none
      integer(kind=4) :: n
      real(kind=8) :: fvn_d_factorial

      integer(kind=4) :: j

      fvn_d_factorial=1.

      if (n < 0) then
            write(*,*) "Factorial of a negative integer"
            stop
      end if

      if (n == 0) then
            return
      end if

      if (n <= 32) then 
            do j=1,n
                  fvn_d_factorial=fvn_d_factorial*j
            end do
            return
      else
            fvn_d_factorial=exp(fvn_d_lngamma(dble(n)+1.))
            return
      end if
end function

function fvn_d_csevl(x,a,n)
      implicit none
      ! This function evaluate the n-term chebyshev series a at x
      ! directly adapted from http://www.netlib.org/fn
      real(kind=8) :: x
      real(kind=8), dimension(n) :: a
      integer(kind=4) :: n
      real(kind=8) :: fvn_d_csevl

      real(kind=8) :: twox, b0, b1, b2
      integer(kind=4) :: i,ni

      twox = 2.0d0*x
      b1 = 0.d0
      b0 = 0.d0
      do i=1,n
        b2 = b1
        b1 = b0
        ni = n - i + 1
        b0 = twox*b1 - b2 + a(ni)
      end do

      fvn_d_csevl = 0.5d0 * (b0-b2)

end function

function fvn_s_csevl(x,a,n)
      implicit none
      ! This function evaluate the n-term chebyshev series a at x
      ! directly adapted from http://www.netlib.org/fn
      real(kind=4) :: x
      real(kind=4), dimension(n) :: a
      integer(kind=4) :: n
      real(kind=4) :: fvn_s_csevl

      real(kind=4) :: twox, b0, b1, b2
      integer(kind=4) :: i,ni

      twox = 2.0d0*x
      b1 = 0.d0
      b0 = 0.d0
      do i=1,n
        b2 = b1
        b1 = b0
        ni = n - i + 1
        b0 = twox*b1 - b2 + a(ni)
      end do

      fvn_s_csevl = 0.5d0 * (b0-b2)

end function



end module fvn