program test_sparse
  ! test sparse routine zl, complex(8) and integer(8)
  use fvn
  implicit none
  integer(kind=dp_kind), parameter :: nz=12
  integer(kind=dp_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=dp_kind),dimension(nz) :: Ti,Tj
  complex(kind=dp_kind),dimension(n) :: B,x

  real(kind=dp_kind),dimension(n) :: Bx,Bz

  integer(kind=dp_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_zl_sparse_det
  ! either generic one fvn_sparse_det

  call fvn_zl_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_zl_sparse_solve
  ! either generic one fvn_sparse_solve
  ! parameter det is optional

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