module fvn_sparse
use fvn_common
implicit none

! 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

interface fvn_sparse_det
    module procedure fvn_zl_sparse_det,fvn_zi_sparse_det,fvn_dl_sparse_det,fvn_di_sparse_det
end interface fvn_sparse_det
contains

function umfpack_return_code(c)
  implicit none
  integer(kind=sp_kind), intent(in) :: c
  character(len=80) :: umfpack_return_code
  select case(c)
    case(0)
      umfpack_return_code="UMFPACK_OK"
    case(1)
      umfpack_return_code="UMFPACK_WARNING_singular_matrix"
    case(2)
      umfpack_return_code="UMFPACK_WARNING_determinant_underflow"
    case(3)
      umfpack_return_code="UMFPACK_WARNING_determinant_overflow"
    case(-1)
      umfpack_return_code="UMFPACK_ERROR_out_of_memory"
    case(-3)
      umfpack_return_code="UMFPACK_ERROR_invalid_Numeric_object"
    case(-4)
      umfpack_return_code="UMFPACK_ERROR_invalid_Symbolic_object"
    case(-5)
      umfpack_return_code="UMFPACK_ERROR_argument_missing"
    case(-6)
      umfpack_return_code="UMFPACK_ERROR_n_nonpositive"
    case(-8)
      umfpack_return_code="UMFPACK_ERROR_invalid_matrix"
    case(-11)
      umfpack_return_code="UMFPACK_ERROR_different_pattern"
    case(-13)
      umfpack_return_code="UMFPACK_ERROR_invalid_system"
    case(-15)
      umfpack_return_code="UMFPACK_ERROR_invalid_permutation"
    case(-911)
      umfpack_return_code="UMFPACK_ERROR_internal_error"
    case(-17)
      umfpack_return_code="UMFPACK_ERROR_file_IO"
    case default
      umfpack_return_code="Unknown return code"
  end select
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(kind=dp_kind)
! * = zi : double complex + integer(kind=sp_kind)
!
subroutine fvn_zl_sparse_solve(n,nz,T,Ti,Tj,B,x,status,det)
implicit none
integer(kind=dp_kind), intent(in) :: n,nz
complex(kind=dp_kind),dimension(nz),intent(in) :: T
integer(kind=dp_kind),dimension(nz),intent(in)  :: Ti,Tj
complex(kind=dp_kind),dimension(n),intent(in) :: B
complex(kind=dp_kind),dimension(n),intent(out) :: x
integer(kind=dp_kind), intent(out) :: status
real(kind=dp_kind), dimension(3), optional, intent(out) :: det

integer(kind=dp_kind),dimension(:),allocatable :: wTi,wTj
real(kind=dp_kind),dimension(:),allocatable :: Tx,Tz
real(kind=dp_kind),dimension(:),allocatable :: Ax,Az
integer(kind=dp_kind),dimension(:),allocatable :: Ap,Ai
integer(kind=dp_kind) :: symbolic,numeric
real(kind=dp_kind),dimension(:),allocatable :: xx,xz,bx,bz
real(kind=dp_kind),dimension(90) :: info
real(kind=dp_kind),dimension(20) :: control
integer(kind=dp_kind) :: 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 : ",trim(umfpack_return_code(int(status,kind=sp_kind)))
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 : ",trim(umfpack_return_code(int(info(1),kind=sp_kind)))
    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 : ",trim(umfpack_return_code(int(info(1),kind=sp_kind)))
    status=info(1)
endif

! free the C symbolic pointer
call umfpack_zl_free_symbolic (symbolic)

! if parameter det is present, the determinant of the matrix is calculated
if (present(det) ) then
  call umfpack_zl_get_determinant(det(1),det(2),det(3),numeric,info,status)  
  ! info(1) should be zero
  if (info(1) /= 0) then
      if ( (info(1) < 1) .or. (info(1) >3) ) then ! not a warning
          write(*,*) "Problem during sparse determinant : ",trim(umfpack_return_code(int(info(1),kind=sp_kind)))
      endif
      status=info(1)
  endif
endif



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 : ",trim(umfpack_return_code(int(info(1),kind=sp_kind)))
    status=info(1)
endif


! free the C numeric pointer
call umfpack_zl_free_numeric (numeric)

x=cmplx(xx,xz,dp_kind)

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,det)
implicit none
integer(kind=sp_kind), intent(in) :: n,nz
complex(kind=dp_kind),dimension(nz),intent(in) :: T
integer(kind=sp_kind),dimension(nz),intent(in)  :: Ti,Tj
complex(kind=dp_kind),dimension(n),intent(in) :: B
complex(kind=dp_kind),dimension(n),intent(out) :: x
integer(kind=sp_kind), intent(out) :: status
real(kind=dp_kind), dimension(3), optional, intent(out) :: det

integer(kind=sp_kind),dimension(:),allocatable :: wTi,wTj
real(kind=dp_kind),dimension(:),allocatable :: Tx,Tz
real(kind=dp_kind),dimension(:),allocatable :: Ax,Az
integer(kind=sp_kind),dimension(:),allocatable :: Ap,Ai
!integer(kind=dp_kind) :: symbolic,numeric
integer(kind=sp_kind),dimension(2) :: symbolic,numeric
! As symbolic and numeric are used to store a C pointer, it is necessary to
! still use an integer(kind=dp_kind) for 64bits machines
! An other possibility : integer(kind=sp_kind),dimension(2) :: symbolic,numeric
real(kind=dp_kind),dimension(:),allocatable :: xx,xz,bx,bz
real(kind=dp_kind),dimension(90) :: info
real(kind=dp_kind),dimension(20) :: control
integer(kind=sp_kind) :: 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 : ",trim(umfpack_return_code(status))
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 : ",trim(umfpack_return_code(int(info(1),kind=sp_kind)))
    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 : ",trim(umfpack_return_code(int(info(1),kind=sp_kind)))
    status=info(1)
endif

! free the C symbolic pointer
call umfpack_zi_free_symbolic (symbolic)

! if parameter det is present, the determinant of the matrix is calculated
if (present(det) ) then
  call umfpack_zi_get_determinant(det(1),det(2),det(3),numeric,info,status)
  ! info(1) should be zero
  if (info(1) /= 0) then
      if ( (info(1) < 1) .or. (info(1) >3) ) then ! not a warning
          write(*,*) "Problem during sparse determinant : ",trim(umfpack_return_code(int(info(1),kind=sp_kind)))
      endif
      status=info(1)
  endif
endif




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 : ",trim(umfpack_return_code(int(info(1),kind=sp_kind)))
    status=info(1)
endif

! free the C numeric pointer
call umfpack_zi_free_numeric (numeric)

x=cmplx(xx,xz,dp_kind)

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,det)
implicit none
integer(kind=dp_kind), intent(in) :: n,nz
real(kind=dp_kind),dimension(nz),intent(in) :: T
integer(kind=dp_kind),dimension(nz),intent(in)  :: Ti,Tj
real(kind=dp_kind),dimension(n),intent(in) :: B
real(kind=dp_kind),dimension(n),intent(out) :: x
integer(kind=dp_kind), intent(out) :: status
real(kind=dp_kind), dimension(2), optional, intent(out) :: det

integer(kind=dp_kind),dimension(:),allocatable :: wTi,wTj
real(kind=dp_kind),dimension(:),allocatable :: A
integer(kind=dp_kind),dimension(:),allocatable :: Ap,Ai
!integer(kind=dp_kind) :: symbolic,numeric
integer(kind=dp_kind) :: symbolic,numeric
real(kind=dp_kind),dimension(90) :: info
real(kind=dp_kind),dimension(20) :: control
integer(kind=dp_kind) :: 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 : ",trim(umfpack_return_code(int(status,kind=sp_kind)))
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 : ",trim(umfpack_return_code(int(info(1),kind=sp_kind)))
    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 : ",trim(umfpack_return_code(int(info(1),kind=sp_kind)))
    status=info(1)
endif

! free the C symbolic pointer
call umfpack_dl_free_symbolic (symbolic)

! if parameter det is present, the determinant of the matrix is calculated
if (present(det) ) then
  call umfpack_dl_get_determinant(det(1),det(2),numeric,info,status)
  ! info(1) should be zero
  if (info(1) /= 0) then
      if ( (info(1) < 1) .or. (info(1) >3) ) then ! not a warning
          write(*,*) "Problem during sparse determinant : ",trim(umfpack_return_code(int(info(1),kind=sp_kind)))
      endif
      status=info(1)
  endif
endif


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 : ",trim(umfpack_return_code(int(info(1),kind=sp_kind)))
    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,det)
implicit none
integer(kind=sp_kind), intent(in) :: n,nz
real(kind=dp_kind),dimension(nz),intent(in) :: T
integer(kind=sp_kind),dimension(nz),intent(in)  :: Ti,Tj
real(kind=dp_kind),dimension(n),intent(in) :: B
real(kind=dp_kind),dimension(n),intent(out) :: x
integer(kind=sp_kind), intent(out) :: status
real(kind=dp_kind), dimension(2), optional, intent(out) :: det

integer(kind=sp_kind),dimension(:),allocatable :: wTi,wTj
real(kind=dp_kind),dimension(:),allocatable :: A
integer(kind=sp_kind),dimension(:),allocatable :: Ap,Ai
!integer(kind=dp_kind) :: symbolic,numeric
integer(kind=sp_kind),dimension(2) :: symbolic,numeric
! As symbolic and numeric are used to store a C pointer, it is necessary to
! still use an integer(kind=dp_kind) for 64bits machines
! An other possibility : integer(kind=sp_kind),dimension(2) :: symbolic,numeric
real(kind=dp_kind),dimension(90) :: info
real(kind=dp_kind),dimension(20) :: control
integer(kind=sp_kind) :: 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 : ",trim(umfpack_return_code(status))
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 : ",trim(umfpack_return_code(int(info(1),kind=sp_kind)))
    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 : ",trim(umfpack_return_code(int(info(1),kind=sp_kind)))
    status=info(1)
endif

! free the C symbolic pointer
call umfpack_di_free_symbolic (symbolic)

! if parameter det is present, the determinant of the matrix is calculated
if (present(det) ) then
  call umfpack_di_get_determinant(det(1),det(2),numeric,info,status)  
  ! info(1) should be zero
  if (info(1) /= 0) then
      if ( (info(1) < 1) .or. (info(1) >3) ) then ! not a warning
          write(*,*) "Problem during sparse determinant : ",trim(umfpack_return_code(int(info(1),kind=sp_kind)))
      endif
      status=info(1)
  endif
endif

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 : ",trim(umfpack_return_code(int(info(1),kind=sp_kind)))
    status=info(1)
endif

! free the C numeric pointer
call umfpack_di_free_numeric (numeric)

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


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! SPARSE DETERMINANT
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine fvn_zl_sparse_det(n,nz,T,Ti,Tj,det,status)
implicit none
integer(kind=dp_kind), intent(in) :: n,nz
complex(kind=dp_kind),dimension(nz),intent(in) :: T
integer(kind=dp_kind),dimension(nz),intent(in)  :: Ti,Tj
real(kind=dp_kind), dimension(3), intent(out) :: det
integer(kind=dp_kind), intent(out) :: status

integer(kind=dp_kind),dimension(:),allocatable :: wTi,wTj
real(kind=dp_kind),dimension(:),allocatable :: Tx,Tz
real(kind=dp_kind),dimension(:),allocatable :: Ax,Az
integer(kind=dp_kind),dimension(:),allocatable :: Ap,Ai
integer(kind=dp_kind) :: symbolic,numeric
real(kind=dp_kind),dimension(90) :: info
real(kind=dp_kind),dimension(20) :: control
real(kind=dp_kind) :: Mx,Mz,Ex


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 ",trim(umfpack_return_code(int(status,kind=sp_kind)))
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 : ",trim(umfpack_return_code(int(info(1),kind=sp_kind)))
    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 : ",trim(umfpack_return_code(int(info(1),kind=sp_kind)))
    status=info(1)
endif

! free the C symbolic pointer
call umfpack_zl_free_symbolic (symbolic)

call umfpack_zl_get_determinant(det(1),det(2),det(3),numeric,info,status)
! info(1) should be zero
if (info(1) /= 0) then
      if ( (info(1) < 1) .or. (info(1) >3) ) then ! not a warning
          write(*,*) "Problem during sparse determinant : ",trim(umfpack_return_code(int(info(1),kind=sp_kind)))
      endif
    status=info(1)
endif

! free the C numeric pointer
call umfpack_zl_free_numeric (numeric)

deallocate(Ax,Az)
deallocate(Tx,Tz)
deallocate(wTi,wTj)
end subroutine


subroutine fvn_zi_sparse_det(n,nz,T,Ti,Tj,det,status)
implicit none
integer(kind=sp_kind), intent(in) :: n,nz
complex(kind=dp_kind),dimension(nz),intent(in) :: T
integer(kind=sp_kind),dimension(nz),intent(in)  :: Ti,Tj
integer(kind=sp_kind), intent(out) :: status
real(kind=dp_kind), dimension(3), intent(out) :: det

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

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 : ",trim(umfpack_return_code(status))
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 : ",trim(umfpack_return_code(int(info(1),kind=sp_kind)))
    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 : ",trim(umfpack_return_code(int(info(1),kind=sp_kind)))
    status=info(1)
endif

! free the C symbolic pointer
call umfpack_zi_free_symbolic (symbolic)

call umfpack_zi_get_determinant(det(1),det(2),det(3),numeric,info,status)
! info(1) should be zero
if (info(1) /= 0) then
      if ( (info(1) < 1) .or. (info(1) >3) ) then ! not a warning
          write(*,*) "Problem during sparse determinant : ",trim(umfpack_return_code(int(info(1),kind=sp_kind)))
      endif
    status=info(1)
endif

! free the C numeric pointer
call umfpack_zi_free_numeric (numeric)

deallocate(Ax,Az)
deallocate(Tx,Tz)
deallocate(wTi,wTj)
end subroutine

subroutine fvn_dl_sparse_det(n,nz,T,Ti,Tj,det,status)
implicit none
integer(kind=dp_kind), intent(in) :: n,nz
real(kind=dp_kind),dimension(nz),intent(in) :: T
integer(kind=dp_kind),dimension(nz),intent(in)  :: Ti,Tj
integer(kind=dp_kind), intent(out) :: status
real(kind=dp_kind), dimension(2), intent(out) :: det

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

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 : ",trim(umfpack_return_code(int(status,kind=sp_kind)))
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 : ",trim(umfpack_return_code(int(info(1),kind=sp_kind)))
    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 : ",trim(umfpack_return_code(int(info(1),kind=sp_kind)))
    status=info(1)
endif

! free the C symbolic pointer
call umfpack_dl_free_symbolic (symbolic)

call umfpack_dl_get_determinant(det(1),det(2),numeric,info,status)
! info(1) should be zero
if (info(1) /= 0) then
      if ( (info(1) < 1) .or. (info(1) >3) ) then ! not a warning
          write(*,*) "Problem during sparse determinant : ",trim(umfpack_return_code(int(info(1),kind=sp_kind)))
      endif
    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_det(n,nz,T,Ti,Tj,det,status)
implicit none
integer(kind=sp_kind), intent(in) :: n,nz
real(kind=dp_kind),dimension(nz),intent(in) :: T
integer(kind=sp_kind),dimension(nz),intent(in)  :: Ti,Tj
integer(kind=sp_kind), intent(out) :: status
real(kind=dp_kind), dimension(2), intent(out) :: det

integer(kind=sp_kind),dimension(:),allocatable :: wTi,wTj
real(kind=dp_kind),dimension(:),allocatable :: A
integer(kind=sp_kind),dimension(:),allocatable :: Ap,Ai
!integer(kind=dp_kind) :: symbolic,numeric
integer(kind=sp_kind),dimension(2) :: symbolic,numeric
! As symbolic and numeric are used to store a C pointer, it is necessary to
! still use an integer(kind=dp_kind) for 64bits machines
! An other possibility : integer(kind=sp_kind),dimension(2) :: symbolic,numeric
real(kind=dp_kind),dimension(90) :: info
real(kind=dp_kind),dimension(20) :: control
integer(kind=sp_kind) :: 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 : ",trim(umfpack_return_code(status))
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 : ",trim(umfpack_return_code(int(info(1),kind=sp_kind)))
    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 : ",trim(umfpack_return_code(int(info(1),kind=sp_kind)))
    status=info(1)
endif


! free the C symbolic pointer
call umfpack_di_free_symbolic (symbolic)

call umfpack_di_get_determinant(det(1),det(2),numeric,info,status)
! info(1) should be zero
if (info(1) /= 0) then
      if ( (info(1) < 1) .or. (info(1) >3) ) then ! not a warning
          write(*,*) "Problem during sparse determinant : ",trim(umfpack_return_code(int(info(1),kind=sp_kind)))
      endif
    status=info(1)
endif

! free the C numeric pointer
call umfpack_di_free_numeric (numeric)

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





end module fvn_sparse