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