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