diff --git a/fvn_sparse/fvn_sparse.f90 b/fvn_sparse/fvn_sparse.f90 index 1fc1ad3..6e32712 100644 --- a/fvn_sparse/fvn_sparse.f90 +++ b/fvn_sparse/fvn_sparse.f90 @@ -41,7 +41,7 @@ integer(kind=dp_kind),dimension(nz),intent(in) :: Ti,Tj 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 +real(kind=dp_kind), dimension(3), optional, intent(out) :: det integer(kind=dp_kind),dimension(:),allocatable :: wTi,wTj real(kind=dp_kind),dimension(:),allocatable :: Tx,Tz @@ -52,7 +52,6 @@ real(kind=dp_kind),dimension(:),allocatable :: xx,xz,bx,bz 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 @@ -99,8 +98,7 @@ 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 + call umfpack_zl_get_determinant(det(1),det(2),det(3),numeric,info,status) ! info(1) should be zero if (info(1) /= 0) then write(*,*) "Problem during sparse determinant, returned code : ",info(1) @@ -148,7 +146,7 @@ integer(kind=sp_kind),dimension(nz),intent(in) :: Ti,Tj 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 +real(kind=dp_kind), dimension(3), optional, intent(out) :: det integer(kind=sp_kind),dimension(:),allocatable :: wTi,wTj real(kind=dp_kind),dimension(:),allocatable :: Tx,Tz @@ -163,7 +161,6 @@ real(kind=dp_kind),dimension(:),allocatable :: xx,xz,bx,bz 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 @@ -208,8 +205,7 @@ 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 + call umfpack_zi_get_determinant(det(1),det(2),det(3),numeric,info,status) ! info(1) should be zero if (info(1) /= 0) then write(*,*) "Problem during sparse determinant, returned code : ",info(1) @@ -258,7 +254,7 @@ integer(kind=dp_kind),dimension(nz),intent(in) :: Ti,Tj 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 +real(kind=dp_kind), dimension(2), optional, intent(out) :: det integer(kind=dp_kind),dimension(:),allocatable :: wTi,wTj real(kind=dp_kind),dimension(:),allocatable :: A @@ -268,7 +264,6 @@ integer(kind=dp_kind) :: symbolic,numeric 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 @@ -309,8 +304,7 @@ 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 + call umfpack_dl_get_determinant(det(1),det(2),numeric,info,status) ! info(1) should be zero if (info(1) /= 0) then write(*,*) "Problem during sparse determinant, returned code : ",info(1) @@ -350,7 +344,7 @@ integer(kind=sp_kind),dimension(nz),intent(in) :: Ti,Tj 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 +real(kind=dp_kind), dimension(2), optional, intent(out) :: det integer(kind=sp_kind),dimension(:),allocatable :: wTi,wTj real(kind=dp_kind),dimension(:),allocatable :: A @@ -363,7 +357,6 @@ 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 @@ -404,8 +397,7 @@ 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 + call umfpack_di_get_determinant(det(1),det(2),numeric,info,status) ! info(1) should be zero if (info(1) /= 0) then write(*,*) "Problem during sparse determinant, returned code : ",info(1) @@ -441,7 +433,7 @@ 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 +real(kind=dp_kind), dimension(3), intent(out) :: det integer(kind=dp_kind), intent(out) :: status integer(kind=dp_kind),dimension(:),allocatable :: wTi,wTj @@ -496,13 +488,12 @@ endif ! free the C symbolic pointer call umfpack_zl_free_symbolic (symbolic) -call umfpack_zl_get_determinant(Mx,Mz,Ex,numeric,info,status) +call umfpack_zl_get_determinant(det(1),det(2),det(3),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) @@ -519,7 +510,7 @@ 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 +real(kind=dp_kind), dimension(3), intent(out) :: det integer(kind=sp_kind),dimension(:),allocatable :: wTi,wTj real(kind=dp_kind),dimension(:),allocatable :: Tx,Tz @@ -532,7 +523,6 @@ integer(kind=sp_kind),dimension(2) :: symbolic,numeric ! 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 @@ -575,13 +565,12 @@ endif ! free the C symbolic pointer call umfpack_zi_free_symbolic (symbolic) -call umfpack_zi_get_determinant(Mx,Mz,Ex,numeric,info,status) +call umfpack_zi_get_determinant(det(1),det(2),det(3),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) @@ -597,7 +586,7 @@ 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 +real(kind=dp_kind), dimension(2), intent(out) :: det integer(kind=dp_kind),dimension(:),allocatable :: wTi,wTj real(kind=dp_kind),dimension(:),allocatable :: A @@ -606,7 +595,6 @@ 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,Ex status=0 ! we use a working copy of Ti and Tj to perform 1-based to 0-based translation @@ -645,13 +633,12 @@ endif ! free the C symbolic pointer call umfpack_dl_free_symbolic (symbolic) -call umfpack_dl_get_determinant(Mx,Ex,numeric,info,status) +call umfpack_dl_get_determinant(det(1),det(2),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) @@ -667,7 +654,7 @@ 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 +real(kind=dp_kind), dimension(2), intent(out) :: det integer(kind=sp_kind),dimension(:),allocatable :: wTi,wTj real(kind=dp_kind),dimension(:),allocatable :: A @@ -680,7 +667,6 @@ 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 @@ -721,8 +707,7 @@ 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 +call umfpack_di_get_determinant(det(1),det(2),numeric,info,status) ! info(1) should be zero if (info(1) /= 0) then write(*,*) "Problem during sparse determinant, returned code : ",info(1) diff --git a/fvn_test/test_sparse_di.f90 b/fvn_test/test_sparse_di.f90 index 868446b..2c2fe0e 100644 --- a/fvn_test/test_sparse_di.f90 +++ b/fvn_test/test_sparse_di.f90 @@ -9,7 +9,7 @@ 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 +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./) @@ -35,13 +35,13 @@ write(*,'("Right hand side :",5f8.4)') B ! either generic one fvn_sparse_det call fvn_di_sparse_det(n,nz,A,Ti,Tj,det,status) write(*,*) -write(*,*) "Sparse Det = ",det +write(*,*) "Sparse Det = ",det(1)*10**det(2) ! 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(*,*) "Sparse Det as solve option = ",det(1)*10**det(2) write(*,*) write(*,'("Solution :",5f8.4)') x write(*,*) diff --git a/fvn_test/test_sparse_dl.f90 b/fvn_test/test_sparse_dl.f90 index c683661..d53c10d 100644 --- a/fvn_test/test_sparse_dl.f90 +++ b/fvn_test/test_sparse_dl.f90 @@ -9,7 +9,7 @@ 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 +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./) @@ -35,13 +35,13 @@ write(*,'("Right hand side :",5f8.4)') B ! either generic one fvn_sparse_det call fvn_dl_sparse_det(n,nz,A,Ti,Tj,det,status) write(*,*) -write(*,*) "Sparse Det = ",det +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 +write(*,*) "Sparse Det as solve option = ",det(1)*10**det(2) write(*,*) write(*,'("Solution :",5f8.4)') x write(*,*) diff --git a/fvn_test/test_sparse_zi.f90 b/fvn_test/test_sparse_zi.f90 index 1d5ac67..58d6f8e 100644 --- a/fvn_test/test_sparse_zi.f90 +++ b/fvn_test/test_sparse_zi.f90 @@ -8,7 +8,7 @@ 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 +real(kind=dp_kind),dimension(3) :: det character(len=80) :: fmcmplx fmcmplx='(5("(",f8.5,",",f8.5,") "))' @@ -39,13 +39,13 @@ write(*,fmcmplx) B ! either generic one fvn_sparse_det call fvn_zi_sparse_det(n,nz,A,Ti,Tj,det,status) write(*,*) -write(*,*) "Sparse Det = ",det +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,A,Ti,Tj,B,x,status,det) write(*,*) -write(*,*) "Sparse Det as solve option= ",det +write(*,*) "Sparse Det as solve option= ",cmplx(det(1),det(2),kind=dp_kind)*10**det(3) write(*,*) write(*,*) "Solution :" write(*,fmcmplx) x diff --git a/fvn_test/test_sparse_zl.f90 b/fvn_test/test_sparse_zl.f90 index 24e7867..0038012 100644 --- a/fvn_test/test_sparse_zl.f90 +++ b/fvn_test/test_sparse_zl.f90 @@ -9,7 +9,7 @@ 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 +real(kind=dp_kind),dimension(3) :: det character(len=80) :: fmcmplx fmcmplx='(5("(",f8.5,",",f8.5,") "))' @@ -39,13 +39,13 @@ write(*,fmcmplx) B ! either generic one fvn_sparse_det call fvn_zl_sparse_det(n,nz,A,Ti,Tj,det,status) write(*,*) -write(*,*) "Sparse Det = ",det +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,A,Ti,Tj,B,x,status,det) write(*,*) -write(*,*) "Sparse Det as solve option= ",det +write(*,*) "Sparse Det as solve option= ",cmplx(det(1),det(2),kind=dp_kind)*10**det(3) write(*,*) write(*,*) "Solution :" write(*,fmcmplx) x