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
! 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 ;-) Nevertheless, you
! may give credits to quadpack authors. 
!
! svn version
! 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
            
            
contains 

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! 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
    !
    integer, intent(in) :: d
    real, intent(in) :: a(d,d)
    real, intent(out) :: inva(d,d)
    integer, intent(out) :: status

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

    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
        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
        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
    !
    integer, intent(in) :: d
    double precision, intent(in) :: a(d,d)
    double precision, intent(out) :: inva(d,d)
    integer, intent(out) :: status

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

    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
        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
        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
    !
    integer, intent(in) :: d
    complex, intent(in) :: a(d,d)
    complex, intent(out) :: inva(d,d)
    integer, intent(out) :: status

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

    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
        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
        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
    !
    integer, intent(in) :: d
    double complex, intent(in) :: a(d,d)
    double complex, intent(out) :: inva(d,d)
    integer, intent(out) :: status

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

    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
        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
        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
    !
    integer, intent(in) :: d
    real, intent(in) :: a(d,d)
    integer, intent(out) :: status
    real :: fvn_s_det
    
    real, allocatable :: wc_a(:,:)
    integer, allocatable :: ipiv(:)
    integer :: info,i
    
    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
        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
    !
    integer, intent(in) :: d
    double precision, intent(in) :: a(d,d)
    integer, intent(out) :: status
    double precision :: fvn_d_det
    
    double precision, allocatable :: wc_a(:,:)
    integer, allocatable :: ipiv(:)
    integer :: info,i
    
    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
        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
    !
    integer, intent(in) :: d
    complex, intent(in) :: a(d,d)
    integer, intent(out) :: status
    complex :: fvn_c_det
    
    complex, allocatable :: wc_a(:,:)
    integer, allocatable :: ipiv(:)
    integer :: info,i
    
    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
        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
    !
    integer, intent(in) :: d
    double complex, intent(in) :: a(d,d)
    integer, intent(out) :: status
    double complex :: fvn_z_det
    
    double complex, allocatable :: wc_a(:,:)
    integer, allocatable :: ipiv(:)
    integer :: info,i
    
    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
        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
    !
    integer, intent(in) :: d
    real, intent(in) :: a(d,d)
    real, intent(out) :: rcond
    integer, intent(out) :: 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
    
    
    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
        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
        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
    !
    integer, intent(in) :: d
    double precision, intent(in) :: a(d,d)
    double precision, intent(out) :: rcond
    integer, intent(out) :: 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
    
    
    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
        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
        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
    !
    integer, intent(in) :: d
    complex, intent(in) :: a(d,d)
    real, intent(out) :: rcond
    integer, intent(out) :: 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
    
    
    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
        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
        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
    !
    integer, intent(in) :: d
    double complex, intent(in) :: a(d,d)
    double precision, intent(out) :: rcond
    integer, intent(out) :: 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
    
    
    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
        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
        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

    integer, intent(in) :: d
    real, intent(in) :: a(d,d)
    complex, intent(out) :: evala(d)
    complex, intent(out) :: eveca(d,d)
    integer, intent(out) :: 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
    
    ! 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
        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
    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) :: 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
    
    ! 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
        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

    integer, intent(in) :: d
    complex, intent(in) :: a(d,d)
    complex, intent(out) :: evala(d)
    complex, intent(out) :: eveca(d,d)
    integer, intent(out) :: 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

    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
        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

    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) :: 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
    
    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
        status=0
    end if
    deallocate(rwork)
    deallocate(work)
    deallocate(wc_a)

end subroutine

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! 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


!
! 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
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
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=dble(z)
    iz=aimag(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=dble(z)
    iz=aimag(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




end module fvn