program test_sparse
use fvn
implicit none
integer(kind=sp_kind), parameter :: nz=12
integer(kind=sp_kind), parameter :: n=5
complex(kind=dp_kind),dimension(nz) :: A
real(kind=dp_kind),dimension(nz) :: Ax,Az
complex(kind=dp_kind),dimension(n,n) :: As
integer(kind=sp_kind),dimension(nz) :: Ti,Tj
complex(kind=dp_kind),dimension(n) :: B,x
real(kind=dp_kind),dimension(n) :: Bx,Bz
integer(kind=sp_kind) :: status,i
real(kind=dp_kind),dimension(3) :: det
character(len=80) :: fmcmplx

fmcmplx='(5("(",f8.5,",",f8.5,")  "))'

! Description of the matrix in triplet form
A = (/ (2.,-1.),(3.,2.),(3.,1.),(-1.,5.),(4.,-7.),(4.,0.),(-3.,-4.),(1.,3.),(2.,0.),(2.,-2.),(6.,4.),(1.,0.) /)
B = (/ (8.,3.), (45.,1.), (-3.,-2.), (3.,0.), (19.,2.) /)
Ti = (/ 1,2,1,3,5,2,3,4,5,3,2,5 /)
Tj = (/ 1,1,2,2,2,3,3,3,3,4,5,5 /)

Ax=real(A)
Az=aimag(A)
Bx=real(B)
Bz=aimag(B)

! Reconstruction of the matrix in standard form
As=0.
do i=1,nz
    As(Ti(i),Tj(i))=A(i)
end do

! sparse routines must be fed up with 0-based indices
Ti=Ti-1
Tj=Tj-1

write(*,*) "Matrix in standard representation :"
do i=1,5
    write(*,fmcmplx) As(i,:)
end do
write(*,*)
write(*,*) "Standard determinant : ",fvn_det(5,As)
write(*,*)
write(*,*) "Right hand side :"
write(*,fmcmplx) B

! can use either specific interface, fvn_zi_sparse_det
! either generic one fvn_sparse_det
call fvn_zi_sparse_det(n,nz,Ax,Az,Ti,Tj,det,status)
write(*,*)
write(*,*) "Sparse Det = ",cmplx(det(1),det(2),kind=dp_kind)*10**det(3)
! can use either specific interface fvn_zi_sparse_solve
! either generic one fvn_sparse_solve
! parameter det is optional
call fvn_zi_sparse_solve(n,nz,Ax,Az,Ti,Tj,Bx,Bz,x,status,det)
write(*,*)
write(*,*) "Sparse Det as solve option= ",cmplx(det(1),det(2),kind=dp_kind)*10**det(3)
write(*,*)
write(*,*) "Solution :"
write(*,fmcmplx) x
write(*,*)
write(*,*) "Product matrix Solution :"
write(*,fmcmplx) matmul(As,x)
end program