program test_sparse
  ! test sparse routine dl, real(8) and integer(8)
  use fvn
  implicit none
  integer(kind=dp_kind), parameter :: nz=12
  integer(kind=dp_kind), parameter :: n=5
  real(kind=dp_kind),dimension(nz) :: A
  real(kind=dp_kind),dimension(n,n) :: As
  integer(kind=dp_kind),dimension(nz) :: Ti,Tj
  real(kind=dp_kind),dimension(n) :: B,x
  integer(kind=dp_kind) :: status,i

  real(kind=dp_kind),dimension(2) :: det

  ! Description of the matrix in triplet form
  A = (/ 2.,3.,3.,-1.,4.,4.,-3.,1.,2.,2.,6.,1. /)
  B = (/ 8., 45., -3., 3., 19./)
  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 /)
  
  ! Reconstruction of the matrix in standard form
  ! just needed for printing the matrix here
  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(*,'(5f8.4)') As(i,:)
  end do
  write(*,*)
  write(*,*) "Standard determinant =", fvn_det(5,As)
  write(*,*)
  write(*,'("Right hand side :",5f8.4)') B
  
  ! can use either specific interface, fvn_dl_sparse_det
  ! either generic one fvn_sparse_det
  call fvn_dl_sparse_det(n,nz,A,Ti,Tj,det,status)
  write(*,*)

  write(*,*) "Sparse Det = ",det(1)*10**det(2)

  ! can use either specific interface fvn_dl_sparse_solve
  ! either generic one fvn_sparse_solve
  ! parameter det is optional
  call fvn_dl_sparse_solve(n,nz,A,Ti,Tj,B,x,status,det)
  write(*,*)

  write(*,*) "Sparse Det as solve option = ",det(1)*10**det(2)

  write(*,*)
  write(*,'("Solution :",5f8.4)') x
  write(*,*)
  write(*,'("Product matrix Solution :",5f8.4)') matmul(As,x)
  end program