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

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

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

! Define defaults control values
call umfpack_zl_defaults(control)

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

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

! free the C symbolic pointer
call umfpack_zl_free_symbolic (symbolic)

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

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


! free the C numeric pointer
call umfpack_zl_free_numeric (numeric)

x=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)
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

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

! Define defaults control values
call umfpack_zi_defaults(control)

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

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

! free the C symbolic pointer
call umfpack_zi_free_symbolic (symbolic)

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

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

! free the C numeric pointer
call umfpack_zi_free_numeric (numeric)

x=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)
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

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

! Define defaults control values
call umfpack_dl_defaults(control)

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

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

! free the C symbolic pointer
call umfpack_dl_free_symbolic (symbolic)

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

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

! free the C numeric pointer
call umfpack_dl_free_numeric (numeric)

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






subroutine fvn_di_sparse_solve(n,nz,T,Ti,Tj,B,x,status)
implicit none
integer(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

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

! Define defaults control values
call umfpack_di_defaults(control)

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

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

! free the C symbolic pointer
call umfpack_di_free_symbolic (symbolic)

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

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

! free the C numeric pointer
call umfpack_di_free_numeric (numeric)

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


end module fvn_sparse