1) added fvn_sparse_det routine to compute determinant of a sparse matrix

2) added optional parameter det to fvn_sparse_solve to compute determinant during solving 3) split the test_sparse.f90 into 4 files to test each specific interface git-svn-id: https://lxsd.femto-st.fr/svn/fvn@64 b657c933-2333-4658-acf2-d3c7c2708721

1) added fvn_sparse_det routine to compute determinant of a sparse matrix
2) added optional parameter det to fvn_sparse_solve to compute determinant during solving 3) split the test_sparse.f90 into 4 files to test each specific interface git-svn-id: https://lxsd.femto-st.fr/svn/fvn@64 b657c933-2333-4658-acf2-d3c7c2708721
wdaniau
1 parent 8ba5c9c788
Showing 8 changed files with 612 additions and 42 deletions Side-by-side Diff
fvn_sparse/fvn_sparse.f90
fvn_sparse/umfpack_wrapper.c
fvn_test/Makefile
fvn_test/test_sparse.f90
fvn_test/test_sparse_di.f90
fvn_test/test_sparse_dl.f90
fvn_test/test_sparse_zi.f90
fvn_test/test_sparse_zl.f90
@@ -7,6 +7,9 @@
     module procedure fvn_zl_sparse_solve,fvn_zi_sparse_solve,fvn_dl_sparse_solve,fvn_di_sparse_solve
 end interface fvn_sparse_solve
  
+interface fvn_sparse_det
+    module procedure fvn_zl_sparse_det,fvn_zi_sparse_det,fvn_dl_sparse_det,fvn_di_sparse_det
+end interface fvn_sparse_det
 contains
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 !
@@ -30,7 +33,7 @@
 ! * = 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)
+subroutine fvn_zl_sparse_solve(n,nz,T,Ti,Tj,B,x,status,det)
 implicit none
 integer(kind=dp_kind), intent(in) :: n,nz
 complex(kind=dp_kind),dimension(nz),intent(in) :: T
@@ -38,6 +41,7 @@
 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
+complex(kind=dp_kind), optional, intent(out) :: det
  
 integer(kind=dp_kind),dimension(:),allocatable :: wTi,wTj
 real(kind=dp_kind),dimension(:),allocatable :: Tx,Tz
@@ -48,6 +52,7 @@
 real(kind=dp_kind),dimension(90) :: info
 real(kind=dp_kind),dimension(20) :: control
 integer(kind=dp_kind) :: sys
+real(kind=dp_kind) :: Mx,Mz,Ex
  
  
 status=0
@@ -92,6 +97,19 @@
 ! free the C symbolic pointer
 call umfpack_zl_free_symbolic (symbolic)
  
+! if parameter det is present, the determinant of the matrix is calculated
+if (present(det) ) then
+  call umfpack_zl_get_determinant(Mx,Mz,Ex,numeric,info,status)
+  det=cmplx(Mx,Mz,kind=dp_kind)*10**Ex
+  ! info(1) should be zero
+  if (info(1) /= 0) then
+      write(*,*) "Problem during sparse determinant, returned code : ",info(1)
+      status=info(1)
+  endif
+endif
+
+
+
 allocate(bx(n),bz(n),xx(n),xz(n))
 bx=dble(B)
 bz=aimag(B)
@@ -122,7 +140,7 @@
  
  
  
-subroutine fvn_zi_sparse_solve(n,nz,T,Ti,Tj,B,x,status)
+subroutine fvn_zi_sparse_solve(n,nz,T,Ti,Tj,B,x,status,det)
 implicit none
 integer(kind=sp_kind), intent(in) :: n,nz
 complex(kind=dp_kind),dimension(nz),intent(in) :: T
@@ -130,6 +148,7 @@
 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
+complex(kind=dp_kind), optional, intent(out) :: det
  
 integer(kind=sp_kind),dimension(:),allocatable :: wTi,wTj
 real(kind=dp_kind),dimension(:),allocatable :: Tx,Tz
@@ -144,6 +163,7 @@
 real(kind=dp_kind),dimension(90) :: info
 real(kind=dp_kind),dimension(20) :: control
 integer(kind=sp_kind) :: sys
+real(kind=dp_kind) :: Mx,Mz,Ex
  
 status=0
 ! we use a working copy of Ti and Tj to perform 1-based to 0-based translation
@@ -186,6 +206,20 @@
 ! free the C symbolic pointer
 call umfpack_zi_free_symbolic (symbolic)
  
+! if parameter det is present, the determinant of the matrix is calculated
+if (present(det) ) then
+  call umfpack_zi_get_determinant(Mx,Mz,Ex,numeric,info,status)
+  det=cmplx(Mx,Mz,kind=dp_kind)*10**Ex
+  ! info(1) should be zero
+  if (info(1) /= 0) then
+      write(*,*) "Problem during sparse determinant, returned code : ",info(1)
+      status=info(1)
+  endif
+endif
+
+
+
+
 allocate(bx(n),bz(n),xx(n),xz(n))
 bx=dble(B)
 bz=aimag(B)
@@ -216,7 +250,7 @@
  
  
  
-subroutine fvn_dl_sparse_solve(n,nz,T,Ti,Tj,B,x,status)
+subroutine fvn_dl_sparse_solve(n,nz,T,Ti,Tj,B,x,status,det)
 implicit none
 integer(kind=dp_kind), intent(in) :: n,nz
 real(kind=dp_kind),dimension(nz),intent(in) :: T
@@ -224,6 +258,7 @@
 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
+real(kind=dp_kind), optional, intent(out) :: det
  
 integer(kind=dp_kind),dimension(:),allocatable :: wTi,wTj
 real(kind=dp_kind),dimension(:),allocatable :: A
@@ -233,6 +268,7 @@
 real(kind=dp_kind),dimension(90) :: info
 real(kind=dp_kind),dimension(20) :: control
 integer(kind=dp_kind) :: sys
+real(kind=dp_kind) :: Mx,Ex
  
 status=0
 ! we use a working copy of Ti and Tj to perform 1-based to 0-based translation
@@ -271,6 +307,18 @@
 ! free the C symbolic pointer
 call umfpack_dl_free_symbolic (symbolic)
  
+! if parameter det is present, the determinant of the matrix is calculated
+if (present(det) ) then
+  call umfpack_dl_get_determinant(Mx,Ex,numeric,info,status)
+  det=Mx*10**Ex
+  ! info(1) should be zero
+  if (info(1) /= 0) then
+      write(*,*) "Problem during sparse determinant, returned code : ",info(1)
+      status=info(1)
+  endif
+endif
+
+
 sys=0
 ! sys may be used to define type of solving -> see umfpack.h
  
@@ -294,7 +342,7 @@
  
  
  
-subroutine fvn_di_sparse_solve(n,nz,T,Ti,Tj,B,x,status)
+subroutine fvn_di_sparse_solve(n,nz,T,Ti,Tj,B,x,status,det)
 implicit none
 integer(kind=sp_kind), intent(in) :: n,nz
 real(kind=dp_kind),dimension(nz),intent(in) :: T
@@ -302,6 +350,7 @@
 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
+real(kind=dp_kind), optional, intent(out) :: det
  
 integer(kind=sp_kind),dimension(:),allocatable :: wTi,wTj
 real(kind=dp_kind),dimension(:),allocatable :: A
@@ -314,6 +363,7 @@
 real(kind=dp_kind),dimension(90) :: info
 real(kind=dp_kind),dimension(20) :: control
 integer(kind=sp_kind) :: sys
+real(kind=dp_kind) :: Mx,Ex
  
 status=0
 ! we use a working copy of Ti and Tj to perform 1-based to 0-based translation
  
@@ -352,9 +402,19 @@
 ! free the C symbolic pointer
 call umfpack_di_free_symbolic (symbolic)
  
+! if parameter det is present, the determinant of the matrix is calculated
+if (present(det) ) then
+  call umfpack_di_get_determinant(Mx,Ex,numeric,info,status)
+  det=Mx*10**Ex
+  ! info(1) should be zero
+  if (info(1) /= 0) then
+      write(*,*) "Problem during sparse determinant, returned code : ",info(1)
+      status=info(1)
+  endif
+endif
+
 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
@@ -369,6 +429,315 @@
 deallocate(A)
 deallocate(wTi,wTj)
 end subroutine
+
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!
+! SPARSE DETERMINANT
+!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+subroutine fvn_zl_sparse_det(n,nz,T,Ti,Tj,det,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),intent(out) :: det
+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(90) :: info
+real(kind=dp_kind),dimension(20) :: control
+real(kind=dp_kind) :: Mx,Mz,Ex
+
+
+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)
+
+call umfpack_zl_get_determinant(Mx,Mz,Ex,numeric,info,status)
+! info(1) should be zero
+if (info(1) /= 0) then
+    write(*,*) "Problem during sparse determinant, returned code : ",info(1)
+    status=info(1)
+endif
+det=cmplx(Mx,Mz,kind=dp_kind)*10**Ex
+
+! free the C numeric pointer
+call umfpack_zl_free_numeric (numeric)
+
+deallocate(Ax,Az)
+deallocate(Tx,Tz)
+deallocate(wTi,wTj)
+end subroutine
+
+
+subroutine fvn_zi_sparse_det(n,nz,T,Ti,Tj,det,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
+integer(kind=sp_kind), intent(out) :: status
+complex(kind=dp_kind), intent(out) :: det
+
+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(90) :: info
+real(kind=dp_kind),dimension(20) :: control
+real(kind=dp_kind) :: Mx,Mz,Ex
+
+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)
+
+call umfpack_zi_get_determinant(Mx,Mz,Ex,numeric,info,status)
+! info(1) should be zero
+if (info(1) /= 0) then
+    write(*,*) "Problem during sparse determinant, returned code : ",info(1)
+    status=info(1)
+endif
+det=cmplx(Mx,Mz,kind=dp_kind)*10**Ex
+
+! free the C numeric pointer
+call umfpack_zi_free_numeric (numeric)
+
+deallocate(Ax,Az)
+deallocate(Tx,Tz)
+deallocate(wTi,wTj)
+end subroutine
+
+subroutine fvn_dl_sparse_det(n,nz,T,Ti,Tj,det,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
+integer(kind=dp_kind), intent(out) :: status
+real(kind=dp_kind), intent(out) :: det
+
+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
+real(kind=dp_kind) :: Mx,Ex
+
+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)
+
+call umfpack_dl_get_determinant(Mx,Ex,numeric,info,status)
+! info(1) should be zero
+if (info(1) /= 0) then
+    write(*,*) "Problem during sparse determinant, returned code : ",info(1)
+    status=info(1)
+endif
+det=Mx*10**Ex
+
+! free the C numeric pointer
+call umfpack_dl_free_numeric (numeric)
+
+deallocate(A)
+deallocate(wTi,wTj)
+end subroutine
+
+
+subroutine fvn_di_sparse_det(n,nz,T,Ti,Tj,det,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
+integer(kind=sp_kind), intent(out) :: status
+real(kind=dp_kind), intent(out) :: det
+
+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
+real(kind=dp_kind) :: Mx,Ex
+
+
+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)
+
+call umfpack_di_get_determinant(Mx,Ex,numeric,info,status)
+det=Mx*10**Ex
+! info(1) should be zero
+if (info(1) /= 0) then
+    write(*,*) "Problem during sparse determinant, returned code : ",info(1)
+    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
@@ -72,6 +72,11 @@
 }
  
  
+/* get determinant */
+void umfpack_zl_get_determinant_ ( double *Mx, double *Mz, double *Ex, void **Numeric, double Info [UMFPACK_INFO], UF_long *status )
+{
+    *status = umfpack_zl_get_determinant (Mx,Mz,Ex,*Numeric, Info);
+}
  
  
 /* complex(8) and integer(4) */
@@ -131,6 +136,14 @@
     *status = umfpack_zi_triplet_to_col (*m, *n, *nz, Ti, Tj, Tx, Tz, Ap, Ai, Ax, Az, (int *) NULL);
 }
  
+
+/* get determinant */
+void umfpack_zi_get_determinant_ ( double *Mx, double *Mz, double *Ex, void **Numeric, double Info [UMFPACK_INFO], int *status )
+{
+    *status = umfpack_zi_get_determinant (Mx,Mz,Ex,*Numeric, Info);
+}
+
+
 /* real(8) and integer(8) */
  
 /* defaults  */
@@ -178,6 +191,14 @@
     *status = umfpack_dl_triplet_to_col (*m, *n, *nz, Ti, Tj, T, Ap, Ai, A, (UF_long *) NULL);
 }
  
+/* get determinant */
+void umfpack_dl_get_determinant_ ( double *Mx, double *Ex, void **Numeric, double Info [UMFPACK_INFO], UF_long *status )
+{
+    *status = umfpack_dl_get_determinant (Mx,Ex,*Numeric, Info);
+}
+
+
+
 /* real(8) and integer(4) */
  
 /* defaults  */
@@ -223,5 +244,11 @@
     double A [ ], int *status)
 {
     *status = umfpack_di_triplet_to_col (*m, *n, *nz, Ti, Tj, T, Ap, Ai, A, (int *) NULL);
+}
+
+/* get determinant */
+void umfpack_di_get_determinant_ ( double *Mx, double *Ex, void **Numeric, double Info [UMFPACK_INFO], int *status )
+{
+    *status = umfpack_di_get_determinant (Mx,Ex,*Numeric, Info);
 }
@@ -2,10 +2,11 @@
 include $(BTREE)/Make.inc
  
 programs = test_fac$(exext) test_matinv$(exext) test_specfunc$(exext)  \
-test_det$(exext) test_matcon$(exext) test_matev$(exext) test_sparse$(exext) test_inter1d$(exext) \
+test_det$(exext) test_matcon$(exext) test_matev$(exext) test_inter1d$(exext) \
 test_inter2d$(exext) test_inter3d$(exext) test_akima$(exext) test_lsp$(exext) test_muller$(exext) \
 test_integ$(exext) test_bsyn$(exext) test_bsjn$(exext) test_bskn$(exext) test_bsin$(exext) test_operators$(exext) test_ze1$(exext) \
-test_besri$(exext) test_besrj$(exext) test_bestime$(exext)
+test_besri$(exext) test_besrj$(exext) test_bestime$(exext) \
+test_sparse_zi$(exext) test_sparse_zl$(exext) test_sparse_di$(exext) test_sparse_dl$(exext)
  
 prog:$(programs)
  
-program test_sparse
-use fvn_sparse
-implicit none
-integer(kind=ip_kind), parameter :: nz=12
-integer(kind=ip_kind), parameter :: n=5
-real(kind=dp_kind),dimension(nz) :: A
-real(kind=dp_kind),dimension(n,n) :: As
-integer(kind=ip_kind),dimension(nz) :: Ti,Tj
-real(kind=dp_kind),dimension(n) :: B,x
-integer(kind=ip_kind) :: status,i
-! 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
-write(*,*) "Matrix in standard representation :"
-do i=1,5
-    write(*,'(5f8.4)') As(i,:)
-end do
-write(*,*)
-write(*,'("Right hand side :",5f8.4)') B
-
-!specific routine that will be used here
-!call fvn_di_sparse_solve(n,nz,A,Ti,Tj,B,x,status)
-call fvn_di_sparse_solve(n,nz,A,Ti,Tj,B,x,status)
-write(*,'("Solution :",5f8.4)') x
-write(*,'("Product matrix Solution :",5f8.4)') matmul(As,x)
-end program
+program test_sparse
+! test sparse routine di, real(8) and integer(4)
+use fvn
+implicit none
+integer(kind=sp_kind), parameter :: nz=12
+integer(kind=sp_kind), parameter :: n=5
+real(kind=dp_kind),dimension(nz) :: A
+real(kind=dp_kind),dimension(n,n) :: As
+integer(kind=sp_kind),dimension(nz) :: Ti,Tj
+real(kind=dp_kind),dimension(n) :: B,x
+integer(kind=sp_kind) :: status,i
+real(kind=dp_kind) :: 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
+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_di_sparse_det
+! either generic one fvn_sparse_det
+call fvn_di_sparse_det(n,nz,A,Ti,Tj,det,status)
+write(*,*)
+write(*,*) "Sparse Det = ",det
+! can use either specific interface fvn_di_sparse_solve
+! either generic one fvn_sparse_solve
+! parameter det is optional
+call fvn_di_sparse_solve(n,nz,A,Ti,Tj,B,x,status,det)
+write(*,*)
+write(*,*) "Sparse Det as solve option = ",det
+write(*,*)
+write(*,'("Solution :",5f8.4)') x
+write(*,*)
+write(*,'("Product matrix Solution :",5f8.4)') matmul(As,x)
+end program
+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) :: 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
+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
+! 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
+write(*,*)
+write(*,'("Solution :",5f8.4)') x
+write(*,*)
+write(*,'("Product matrix Solution :",5f8.4)') matmul(As,x)
+end program
+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
+complex(kind=dp_kind),dimension(n,n) :: As
+integer(kind=sp_kind),dimension(nz) :: Ti,Tj
+complex(kind=dp_kind),dimension(n) :: B,x
+integer(kind=sp_kind) :: status,i
+complex(kind=dp_kind) :: 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 /)
+
+! Reconstruction of the matrix in standard form
+As=0.
+do i=1,nz
+    As(Ti(i),Tj(i))=A(i)
+end do
+
+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,A,Ti,Tj,det,status)
+write(*,*)
+write(*,*) "Sparse Det = ",det
+! 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,A,Ti,Tj,B,x,status,det)
+write(*,*)
+write(*,*) "Sparse Det as solve option= ",det
+write(*,*)
+write(*,*) "Solution :"
+write(*,fmcmplx) x
+write(*,*)
+write(*,*) "Product matrix Solution :"
+write(*,fmcmplx) matmul(As,x)
+end program
+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
+complex(kind=dp_kind),dimension(n,n) :: As
+integer(kind=dp_kind),dimension(nz) :: Ti,Tj
+complex(kind=dp_kind),dimension(n) :: B,x
+integer(kind=dp_kind) :: status,i
+complex(kind=dp_kind) :: 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 /)
+
+! Reconstruction of the matrix in standard form
+As=0.
+do i=1,nz
+    As(Ti(i),Tj(i))=A(i)
+end do
+
+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,A,Ti,Tj,det,status)
+write(*,*)
+write(*,*) "Sparse Det = ",det
+! 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,A,Ti,Tj,B,x,status,det)
+write(*,*)
+write(*,*) "Sparse Det as solve option= ",det
+write(*,*)
+write(*,*) "Solution :"
+write(*,fmcmplx) x
+write(*,*)
+write(*,*) "Product matrix Solution :"
+write(*,fmcmplx) matmul(As,x)
+end program
...	...	@@ -7,6 +7,9 @@
7	7	module procedure fvn_zl_sparse_solve,fvn_zi_sparse_solve,fvn_dl_sparse_solve,fvn_di_sparse_solve
8	8	end interface fvn_sparse_solve
9	9
	10	+interface fvn_sparse_det
	11	+ module procedure fvn_zl_sparse_det,fvn_zi_sparse_det,fvn_dl_sparse_det,fvn_di_sparse_det
	12	+end interface fvn_sparse_det
10	13	contains
11	14	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
12	15	!
...	...	@@ -30,7 +33,7 @@
30	33	! * = zl : double complex + integer(kind=dp_kind)
31	34	! * = zi : double complex + integer(kind=sp_kind)
32	35	!
33		-subroutine fvn_zl_sparse_solve(n,nz,T,Ti,Tj,B,x,status)
	36	+subroutine fvn_zl_sparse_solve(n,nz,T,Ti,Tj,B,x,status,det)
34	37	implicit none
35	38	integer(kind=dp_kind), intent(in) :: n,nz
36	39	complex(kind=dp_kind),dimension(nz),intent(in) :: T
...	...	@@ -38,6 +41,7 @@
38	41	complex(kind=dp_kind),dimension(n),intent(in) :: B
39	42	complex(kind=dp_kind),dimension(n),intent(out) :: x
40	43	integer(kind=dp_kind), intent(out) :: status
	44	+complex(kind=dp_kind), optional, intent(out) :: det
41	45
42	46	integer(kind=dp_kind),dimension(:),allocatable :: wTi,wTj
43	47	real(kind=dp_kind),dimension(:),allocatable :: Tx,Tz
...	...	@@ -48,6 +52,7 @@
48	52	real(kind=dp_kind),dimension(90) :: info
49	53	real(kind=dp_kind),dimension(20) :: control
50	54	integer(kind=dp_kind) :: sys
	55	+real(kind=dp_kind) :: Mx,Mz,Ex
51	56
52	57
53	58	status=0
...	...	@@ -92,6 +97,19 @@
92	97	! free the C symbolic pointer
93	98	call umfpack_zl_free_symbolic (symbolic)
94	99
	100	+! if parameter det is present, the determinant of the matrix is calculated
	101	+if (present(det) ) then
	102	+ call umfpack_zl_get_determinant(Mx,Mz,Ex,numeric,info,status)
	103	+ det=cmplx(Mx,Mz,kind=dp_kind)10*Ex
	104	+ ! info(1) should be zero
	105	+ if (info(1) /= 0) then
	106	+ write(,) "Problem during sparse determinant, returned code : ",info(1)
	107	+ status=info(1)
	108	+ endif
	109	+endif
	110	+
	111	+
	112	+
95	113	allocate(bx(n),bz(n),xx(n),xz(n))
96	114	bx=dble(B)
97	115	bz=aimag(B)
...	...	@@ -122,7 +140,7 @@
122	140
123	141
124	142
125		-subroutine fvn_zi_sparse_solve(n,nz,T,Ti,Tj,B,x,status)
	143	+subroutine fvn_zi_sparse_solve(n,nz,T,Ti,Tj,B,x,status,det)
126	144	implicit none
127	145	integer(kind=sp_kind), intent(in) :: n,nz
128	146	complex(kind=dp_kind),dimension(nz),intent(in) :: T
...	...	@@ -130,6 +148,7 @@
130	148	complex(kind=dp_kind),dimension(n),intent(in) :: B
131	149	complex(kind=dp_kind),dimension(n),intent(out) :: x
132	150	integer(kind=sp_kind), intent(out) :: status
	151	+complex(kind=dp_kind), optional, intent(out) :: det
133	152
134	153	integer(kind=sp_kind),dimension(:),allocatable :: wTi,wTj
135	154	real(kind=dp_kind),dimension(:),allocatable :: Tx,Tz
...	...	@@ -144,6 +163,7 @@
144	163	real(kind=dp_kind),dimension(90) :: info
145	164	real(kind=dp_kind),dimension(20) :: control
146	165	integer(kind=sp_kind) :: sys
	166	+real(kind=dp_kind) :: Mx,Mz,Ex
147	167
148	168	status=0
149	169	! we use a working copy of Ti and Tj to perform 1-based to 0-based translation
...	...	@@ -186,6 +206,20 @@
186	206	! free the C symbolic pointer
187	207	call umfpack_zi_free_symbolic (symbolic)
188	208
	209	+! if parameter det is present, the determinant of the matrix is calculated
	210	+if (present(det) ) then
	211	+ call umfpack_zi_get_determinant(Mx,Mz,Ex,numeric,info,status)
	212	+ det=cmplx(Mx,Mz,kind=dp_kind)10*Ex
	213	+ ! info(1) should be zero
	214	+ if (info(1) /= 0) then
	215	+ write(,) "Problem during sparse determinant, returned code : ",info(1)
	216	+ status=info(1)
	217	+ endif
	218	+endif
	219	+
	220	+
	221	+
	222	+
189	223	allocate(bx(n),bz(n),xx(n),xz(n))
190	224	bx=dble(B)
191	225	bz=aimag(B)
...	...	@@ -216,7 +250,7 @@
216	250
217	251
218	252
219		-subroutine fvn_dl_sparse_solve(n,nz,T,Ti,Tj,B,x,status)
	253	+subroutine fvn_dl_sparse_solve(n,nz,T,Ti,Tj,B,x,status,det)
220	254	implicit none
221	255	integer(kind=dp_kind), intent(in) :: n,nz
222	256	real(kind=dp_kind),dimension(nz),intent(in) :: T
...	...	@@ -224,6 +258,7 @@
224	258	real(kind=dp_kind),dimension(n),intent(in) :: B
225	259	real(kind=dp_kind),dimension(n),intent(out) :: x
226	260	integer(kind=dp_kind), intent(out) :: status
	261	+real(kind=dp_kind), optional, intent(out) :: det
227	262
228	263	integer(kind=dp_kind),dimension(:),allocatable :: wTi,wTj
229	264	real(kind=dp_kind),dimension(:),allocatable :: A
...	...	@@ -233,6 +268,7 @@
233	268	real(kind=dp_kind),dimension(90) :: info
234	269	real(kind=dp_kind),dimension(20) :: control
235	270	integer(kind=dp_kind) :: sys
	271	+real(kind=dp_kind) :: Mx,Ex
236	272
237	273	status=0
238	274	! we use a working copy of Ti and Tj to perform 1-based to 0-based translation
...	...	@@ -271,6 +307,18 @@
271	307	! free the C symbolic pointer
272	308	call umfpack_dl_free_symbolic (symbolic)
273	309
	310	+! if parameter det is present, the determinant of the matrix is calculated
	311	+if (present(det) ) then
	312	+ call umfpack_dl_get_determinant(Mx,Ex,numeric,info,status)
	313	+ det=Mx10*Ex
	314	+ ! info(1) should be zero
	315	+ if (info(1) /= 0) then
	316	+ write(,) "Problem during sparse determinant, returned code : ",info(1)
	317	+ status=info(1)
	318	+ endif
	319	+endif
	320	+
	321	+
274	322	sys=0
275	323	! sys may be used to define type of solving -> see umfpack.h
276	324
...	...	@@ -294,7 +342,7 @@
294	342
295	343
296	344
297		-subroutine fvn_di_sparse_solve(n,nz,T,Ti,Tj,B,x,status)
	345	+subroutine fvn_di_sparse_solve(n,nz,T,Ti,Tj,B,x,status,det)
298	346	implicit none
299	347	integer(kind=sp_kind), intent(in) :: n,nz
300	348	real(kind=dp_kind),dimension(nz),intent(in) :: T
...	...	@@ -302,6 +350,7 @@
302	350	real(kind=dp_kind),dimension(n),intent(in) :: B
303	351	real(kind=dp_kind),dimension(n),intent(out) :: x
304	352	integer(kind=sp_kind), intent(out) :: status
	353	+real(kind=dp_kind), optional, intent(out) :: det
305	354
306	355	integer(kind=sp_kind),dimension(:),allocatable :: wTi,wTj
307	356	real(kind=dp_kind),dimension(:),allocatable :: A
...	...	@@ -314,6 +363,7 @@
314	363	real(kind=dp_kind),dimension(90) :: info
315	364	real(kind=dp_kind),dimension(20) :: control
316	365	integer(kind=sp_kind) :: sys
	366	+real(kind=dp_kind) :: Mx,Ex
317	367
318	368	status=0
319	369	! we use a working copy of Ti and Tj to perform 1-based to 0-based translation
320	370
...	...	@@ -352,9 +402,19 @@
352	402	! free the C symbolic pointer
353	403	call umfpack_di_free_symbolic (symbolic)
354	404
	405	+! if parameter det is present, the determinant of the matrix is calculated
	406	+if (present(det) ) then
	407	+ call umfpack_di_get_determinant(Mx,Ex,numeric,info,status)
	408	+ det=Mx10*Ex
	409	+ ! info(1) should be zero
	410	+ if (info(1) /= 0) then
	411	+ write(,) "Problem during sparse determinant, returned code : ",info(1)
	412	+ status=info(1)
	413	+ endif
	414	+endif
	415	+
355	416	sys=0
356	417	! sys may be used to define type of solving -> see umfpack.h
357		-
358	418	! Solving
359	419	call umfpack_di_solve (sys, Ap, Ai, A, x, B, numeric, control, info)
360	420	! info(1) should be zero
...	...	@@ -369,6 +429,315 @@
369	429	deallocate(A)
370	430	deallocate(wTi,wTj)
371	431	end subroutine
	432	+
	433	+
	434	+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	435	+!
	436	+! SPARSE DETERMINANT
	437	+!
	438	+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	439	+subroutine fvn_zl_sparse_det(n,nz,T,Ti,Tj,det,status)
	440	+implicit none
	441	+integer(kind=dp_kind), intent(in) :: n,nz
	442	+complex(kind=dp_kind),dimension(nz),intent(in) :: T
	443	+integer(kind=dp_kind),dimension(nz),intent(in) :: Ti,Tj
	444	+complex(kind=dp_kind),intent(out) :: det
	445	+integer(kind=dp_kind), intent(out) :: status
	446	+
	447	+integer(kind=dp_kind),dimension(:),allocatable :: wTi,wTj
	448	+real(kind=dp_kind),dimension(:),allocatable :: Tx,Tz
	449	+real(kind=dp_kind),dimension(:),allocatable :: Ax,Az
	450	+integer(kind=dp_kind),dimension(:),allocatable :: Ap,Ai
	451	+integer(kind=dp_kind) :: symbolic,numeric
	452	+real(kind=dp_kind),dimension(90) :: info
	453	+real(kind=dp_kind),dimension(20) :: control
	454	+real(kind=dp_kind) :: Mx,Mz,Ex
	455	+
	456	+
	457	+status=0
	458	+
	459	+! we use a working copy of Ti and Tj to perform 1-based to 0-based translation
	460	+! Tx and Tz are the real and imaginary parts of T
	461	+allocate(wTi(nz),wTj(nz))
	462	+allocate(Tx(nz),Tz(nz))
	463	+Tx=dble(T)
	464	+Tz=aimag(T)
	465	+wTi=Ti-1
	466	+wTj=Tj-1
	467	+allocate(Ax(nz),Az(nz))
	468	+allocate(Ap(n+1),Ai(nz))
	469	+
	470	+! perform the triplet to compressed column form -> Ap,Ai,Ax,Az
	471	+call umfpack_zl_triplet_to_col(n,n,nz,wTi,wTj,Tx,Tz,Ap,Ai,Ax,Az,status)
	472	+! if status is not zero a problem has occured
	473	+if (status /= 0) then
	474	+ write(,) "Problem during umfpack_zl_triplet_to_col"
	475	+endif
	476	+
	477	+! Define defaults control values
	478	+call umfpack_zl_defaults(control)
	479	+
	480	+! Symbolic analysis
	481	+call umfpack_zl_symbolic(n,n,Ap,Ai,Ax,Az,symbolic, control, info)
	482	+! info(1) should be zero
	483	+if (info(1) /= 0) then
	484	+ write(,) "Problem during symbolic analysis"
	485	+ status=info(1)
	486	+endif
	487	+
	488	+! Numerical factorization
	489	+call umfpack_zl_numeric (Ap, Ai, Ax, Az, symbolic, numeric, control, info)
	490	+! info(1) should be zero
	491	+if (info(1) /= 0) then
	492	+ write(,) "Problem during numerical factorization"
	493	+ status=info(1)
	494	+endif
	495	+
	496	+! free the C symbolic pointer
	497	+call umfpack_zl_free_symbolic (symbolic)
	498	+
	499	+call umfpack_zl_get_determinant(Mx,Mz,Ex,numeric,info,status)
	500	+! info(1) should be zero
	501	+if (info(1) /= 0) then
	502	+ write(,) "Problem during sparse determinant, returned code : ",info(1)
	503	+ status=info(1)
	504	+endif
	505	+det=cmplx(Mx,Mz,kind=dp_kind)10*Ex
	506	+
	507	+! free the C numeric pointer
	508	+call umfpack_zl_free_numeric (numeric)
	509	+
	510	+deallocate(Ax,Az)
	511	+deallocate(Tx,Tz)
	512	+deallocate(wTi,wTj)
	513	+end subroutine
	514	+
	515	+
	516	+subroutine fvn_zi_sparse_det(n,nz,T,Ti,Tj,det,status)
	517	+implicit none
	518	+integer(kind=sp_kind), intent(in) :: n,nz
	519	+complex(kind=dp_kind),dimension(nz),intent(in) :: T
	520	+integer(kind=sp_kind),dimension(nz),intent(in) :: Ti,Tj
	521	+integer(kind=sp_kind), intent(out) :: status
	522	+complex(kind=dp_kind), intent(out) :: det
	523	+
	524	+integer(kind=sp_kind),dimension(:),allocatable :: wTi,wTj
	525	+real(kind=dp_kind),dimension(:),allocatable :: Tx,Tz
	526	+real(kind=dp_kind),dimension(:),allocatable :: Ax,Az
	527	+integer(kind=sp_kind),dimension(:),allocatable :: Ap,Ai
	528	+!integer(kind=dp_kind) :: symbolic,numeric
	529	+integer(kind=sp_kind),dimension(2) :: symbolic,numeric
	530	+! As symbolic and numeric are used to store a C pointer, it is necessary to
	531	+! still use an integer(kind=dp_kind) for 64bits machines
	532	+! An other possibility : integer(kind=sp_kind),dimension(2) :: symbolic,numeric
	533	+real(kind=dp_kind),dimension(90) :: info
	534	+real(kind=dp_kind),dimension(20) :: control
	535	+real(kind=dp_kind) :: Mx,Mz,Ex
	536	+
	537	+status=0
	538	+! we use a working copy of Ti and Tj to perform 1-based to 0-based translation
	539	+! Tx and Tz are the real and imaginary parts of T
	540	+allocate(wTi(nz),wTj(nz))
	541	+allocate(Tx(nz),Tz(nz))
	542	+Tx=dble(T)
	543	+Tz=aimag(T)
	544	+wTi=Ti-1
	545	+wTj=Tj-1
	546	+allocate(Ax(nz),Az(nz))
	547	+allocate(Ap(n+1),Ai(nz))
	548	+
	549	+! perform the triplet to compressed column form -> Ap,Ai,Ax,Az
	550	+call umfpack_zi_triplet_to_col(n,n,nz,wTi,wTj,Tx,Tz,Ap,Ai,Ax,Az,status)
	551	+! if status is not zero a problem has occured
	552	+if (status /= 0) then
	553	+ write(,) "Problem during umfpack_zl_triplet_to_col"
	554	+endif
	555	+
	556	+! Define defaults control values
	557	+call umfpack_zi_defaults(control)
	558	+
	559	+! Symbolic analysis
	560	+call umfpack_zi_symbolic(n,n,Ap,Ai,Ax,Az,symbolic, control, info)
	561	+! info(1) should be zero
	562	+if (info(1) /= 0) then
	563	+ write(,) "Problem during symbolic analysis"
	564	+ status=info(1)
	565	+endif
	566	+
	567	+! Numerical factorization
	568	+call umfpack_zi_numeric (Ap, Ai, Ax, Az, symbolic, numeric, control, info)
	569	+! info(1) should be zero
	570	+if (info(1) /= 0) then
	571	+ write(,) "Problem during numerical factorization"
	572	+ status=info(1)
	573	+endif
	574	+
	575	+! free the C symbolic pointer
	576	+call umfpack_zi_free_symbolic (symbolic)
	577	+
	578	+call umfpack_zi_get_determinant(Mx,Mz,Ex,numeric,info,status)
	579	+! info(1) should be zero
	580	+if (info(1) /= 0) then
	581	+ write(,) "Problem during sparse determinant, returned code : ",info(1)
	582	+ status=info(1)
	583	+endif
	584	+det=cmplx(Mx,Mz,kind=dp_kind)10*Ex
	585	+
	586	+! free the C numeric pointer
	587	+call umfpack_zi_free_numeric (numeric)
	588	+
	589	+deallocate(Ax,Az)
	590	+deallocate(Tx,Tz)
	591	+deallocate(wTi,wTj)
	592	+end subroutine
	593	+
	594	+subroutine fvn_dl_sparse_det(n,nz,T,Ti,Tj,det,status)
	595	+implicit none
	596	+integer(kind=dp_kind), intent(in) :: n,nz
	597	+real(kind=dp_kind),dimension(nz),intent(in) :: T
	598	+integer(kind=dp_kind),dimension(nz),intent(in) :: Ti,Tj
	599	+integer(kind=dp_kind), intent(out) :: status
	600	+real(kind=dp_kind), intent(out) :: det
	601	+
	602	+integer(kind=dp_kind),dimension(:),allocatable :: wTi,wTj
	603	+real(kind=dp_kind),dimension(:),allocatable :: A
	604	+integer(kind=dp_kind),dimension(:),allocatable :: Ap,Ai
	605	+!integer(kind=dp_kind) :: symbolic,numeric
	606	+integer(kind=dp_kind) :: symbolic,numeric
	607	+real(kind=dp_kind),dimension(90) :: info
	608	+real(kind=dp_kind),dimension(20) :: control
	609	+real(kind=dp_kind) :: Mx,Ex
	610	+
	611	+status=0
	612	+! we use a working copy of Ti and Tj to perform 1-based to 0-based translation
	613	+allocate(wTi(nz),wTj(nz))
	614	+wTi=Ti-1
	615	+wTj=Tj-1
	616	+allocate(A(nz))
	617	+allocate(Ap(n+1),Ai(nz))
	618	+
	619	+! perform the triplet to compressed column form -> Ap,Ai,Ax,Az
	620	+call umfpack_dl_triplet_to_col(n,n,nz,wTi,wTj,T,Ap,Ai,A,status)
	621	+! if status is not zero a problem has occured
	622	+if (status /= 0) then
	623	+ write(,) "Problem during umfpack_dl_triplet_to_col"
	624	+endif
	625	+
	626	+! Define defaults control values
	627	+call umfpack_dl_defaults(control)
	628	+
	629	+! Symbolic analysis
	630	+call umfpack_dl_symbolic(n,n,Ap,Ai,A,symbolic, control, info)
	631	+! info(1) should be zero
	632	+if (info(1) /= 0) then
	633	+ write(,) "Problem during symbolic analysis"
	634	+ status=info(1)
	635	+endif
	636	+
	637	+! Numerical factorization
	638	+call umfpack_dl_numeric (Ap, Ai, A, symbolic, numeric, control, info)
	639	+! info(1) should be zero
	640	+if (info(1) /= 0) then
	641	+ write(,) "Problem during numerical factorization"
	642	+ status=info(1)
	643	+endif
	644	+
	645	+! free the C symbolic pointer
	646	+call umfpack_dl_free_symbolic (symbolic)
	647	+
	648	+call umfpack_dl_get_determinant(Mx,Ex,numeric,info,status)
	649	+! info(1) should be zero
	650	+if (info(1) /= 0) then
	651	+ write(,) "Problem during sparse determinant, returned code : ",info(1)
	652	+ status=info(1)
	653	+endif
	654	+det=Mx10*Ex
	655	+
	656	+! free the C numeric pointer
	657	+call umfpack_dl_free_numeric (numeric)
	658	+
	659	+deallocate(A)
	660	+deallocate(wTi,wTj)
	661	+end subroutine
	662	+
	663	+
	664	+subroutine fvn_di_sparse_det(n,nz,T,Ti,Tj,det,status)
	665	+implicit none
	666	+integer(kind=sp_kind), intent(in) :: n,nz
	667	+real(kind=dp_kind),dimension(nz),intent(in) :: T
	668	+integer(kind=sp_kind),dimension(nz),intent(in) :: Ti,Tj
	669	+integer(kind=sp_kind), intent(out) :: status
	670	+real(kind=dp_kind), intent(out) :: det
	671	+
	672	+integer(kind=sp_kind),dimension(:),allocatable :: wTi,wTj
	673	+real(kind=dp_kind),dimension(:),allocatable :: A
	674	+integer(kind=sp_kind),dimension(:),allocatable :: Ap,Ai
	675	+!integer(kind=dp_kind) :: symbolic,numeric
	676	+integer(kind=sp_kind),dimension(2) :: symbolic,numeric
	677	+! As symbolic and numeric are used to store a C pointer, it is necessary to
	678	+! still use an integer(kind=dp_kind) for 64bits machines
	679	+! An other possibility : integer(kind=sp_kind),dimension(2) :: symbolic,numeric
	680	+real(kind=dp_kind),dimension(90) :: info
	681	+real(kind=dp_kind),dimension(20) :: control
	682	+integer(kind=sp_kind) :: sys
	683	+real(kind=dp_kind) :: Mx,Ex
	684	+
	685	+
	686	+status=0
	687	+! we use a working copy of Ti and Tj to perform 1-based to 0-based translation
	688	+allocate(wTi(nz),wTj(nz))
	689	+wTi=Ti-1
	690	+wTj=Tj-1
	691	+allocate(A(nz))
	692	+allocate(Ap(n+1),Ai(nz))
	693	+
	694	+! perform the triplet to compressed column form -> Ap,Ai,Ax,Az
	695	+call umfpack_di_triplet_to_col(n,n,nz,wTi,wTj,T,Ap,Ai,A,status)
	696	+! if status is not zero a problem has occured
	697	+if (status /= 0) then
	698	+ write(,) "Problem during umfpack_di_triplet_to_col"
	699	+endif
	700	+
	701	+! Define defaults control values
	702	+call umfpack_di_defaults(control)
	703	+
	704	+! Symbolic analysis
	705	+call umfpack_di_symbolic(n,n,Ap,Ai,A,symbolic, control, info)
	706	+! info(1) should be zero
	707	+if (info(1) /= 0) then
	708	+ write(,) "Problem during symbolic analysis"
	709	+ status=info(1)
	710	+endif
	711	+
	712	+! Numerical factorization
	713	+call umfpack_di_numeric (Ap, Ai, A, symbolic, numeric, control, info)
	714	+! info(1) should be zero
	715	+if (info(1) /= 0) then
	716	+ write(,) "Problem during numerical factorization"
	717	+ status=info(1)
	718	+endif
	719	+
	720	+
	721	+! free the C symbolic pointer
	722	+call umfpack_di_free_symbolic (symbolic)
	723	+
	724	+call umfpack_di_get_determinant(Mx,Ex,numeric,info,status)
	725	+det=Mx10*Ex
	726	+! info(1) should be zero
	727	+if (info(1) /= 0) then
	728	+ write(,) "Problem during sparse determinant, returned code : ",info(1)
	729	+ status=info(1)
	730	+endif
	731	+
	732	+! free the C numeric pointer
	733	+call umfpack_di_free_numeric (numeric)
	734	+
	735	+deallocate(A)
	736	+deallocate(wTi,wTj)
	737	+end subroutine
	738	+
	739	+
	740	+
372	741
373	742
374	743	end module fvn_sparse
...	...	@@ -72,6 +72,11 @@
72	72	}
73	73
74	74
	75	+/* get determinant */
	76	+void umfpack_zl_get_determinant_ ( double Mx, double Mz, double Ex, void Numeric, double Info [UMFPACK_INFO], UF_long status )
	77	+{
	78	+ status = umfpack_zl_get_determinant (Mx,Mz,Ex,Numeric, Info);
	79	+}
75	80
76	81
77	82	/* complex(8) and integer(4) */
...	...	@@ -131,6 +136,14 @@
131	136	status = umfpack_zi_triplet_to_col (m, n, nz, Ti, Tj, Tx, Tz, Ap, Ai, Ax, Az, (int *) NULL);
132	137	}
133	138
	139	+
	140	+/* get determinant */
	141	+void umfpack_zi_get_determinant_ ( double Mx, double Mz, double Ex, void Numeric, double Info [UMFPACK_INFO], int status )
	142	+{
	143	+ status = umfpack_zi_get_determinant (Mx,Mz,Ex,Numeric, Info);
	144	+}
	145	+
	146	+
134	147	/* real(8) and integer(8) */
135	148
136	149	/* defaults */
...	...	@@ -178,6 +191,14 @@
178	191	status = umfpack_dl_triplet_to_col (m, n, nz, Ti, Tj, T, Ap, Ai, A, (UF_long *) NULL);
179	192	}
180	193
	194	+/* get determinant */
	195	+void umfpack_dl_get_determinant_ ( double Mx, double Ex, void *Numeric, double Info [UMFPACK_INFO], UF_long status )
	196	+{
	197	+ status = umfpack_dl_get_determinant (Mx,Ex,Numeric, Info);
	198	+}
	199	+
	200	+
	201	+
181	202	/* real(8) and integer(4) */
182	203
183	204	/* defaults */
...	...	@@ -223,5 +244,11 @@
223	244	double A [ ], int *status)
224	245	{
225	246	status = umfpack_di_triplet_to_col (m, n, nz, Ti, Tj, T, Ap, Ai, A, (int *) NULL);
	247	+}
	248	+
	249	+/* get determinant */
	250	+void umfpack_di_get_determinant_ ( double Mx, double Ex, void *Numeric, double Info [UMFPACK_INFO], int status )
	251	+{
	252	+ status = umfpack_di_get_determinant (Mx,Ex,Numeric, Info);
226	253	}
...	...	@@ -2,10 +2,11 @@
2	2	include $(BTREE)/Make.inc
3	3
4	4	programs = test_fac$(exext) test_matinv$(exext) test_specfunc$(exext) \
5		-test_det$(exext) test_matcon$(exext) test_matev$(exext) test_sparse$(exext) test_inter1d$(exext) \
	5	+test_det$(exext) test_matcon$(exext) test_matev$(exext) test_inter1d$(exext) \
6	6	test_inter2d$(exext) test_inter3d$(exext) test_akima$(exext) test_lsp$(exext) test_muller$(exext) \
7	7	test_integ$(exext) test_bsyn$(exext) test_bsjn$(exext) test_bskn$(exext) test_bsin$(exext) test_operators$(exext) test_ze1$(exext) \
8		-test_besri$(exext) test_besrj$(exext) test_bestime$(exext)
	8	+test_besri$(exext) test_besrj$(exext) test_bestime$(exext) \
	9	+test_sparse_zi$(exext) test_sparse_zl$(exext) test_sparse_di$(exext) test_sparse_dl$(exext)
9	10
10	11	prog:$(programs)
11	12
1		-program test_sparse
2		-use fvn_sparse
3		-implicit none
4		-integer(kind=ip_kind), parameter :: nz=12
5		-integer(kind=ip_kind), parameter :: n=5
6		-real(kind=dp_kind),dimension(nz) :: A
7		-real(kind=dp_kind),dimension(n,n) :: As
8		-integer(kind=ip_kind),dimension(nz) :: Ti,Tj
9		-real(kind=dp_kind),dimension(n) :: B,x
10		-integer(kind=ip_kind) :: status,i
11		-! Description of the matrix in triplet form
12		-A = (/ 2.,3.,3.,-1.,4.,4.,-3.,1.,2.,2.,6.,1. /)
13		-B = (/ 8., 45., -3., 3., 19./)
14		-Ti = (/ 1,2,1,3,5,2,3,4,5,3,2,5 /)
15		-Tj = (/ 1,1,2,2,2,3,3,3,3,4,5,5 /)
16		-
17		-! Reconstruction of the matrix in standard form
18		-! just needed for printing the matrix here
19		-As=0.
20		-do i=1,nz
21		- As(Ti(i),Tj(i))=A(i)
22		-end do
23		-write(,) "Matrix in standard representation :"
24		-do i=1,5
25		- write(*,'(5f8.4)') As(i,:)
26		-end do
27		-write(,)
28		-write(*,'("Right hand side :",5f8.4)') B
29		-
30		-!specific routine that will be used here
31		-!call fvn_di_sparse_solve(n,nz,A,Ti,Tj,B,x,status)
32		-call fvn_di_sparse_solve(n,nz,A,Ti,Tj,B,x,status)
33		-write(*,'("Solution :",5f8.4)') x
34		-write(*,'("Product matrix Solution :",5f8.4)') matmul(As,x)
35		-end program
	1	+program test_sparse
	2	+! test sparse routine di, real(8) and integer(4)
	3	+use fvn
	4	+implicit none
	5	+integer(kind=sp_kind), parameter :: nz=12
	6	+integer(kind=sp_kind), parameter :: n=5
	7	+real(kind=dp_kind),dimension(nz) :: A
	8	+real(kind=dp_kind),dimension(n,n) :: As
	9	+integer(kind=sp_kind),dimension(nz) :: Ti,Tj
	10	+real(kind=dp_kind),dimension(n) :: B,x
	11	+integer(kind=sp_kind) :: status,i
	12	+real(kind=dp_kind) :: det
	13	+! Description of the matrix in triplet form
	14	+A = (/ 2.,3.,3.,-1.,4.,4.,-3.,1.,2.,2.,6.,1. /)
	15	+B = (/ 8., 45., -3., 3., 19./)
	16	+Ti = (/ 1,2,1,3,5,2,3,4,5,3,2,5 /)
	17	+Tj = (/ 1,1,2,2,2,3,3,3,3,4,5,5 /)
	18	+
	19	+! Reconstruction of the matrix in standard form
	20	+! just needed for printing the matrix here
	21	+As=0.
	22	+do i=1,nz
	23	+ As(Ti(i),Tj(i))=A(i)
	24	+end do
	25	+write(,) "Matrix in standard representation :"
	26	+do i=1,5
	27	+ write(*,'(5f8.4)') As(i,:)
	28	+end do
	29	+write(,)
	30	+write(,) "Standard determinant =", fvn_det(5,As)
	31	+write(,)
	32	+write(*,'("Right hand side :",5f8.4)') B
	33	+
	34	+! can use either specific interface, fvn_di_sparse_det
	35	+! either generic one fvn_sparse_det
	36	+call fvn_di_sparse_det(n,nz,A,Ti,Tj,det,status)
	37	+write(,)
	38	+write(,) "Sparse Det = ",det
	39	+! can use either specific interface fvn_di_sparse_solve
	40	+! either generic one fvn_sparse_solve
	41	+! parameter det is optional
	42	+call fvn_di_sparse_solve(n,nz,A,Ti,Tj,B,x,status,det)
	43	+write(,)
	44	+write(,) "Sparse Det as solve option = ",det
	45	+write(,)
	46	+write(*,'("Solution :",5f8.4)') x
	47	+write(,)
	48	+write(*,'("Product matrix Solution :",5f8.4)') matmul(As,x)
	49	+end program
	1	+program test_sparse
	2	+! test sparse routine dl, real(8) and integer(8)
	3	+use fvn
	4	+implicit none
	5	+integer(kind=dp_kind), parameter :: nz=12
	6	+integer(kind=dp_kind), parameter :: n=5
	7	+real(kind=dp_kind),dimension(nz) :: A
	8	+real(kind=dp_kind),dimension(n,n) :: As
	9	+integer(kind=dp_kind),dimension(nz) :: Ti,Tj
	10	+real(kind=dp_kind),dimension(n) :: B,x
	11	+integer(kind=dp_kind) :: status,i
	12	+real(kind=dp_kind) :: det
	13	+! Description of the matrix in triplet form
	14	+A = (/ 2.,3.,3.,-1.,4.,4.,-3.,1.,2.,2.,6.,1. /)
	15	+B = (/ 8., 45., -3., 3., 19./)
	16	+Ti = (/ 1,2,1,3,5,2,3,4,5,3,2,5 /)
	17	+Tj = (/ 1,1,2,2,2,3,3,3,3,4,5,5 /)
	18	+
	19	+! Reconstruction of the matrix in standard form
	20	+! just needed for printing the matrix here
	21	+As=0.
	22	+do i=1,nz
	23	+ As(Ti(i),Tj(i))=A(i)
	24	+end do
	25	+write(,) "Matrix in standard representation :"
	26	+do i=1,5
	27	+ write(*,'(5f8.4)') As(i,:)
	28	+end do
	29	+write(,)
	30	+write(,) "Standard determinant =", fvn_det(5,As)
	31	+write(,)
	32	+write(*,'("Right hand side :",5f8.4)') B
	33	+
	34	+! can use either specific interface, fvn_dl_sparse_det
	35	+! either generic one fvn_sparse_det
	36	+call fvn_dl_sparse_det(n,nz,A,Ti,Tj,det,status)
	37	+write(,)
	38	+write(,) "Sparse Det = ",det
	39	+! can use either specific interface fvn_dl_sparse_solve
	40	+! either generic one fvn_sparse_solve
	41	+! parameter det is optional
	42	+call fvn_dl_sparse_solve(n,nz,A,Ti,Tj,B,x,status,det)
	43	+write(,)
	44	+write(,) "Sparse Det as solve option = ",det
	45	+write(,)
	46	+write(*,'("Solution :",5f8.4)') x
	47	+write(,)
	48	+write(*,'("Product matrix Solution :",5f8.4)') matmul(As,x)
	49	+end program
	1	+program test_sparse
	2	+! test sparse routine zl, complex(8) and integer(8)
	3	+use fvn
	4	+implicit none
	5	+integer(kind=dp_kind), parameter :: nz=12
	6	+integer(kind=dp_kind), parameter :: n=5
	7	+complex(kind=dp_kind),dimension(nz) :: A
	8	+complex(kind=dp_kind),dimension(n,n) :: As
	9	+integer(kind=dp_kind),dimension(nz) :: Ti,Tj
	10	+complex(kind=dp_kind),dimension(n) :: B,x
	11	+integer(kind=dp_kind) :: status,i
	12	+complex(kind=dp_kind) :: det
	13	+character(len=80) :: fmcmplx
	14	+
	15	+fmcmplx='(5("(",f8.5,",",f8.5,") "))'
	16	+
	17	+! Description of the matrix in triplet form
	18	+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.) /)
	19	+B = (/ (8.,3.), (45.,1.), (-3.,-2.), (3.,0.), (19.,2.) /)
	20	+Ti = (/ 1,2,1,3,5,2,3,4,5,3,2,5 /)
	21	+Tj = (/ 1,1,2,2,2,3,3,3,3,4,5,5 /)
	22	+
	23	+! Reconstruction of the matrix in standard form
	24	+As=0.
	25	+do i=1,nz
	26	+ As(Ti(i),Tj(i))=A(i)
	27	+end do
	28	+
	29	+write(,) "Matrix in standard representation :"
	30	+do i=1,5
	31	+ write(*,fmcmplx) As(i,:)
	32	+end do
	33	+write(,)
	34	+write(,) "Standard determinant : ",fvn_det(5,As)
	35	+write(,)
	36	+write(,) "Right hand side :"
	37	+write(*,fmcmplx) B
	38	+! can use either specific interface, fvn_zl_sparse_det
	39	+! either generic one fvn_sparse_det
	40	+call fvn_zl_sparse_det(n,nz,A,Ti,Tj,det,status)
	41	+write(,)
	42	+write(,) "Sparse Det = ",det
	43	+! can use either specific interface fvn_zl_sparse_solve
	44	+! either generic one fvn_sparse_solve
	45	+! parameter det is optional
	46	+call fvn_zl_sparse_solve(n,nz,A,Ti,Tj,B,x,status,det)
	47	+write(,)
	48	+write(,) "Sparse Det as solve option= ",det
	49	+write(,)
	50	+write(,) "Solution :"
	51	+write(*,fmcmplx) x
	52	+write(,)
	53	+write(,) "Product matrix Solution :"
	54	+write(*,fmcmplx) matmul(As,x)
	55	+end program