lxsd / fvn

Download as
Browse Code »

Commit 14c923613ee917ad1024b77d5e354f4fa6ddc8ff

Authored by wdaniau
1 parent aee8aaaf18
Exists in master and in 2 other branches geevx_2020, multi

Return explicit umfpack return code when necessary 

git-svn-id: https://lxsd.femto-st.fr/svn/fvn@73 b657c933-2333-4658-acf2-d3c7c2708721

Showing 1 changed file with 30 additions and 29 deletions Inline Diff

fvn_sparse/fvn_sparse.f90
Diff comments   View file @ 14c9236
module fvn_sparse 1 1 module fvn_sparse
use fvn_common 2 2 use fvn_common
implicit none 3 3 implicit none
4 4
! Sparse solving 5 5 ! Sparse solving
interface fvn_sparse_solve 6 6 interface fvn_sparse_solve
module procedure fvn_zl_sparse_solve,fvn_zi_sparse_solve,fvn_dl_sparse_solve,fvn_di_sparse_solve 7 7 module procedure fvn_zl_sparse_solve,fvn_zi_sparse_solve,fvn_dl_sparse_solve,fvn_di_sparse_solve
end interface fvn_sparse_solve 8 8 end interface fvn_sparse_solve
9 9
interface fvn_sparse_det 10 10 interface fvn_sparse_det
module procedure fvn_zl_sparse_det,fvn_zi_sparse_det,fvn_dl_sparse_det,fvn_di_sparse_det 11 11 module procedure fvn_zl_sparse_det,fvn_zi_sparse_det,fvn_dl_sparse_det,fvn_di_sparse_det
end interface fvn_sparse_det 12 12 end interface fvn_sparse_det
contains 13 13 contains
14 14
function umfpack_return_code(c) 15 15 function umfpack_return_code(c)
implicit none 16 16 implicit none
integer(4), intent(in) :: c 17 17 integer(kind=sp_kind), intent(in) :: c
character(len=80) :: umfpack_return_code 18 18 character(len=80) :: umfpack_return_code
select case(c) 19 19 select case(c)
case(0) 20 20 case(0)
umfpack_return_code="UMFPACK_OK" 21 21 umfpack_return_code="UMFPACK_OK"
case(1) 22 22 case(1)
umfpack_return_code="UMFPACK_WARNING_singular_matrix" 23 23 umfpack_return_code="UMFPACK_WARNING_singular_matrix"
case(2) 24 24 case(2)
umfpack_return_code="UMFPACK_WARNING_determinant_underflow" 25 25 umfpack_return_code="UMFPACK_WARNING_determinant_underflow"
case(3) 26 26 case(3)
umfpack_return_code="UMFPACK_WARNING_determinant_overflow" 27 27 umfpack_return_code="UMFPACK_WARNING_determinant_overflow"
case(-1) 28 28 case(-1)
umfpack_return_code="UMFPACK_ERROR_out_of_memory" 29 29 umfpack_return_code="UMFPACK_ERROR_out_of_memory"
case(-3) 30 30 case(-3)
umfpack_return_code="UMFPACK_ERROR_invalid_Numeric_object" 31 31 umfpack_return_code="UMFPACK_ERROR_invalid_Numeric_object"
case(-4) 32 32 case(-4)
umfpack_return_code="UMFPACK_ERROR_invalid_Symbolic_object" 33 33 umfpack_return_code="UMFPACK_ERROR_invalid_Symbolic_object"
case(-5) 34 34 case(-5)
umfpack_return_code="UMFPACK_ERROR_argument_missing" 35 35 umfpack_return_code="UMFPACK_ERROR_argument_missing"
case(-6) 36 36 case(-6)
umfpack_return_code="UMFPACK_ERROR_n_nonpositive" 37 37 umfpack_return_code="UMFPACK_ERROR_n_nonpositive"
case(-8) 38 38 case(-8)
umfpack_return_code="UMFPACK_ERROR_invalid_matrix" 39 39 umfpack_return_code="UMFPACK_ERROR_invalid_matrix"
case(-11) 40 40 case(-11)
umfpack_return_code="UMFPACK_ERROR_different_pattern" 41 41 umfpack_return_code="UMFPACK_ERROR_different_pattern"
case(-13) 42 42 case(-13)
umfpack_return_code="UMFPACK_ERROR_invalid_system" 43 43 umfpack_return_code="UMFPACK_ERROR_invalid_system"
case(-15) 44 44 case(-15)
umfpack_return_code="UMFPACK_ERROR_invalid_permutation" 45 45 umfpack_return_code="UMFPACK_ERROR_invalid_permutation"
case(-911) 46 46 case(-911)
umfpack_return_code="UMFPACK_ERROR_internal_error" 47 47 umfpack_return_code="UMFPACK_ERROR_internal_error"
case(-17) 48 48 case(-17)
umfpack_return_code="UMFPACK_ERROR_file_IO" 49 49 umfpack_return_code="UMFPACK_ERROR_file_IO"
case default 50 50 case default
umfpack_return_code="Unknown return code" 51 51 umfpack_return_code="Unknown return code"
end select 52 52 end select
end function 53 53 end function
54
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 54 55 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! 55 56 !
! SPARSE RESOLUTION 56 57 ! SPARSE RESOLUTION
! 57 58 !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 58 59 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! 59 60 !
! Sparse resolution is done by interfaçing Tim Davi's UMFPACK 60 61 ! Sparse resolution is done by interfaçing Tim Davi's UMFPACK
! http://www.cise.ufl.edu/research/sparse/SuiteSparse/ 61 62 ! http://www.cise.ufl.edu/research/sparse/SuiteSparse/
! Used packages from SuiteSparse : AMD,UMFPACK,UFconfig 62 63 ! Used packages from SuiteSparse : AMD,UMFPACK,UFconfig
! 63 64 !
! Solve Ax=B using UMFPACK 64 65 ! Solve Ax=B using UMFPACK
! 65 66 !
! Where A is a sparse matrix given in its triplet form 66 67 ! Where A is a sparse matrix given in its triplet form
! T -> non zero elements 67 68 ! T -> non zero elements
! Ti,Tj -> row and column index (1-based) of the given elt 68 69 ! Ti,Tj -> row and column index (1-based) of the given elt
! n : rank of matrix A 69 70 ! n : rank of matrix A
! nz : number of non zero elts 70 71 ! nz : number of non zero elts
! 71 72 !
! fvn_*_sparse_solve 72 73 ! fvn_*_sparse_solve
! * = zl : double complex + integer(kind=dp_kind) 73 74 ! * = zl : double complex + integer(kind=dp_kind)
! * = zi : double complex + integer(kind=sp_kind) 74 75 ! * = zi : double complex + integer(kind=sp_kind)
! 75 76 !
subroutine fvn_zl_sparse_solve(n,nz,T,Ti,Tj,B,x,status,det) 76 77 subroutine fvn_zl_sparse_solve(n,nz,T,Ti,Tj,B,x,status,det)
implicit none 77 78 implicit none
integer(kind=dp_kind), intent(in) :: n,nz 78 79 integer(kind=dp_kind), intent(in) :: n,nz
complex(kind=dp_kind),dimension(nz),intent(in) :: T 79 80 complex(kind=dp_kind),dimension(nz),intent(in) :: T
integer(kind=dp_kind),dimension(nz),intent(in) :: Ti,Tj 80 81 integer(kind=dp_kind),dimension(nz),intent(in) :: Ti,Tj
complex(kind=dp_kind),dimension(n),intent(in) :: B 81 82 complex(kind=dp_kind),dimension(n),intent(in) :: B
complex(kind=dp_kind),dimension(n),intent(out) :: x 82 83 complex(kind=dp_kind),dimension(n),intent(out) :: x
integer(kind=dp_kind), intent(out) :: status 83 84 integer(kind=dp_kind), intent(out) :: status
real(kind=dp_kind), dimension(3), optional, intent(out) :: det 84 85 real(kind=dp_kind), dimension(3), optional, intent(out) :: det
85 86
integer(kind=dp_kind),dimension(:),allocatable :: wTi,wTj 86 87 integer(kind=dp_kind),dimension(:),allocatable :: wTi,wTj
real(kind=dp_kind),dimension(:),allocatable :: Tx,Tz 87 88 real(kind=dp_kind),dimension(:),allocatable :: Tx,Tz
real(kind=dp_kind),dimension(:),allocatable :: Ax,Az 88 89 real(kind=dp_kind),dimension(:),allocatable :: Ax,Az
integer(kind=dp_kind),dimension(:),allocatable :: Ap,Ai 89 90 integer(kind=dp_kind),dimension(:),allocatable :: Ap,Ai
integer(kind=dp_kind) :: symbolic,numeric 90 91 integer(kind=dp_kind) :: symbolic,numeric
real(kind=dp_kind),dimension(:),allocatable :: xx,xz,bx,bz 91 92 real(kind=dp_kind),dimension(:),allocatable :: xx,xz,bx,bz
real(kind=dp_kind),dimension(90) :: info 92 93 real(kind=dp_kind),dimension(90) :: info
real(kind=dp_kind),dimension(20) :: control 93 94 real(kind=dp_kind),dimension(20) :: control
integer(kind=dp_kind) :: sys 94 95 integer(kind=dp_kind) :: sys
95 96
96 97
status=0 97 98 status=0
98 99
! we use a working copy of Ti and Tj to perform 1-based to 0-based translation 99 100 ! 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 100 101 ! Tx and Tz are the real and imaginary parts of T
allocate(wTi(nz),wTj(nz)) 101 102 allocate(wTi(nz),wTj(nz))
allocate(Tx(nz),Tz(nz)) 102 103 allocate(Tx(nz),Tz(nz))
Tx=dble(T) 103 104 Tx=dble(T)
Tz=aimag(T) 104 105 Tz=aimag(T)
wTi=Ti-1 105 106 wTi=Ti-1
wTj=Tj-1 106 107 wTj=Tj-1
allocate(Ax(nz),Az(nz)) 107 108 allocate(Ax(nz),Az(nz))
allocate(Ap(n+1),Ai(nz)) 108 109 allocate(Ap(n+1),Ai(nz))
109 110
! perform the triplet to compressed column form -> Ap,Ai,Ax,Az 110 111 ! 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) 111 112 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 112 113 ! if status is not zero a problem has occured
if (status /= 0) then 113 114 if (status /= 0) then
write(*,*) "Problem during umfpack_zl_triplet_to_col" 114 115 write(*,*) "Problem during umfpack_zl_triplet_to_col"
endif 115 116 endif
116 117
! Define defaults control values 117 118 ! Define defaults control values
call umfpack_zl_defaults(control) 118 119 call umfpack_zl_defaults(control)
119 120
! Symbolic analysis 120 121 ! Symbolic analysis
call umfpack_zl_symbolic(n,n,Ap,Ai,Ax,Az,symbolic, control, info) 121 122 call umfpack_zl_symbolic(n,n,Ap,Ai,Ax,Az,symbolic, control, info)
! info(1) should be zero 122 123 ! info(1) should be zero
if (info(1) /= 0) then 123 124 if (info(1) /= 0) then
write(*,*) "Problem during symbolic analysis" 124 125 write(*,*) "Problem during symbolic analysis : ",trim(umfpack_return_code(int(info(1))))
status=info(1) 125 126 status=info(1)
endif 126 127 endif
127 128
! Numerical factorization 128 129 ! Numerical factorization
call umfpack_zl_numeric (Ap, Ai, Ax, Az, symbolic, numeric, control, info) 129 130 call umfpack_zl_numeric (Ap, Ai, Ax, Az, symbolic, numeric, control, info)
! info(1) should be zero 130 131 ! info(1) should be zero
if (info(1) /= 0) then 131 132 if (info(1) /= 0) then
write(*,*) "Problem during numerical factorization" 132 133 write(*,*) "Problem during numerical factorization : ",trim(umfpack_return_code(int(info(1))))
status=info(1) 133 134 status=info(1)
endif 134 135 endif
135 136
! free the C symbolic pointer 136 137 ! free the C symbolic pointer
call umfpack_zl_free_symbolic (symbolic) 137 138 call umfpack_zl_free_symbolic (symbolic)
138 139
! if parameter det is present, the determinant of the matrix is calculated 139 140 ! if parameter det is present, the determinant of the matrix is calculated
if (present(det) ) then 140 141 if (present(det) ) then
call umfpack_zl_get_determinant(det(1),det(2),det(3),numeric,info,status) 141 142 call umfpack_zl_get_determinant(det(1),det(2),det(3),numeric,info,status)
! info(1) should be zero 142 143 ! info(1) should be zero
if (info(1) /= 0) then 143 144 if (info(1) /= 0) then
if ( (info(1) < 1) .or. (info(1) >3) ) then ! not a warning 144 145 if ( (info(1) < 1) .or. (info(1) >3) ) then ! not a warning
write(*,*) "Problem during sparse determinant, returned code : ",info(1) 145 146 write(*,*) "Problem during sparse determinant : ",trim(umfpack_return_code(int(info(1))))
endif 146 147 endif
status=info(1) 147 148 status=info(1)
endif 148 149 endif
endif 149 150 endif
150 151
151 152
152 153
allocate(bx(n),bz(n),xx(n),xz(n)) 153 154 allocate(bx(n),bz(n),xx(n),xz(n))
bx=dble(B) 154 155 bx=dble(B)
bz=aimag(B) 155 156 bz=aimag(B)
sys=0 156 157 sys=0
! sys may be used to define type of solving -> see umfpack.h 157 158 ! sys may be used to define type of solving -> see umfpack.h
158 159
! Solving 159 160 ! Solving
call umfpack_zl_solve (sys, Ap, Ai, Ax,Az, xx,xz, bx,bz, numeric, control, info) 160 161 call umfpack_zl_solve (sys, Ap, Ai, Ax,Az, xx,xz, bx,bz, numeric, control, info)
! info(1) should be zero 161 162 ! info(1) should be zero
if (info(1) /= 0) then 162 163 if (info(1) /= 0) then
write(*,*) "Problem during solving" 163 164 write(*,*) "Problem during solving : ",trim(umfpack_return_code(int(info(1))))
status=info(1) 164 165 status=info(1)
endif 165 166 endif
166 167
167 168
! free the C numeric pointer 168 169 ! free the C numeric pointer
call umfpack_zl_free_numeric (numeric) 169 170 call umfpack_zl_free_numeric (numeric)
170 171
x=cmplx(xx,xz,dp_kind) 171 172 x=cmplx(xx,xz,dp_kind)
172 173
deallocate(bx,bz,xx,xz) 173 174 deallocate(bx,bz,xx,xz)
deallocate(Ax,Az) 174 175 deallocate(Ax,Az)
deallocate(Tx,Tz) 175 176 deallocate(Tx,Tz)
deallocate(wTi,wTj) 176 177 deallocate(wTi,wTj)
end subroutine 177 178 end subroutine
178 179
179 180
180 181
181 182
182 183
subroutine fvn_zi_sparse_solve(n,nz,T,Ti,Tj,B,x,status,det) 183 184 subroutine fvn_zi_sparse_solve(n,nz,T,Ti,Tj,B,x,status,det)
implicit none 184 185 implicit none
integer(kind=sp_kind), intent(in) :: n,nz 185 186 integer(kind=sp_kind), intent(in) :: n,nz
complex(kind=dp_kind),dimension(nz),intent(in) :: T 186 187 complex(kind=dp_kind),dimension(nz),intent(in) :: T
integer(kind=sp_kind),dimension(nz),intent(in) :: Ti,Tj 187 188 integer(kind=sp_kind),dimension(nz),intent(in) :: Ti,Tj
complex(kind=dp_kind),dimension(n),intent(in) :: B 188 189 complex(kind=dp_kind),dimension(n),intent(in) :: B
complex(kind=dp_kind),dimension(n),intent(out) :: x 189 190 complex(kind=dp_kind),dimension(n),intent(out) :: x
integer(kind=sp_kind), intent(out) :: status 190 191 integer(kind=sp_kind), intent(out) :: status
real(kind=dp_kind), dimension(3), optional, intent(out) :: det 191 192 real(kind=dp_kind), dimension(3), optional, intent(out) :: det
192 193
integer(kind=sp_kind),dimension(:),allocatable :: wTi,wTj 193 194 integer(kind=sp_kind),dimension(:),allocatable :: wTi,wTj
real(kind=dp_kind),dimension(:),allocatable :: Tx,Tz 194 195 real(kind=dp_kind),dimension(:),allocatable :: Tx,Tz
real(kind=dp_kind),dimension(:),allocatable :: Ax,Az 195 196 real(kind=dp_kind),dimension(:),allocatable :: Ax,Az
integer(kind=sp_kind),dimension(:),allocatable :: Ap,Ai 196 197 integer(kind=sp_kind),dimension(:),allocatable :: Ap,Ai
!integer(kind=dp_kind) :: symbolic,numeric 197 198 !integer(kind=dp_kind) :: symbolic,numeric
integer(kind=sp_kind),dimension(2) :: symbolic,numeric 198 199 integer(kind=sp_kind),dimension(2) :: symbolic,numeric
! As symbolic and numeric are used to store a C pointer, it is necessary to 199 200 ! 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 200 201 ! still use an integer(kind=dp_kind) for 64bits machines
! An other possibility : integer(kind=sp_kind),dimension(2) :: symbolic,numeric 201 202 ! An other possibility : integer(kind=sp_kind),dimension(2) :: symbolic,numeric
real(kind=dp_kind),dimension(:),allocatable :: xx,xz,bx,bz 202 203 real(kind=dp_kind),dimension(:),allocatable :: xx,xz,bx,bz
real(kind=dp_kind),dimension(90) :: info 203 204 real(kind=dp_kind),dimension(90) :: info
real(kind=dp_kind),dimension(20) :: control 204 205 real(kind=dp_kind),dimension(20) :: control
integer(kind=sp_kind) :: sys 205 206 integer(kind=sp_kind) :: sys
206 207
status=0 207 208 status=0
! we use a working copy of Ti and Tj to perform 1-based to 0-based translation 208 209 ! 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 209 210 ! Tx and Tz are the real and imaginary parts of T
allocate(wTi(nz),wTj(nz)) 210 211 allocate(wTi(nz),wTj(nz))
allocate(Tx(nz),Tz(nz)) 211 212 allocate(Tx(nz),Tz(nz))
Tx=dble(T) 212 213 Tx=dble(T)
Tz=aimag(T) 213 214 Tz=aimag(T)
wTi=Ti-1 214 215 wTi=Ti-1
wTj=Tj-1 215 216 wTj=Tj-1
allocate(Ax(nz),Az(nz)) 216 217 allocate(Ax(nz),Az(nz))
allocate(Ap(n+1),Ai(nz)) 217 218 allocate(Ap(n+1),Ai(nz))
218 219
! perform the triplet to compressed column form -> Ap,Ai,Ax,Az 219 220 ! 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) 220 221 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 221 222 ! if status is not zero a problem has occured
if (status /= 0) then 222 223 if (status /= 0) then
write(*,*) "Problem during umfpack_zl_triplet_to_col" 223 224 write(*,*) "Problem during umfpack_zl_triplet_to_col"
endif 224 225 endif
225 226
! Define defaults control values 226 227 ! Define defaults control values
call umfpack_zi_defaults(control) 227 228 call umfpack_zi_defaults(control)
228 229
! Symbolic analysis 229 230 ! Symbolic analysis
call umfpack_zi_symbolic(n,n,Ap,Ai,Ax,Az,symbolic, control, info) 230 231 call umfpack_zi_symbolic(n,n,Ap,Ai,Ax,Az,symbolic, control, info)
! info(1) should be zero 231 232 ! info(1) should be zero
if (info(1) /= 0) then 232 233 if (info(1) /= 0) then
write(*,*) "Problem during symbolic analysis" 233 234 write(*,*) "Problem during symbolic analysis : ",trim(umfpack_return_code(int(info(1))))
status=info(1) 234 235 status=info(1)
endif 235 236 endif
236 237
! Numerical factorization 237 238 ! Numerical factorization
call umfpack_zi_numeric (Ap, Ai, Ax, Az, symbolic, numeric, control, info) 238 239 call umfpack_zi_numeric (Ap, Ai, Ax, Az, symbolic, numeric, control, info)
! info(1) should be zero 239 240 ! info(1) should be zero
if (info(1) /= 0) then 240 241 if (info(1) /= 0) then
write(*,*) "Problem during numerical factorization" 241 242 write(*,*) "Problem during numerical factorization : ",trim(umfpack_return_code(int(info(1))))
status=info(1) 242 243 status=info(1)
endif 243 244 endif
244 245
! free the C symbolic pointer 245 246 ! free the C symbolic pointer
call umfpack_zi_free_symbolic (symbolic) 246 247 call umfpack_zi_free_symbolic (symbolic)
247 248
! if parameter det is present, the determinant of the matrix is calculated 248 249 ! if parameter det is present, the determinant of the matrix is calculated
if (present(det) ) then 249 250 if (present(det) ) then
call umfpack_zi_get_determinant(det(1),det(2),det(3),numeric,info,status) 250 251 call umfpack_zi_get_determinant(det(1),det(2),det(3),numeric,info,status)
! info(1) should be zero 251 252 ! info(1) should be zero
if (info(1) /= 0) then 252 253 if (info(1) /= 0) then
if ( (info(1) < 1) .or. (info(1) >3) ) then ! not a warning 253 254 if ( (info(1) < 1) .or. (info(1) >3) ) then ! not a warning
write(*,*) "Problem during sparse determinant, returned code : ",info(1) 254 255 write(*,*) "Problem during sparse determinant : ",trim(umfpack_return_code(int(info(1))))
endif 255 256 endif
status=info(1) 256 257 status=info(1)
endif 257 258 endif
endif 258 259 endif
259 260
260 261
261 262
262 263
allocate(bx(n),bz(n),xx(n),xz(n)) 263 264 allocate(bx(n),bz(n),xx(n),xz(n))
bx=dble(B) 264 265 bx=dble(B)
bz=aimag(B) 265 266 bz=aimag(B)
sys=0 266 267 sys=0
! sys may be used to define type of solving -> see umfpack.h 267 268 ! sys may be used to define type of solving -> see umfpack.h
268 269
! Solving 269 270 ! Solving
call umfpack_zi_solve (sys, Ap, Ai, Ax,Az, xx,xz, bx,bz, numeric, control, info) 270 271 call umfpack_zi_solve (sys, Ap, Ai, Ax,Az, xx,xz, bx,bz, numeric, control, info)
! info(1) should be zero 271 272 ! info(1) should be zero
if (info(1) /= 0) then 272 273 if (info(1) /= 0) then
write(*,*) "Problem during solving" 273 274 write(*,*) "Problem during solving : ",trim(umfpack_return_code(int(info(1))))
status=info(1) 274 275 status=info(1)
endif 275 276 endif
276 277
! free the C numeric pointer 277 278 ! free the C numeric pointer
call umfpack_zi_free_numeric (numeric) 278 279 call umfpack_zi_free_numeric (numeric)
279 280
x=cmplx(xx,xz,dp_kind) 280 281 x=cmplx(xx,xz,dp_kind)
281 282
deallocate(bx,bz,xx,xz) 282 283 deallocate(bx,bz,xx,xz)
deallocate(Ax,Az) 283 284 deallocate(Ax,Az)
deallocate(Tx,Tz) 284 285 deallocate(Tx,Tz)
deallocate(wTi,wTj) 285 286 deallocate(wTi,wTj)
end subroutine 286 287 end subroutine
287 288
288 289
289 290
290 291
291 292
292 293
subroutine fvn_dl_sparse_solve(n,nz,T,Ti,Tj,B,x,status,det) 293 294 subroutine fvn_dl_sparse_solve(n,nz,T,Ti,Tj,B,x,status,det)
implicit none 294 295 implicit none
integer(kind=dp_kind), intent(in) :: n,nz 295 296 integer(kind=dp_kind), intent(in) :: n,nz
real(kind=dp_kind),dimension(nz),intent(in) :: T 296 297 real(kind=dp_kind),dimension(nz),intent(in) :: T
integer(kind=dp_kind),dimension(nz),intent(in) :: Ti,Tj 297 298 integer(kind=dp_kind),dimension(nz),intent(in) :: Ti,Tj
real(kind=dp_kind),dimension(n),intent(in) :: B 298 299 real(kind=dp_kind),dimension(n),intent(in) :: B
real(kind=dp_kind),dimension(n),intent(out) :: x 299 300 real(kind=dp_kind),dimension(n),intent(out) :: x
integer(kind=dp_kind), intent(out) :: status 300 301 integer(kind=dp_kind), intent(out) :: status
real(kind=dp_kind), dimension(2), optional, intent(out) :: det 301 302 real(kind=dp_kind), dimension(2), optional, intent(out) :: det
302 303
integer(kind=dp_kind),dimension(:),allocatable :: wTi,wTj 303 304 integer(kind=dp_kind),dimension(:),allocatable :: wTi,wTj
real(kind=dp_kind),dimension(:),allocatable :: A 304 305 real(kind=dp_kind),dimension(:),allocatable :: A
integer(kind=dp_kind),dimension(:),allocatable :: Ap,Ai 305 306 integer(kind=dp_kind),dimension(:),allocatable :: Ap,Ai
!integer(kind=dp_kind) :: symbolic,numeric 306 307 !integer(kind=dp_kind) :: symbolic,numeric
integer(kind=dp_kind) :: symbolic,numeric 307 308 integer(kind=dp_kind) :: symbolic,numeric
real(kind=dp_kind),dimension(90) :: info 308 309 real(kind=dp_kind),dimension(90) :: info
real(kind=dp_kind),dimension(20) :: control 309 310 real(kind=dp_kind),dimension(20) :: control
integer(kind=dp_kind) :: sys 310 311 integer(kind=dp_kind) :: sys
311 312
status=0 312 313 status=0
! we use a working copy of Ti and Tj to perform 1-based to 0-based translation 313 314 ! we use a working copy of Ti and Tj to perform 1-based to 0-based translation
allocate(wTi(nz),wTj(nz)) 314 315 allocate(wTi(nz),wTj(nz))
wTi=Ti-1 315 316 wTi=Ti-1
wTj=Tj-1 316 317 wTj=Tj-1
allocate(A(nz)) 317 318 allocate(A(nz))
allocate(Ap(n+1),Ai(nz)) 318 319 allocate(Ap(n+1),Ai(nz))
319 320
! perform the triplet to compressed column form -> Ap,Ai,Ax,Az 320 321 ! 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) 321 322 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 322 323 ! if status is not zero a problem has occured
if (status /= 0) then 323 324 if (status /= 0) then
write(*,*) "Problem during umfpack_dl_triplet_to_col" 324 325 write(*,*) "Problem during umfpack_dl_triplet_to_col"
endif 325 326 endif
326 327
! Define defaults control values 327 328 ! Define defaults control values
call umfpack_dl_defaults(control) 328 329 call umfpack_dl_defaults(control)
329 330
! Symbolic analysis 330 331 ! Symbolic analysis
call umfpack_dl_symbolic(n,n,Ap,Ai,A,symbolic, control, info) 331 332 call umfpack_dl_symbolic(n,n,Ap,Ai,A,symbolic, control, info)
! info(1) should be zero 332 333 ! info(1) should be zero
if (info(1) /= 0) then 333 334 if (info(1) /= 0) then
write(*,*) "Problem during symbolic analysis" 334 335 write(*,*) "Problem during symbolic analysis : ",trim(umfpack_return_code(int(info(1))))
status=info(1) 335 336 status=info(1)
endif 336 337 endif
337 338
! Numerical factorization 338 339 ! Numerical factorization
call umfpack_dl_numeric (Ap, Ai, A, symbolic, numeric, control, info) 339 340 call umfpack_dl_numeric (Ap, Ai, A, symbolic, numeric, control, info)
! info(1) should be zero 340 341 ! info(1) should be zero
if (info(1) /= 0) then 341 342 if (info(1) /= 0) then
write(*,*) "Problem during numerical factorization" 342 343 write(*,*) "Problem during numerical factorization : ",trim(umfpack_return_code(int(info(1))))
status=info(1) 343 344 status=info(1)
endif 344 345 endif
345 346
! free the C symbolic pointer 346 347 ! free the C symbolic pointer
call umfpack_dl_free_symbolic (symbolic) 347 348 call umfpack_dl_free_symbolic (symbolic)
348 349
! if parameter det is present, the determinant of the matrix is calculated 349 350 ! if parameter det is present, the determinant of the matrix is calculated
if (present(det) ) then 350 351 if (present(det) ) then
call umfpack_dl_get_determinant(det(1),det(2),numeric,info,status) 351 352 call umfpack_dl_get_determinant(det(1),det(2),numeric,info,status)
! info(1) should be zero 352 353 ! info(1) should be zero
if (info(1) /= 0) then 353 354 if (info(1) /= 0) then
if ( (info(1) < 1) .or. (info(1) >3) ) then ! not a warning 354 355 if ( (info(1) < 1) .or. (info(1) >3) ) then ! not a warning
write(*,*) "Problem during sparse determinant, returned code : ",info(1) 355 356 write(*,*) "Problem during sparse determinant : ",trim(umfpack_return_code(int(info(1))))
endif 356 357 endif
status=info(1) 357 358 status=info(1)
endif 358 359 endif
endif 359 360 endif
360 361
361 362
sys=0 362 363 sys=0
! sys may be used to define type of solving -> see umfpack.h 363 364 ! sys may be used to define type of solving -> see umfpack.h
364 365
! Solving 365 366 ! Solving
call umfpack_dl_solve (sys, Ap, Ai, A, x, B, numeric, control, info) 366 367 call umfpack_dl_solve (sys, Ap, Ai, A, x, B, numeric, control, info)
! info(1) should be zero 367 368 ! info(1) should be zero
if (info(1) /= 0) then 368 369 if (info(1) /= 0) then
write(*,*) "Problem during solving" 369 370 write(*,*) "Problem during solving : ",trim(umfpack_return_code(int(info(1))))
status=info(1) 370 371 status=info(1)
endif 371 372 endif
372 373
! free the C numeric pointer 373 374 ! free the C numeric pointer
call umfpack_dl_free_numeric (numeric) 374 375 call umfpack_dl_free_numeric (numeric)
375 376
deallocate(A) 376 377 deallocate(A)
deallocate(wTi,wTj) 377 378 deallocate(wTi,wTj)
end subroutine 378 379 end subroutine
379 380
380 381
381 382
382 383
383 384
384 385
subroutine fvn_di_sparse_solve(n,nz,T,Ti,Tj,B,x,status,det) 385 386 subroutine fvn_di_sparse_solve(n,nz,T,Ti,Tj,B,x,status,det)
implicit none 386 387 implicit none
integer(kind=sp_kind), intent(in) :: n,nz 387 388 integer(kind=sp_kind), intent(in) :: n,nz
real(kind=dp_kind),dimension(nz),intent(in) :: T 388 389 real(kind=dp_kind),dimension(nz),intent(in) :: T
integer(kind=sp_kind),dimension(nz),intent(in) :: Ti,Tj 389 390 integer(kind=sp_kind),dimension(nz),intent(in) :: Ti,Tj
real(kind=dp_kind),dimension(n),intent(in) :: B 390 391 real(kind=dp_kind),dimension(n),intent(in) :: B
real(kind=dp_kind),dimension(n),intent(out) :: x 391 392 real(kind=dp_kind),dimension(n),intent(out) :: x
integer(kind=sp_kind), intent(out) :: status 392 393 integer(kind=sp_kind), intent(out) :: status
real(kind=dp_kind), dimension(2), optional, intent(out) :: det 393 394 real(kind=dp_kind), dimension(2), optional, intent(out) :: det
394 395
integer(kind=sp_kind),dimension(:),allocatable :: wTi,wTj 395 396 integer(kind=sp_kind),dimension(:),allocatable :: wTi,wTj
real(kind=dp_kind),dimension(:),allocatable :: A 396 397 real(kind=dp_kind),dimension(:),allocatable :: A
integer(kind=sp_kind),dimension(:),allocatable :: Ap,Ai 397 398 integer(kind=sp_kind),dimension(:),allocatable :: Ap,Ai
!integer(kind=dp_kind) :: symbolic,numeric 398 399 !integer(kind=dp_kind) :: symbolic,numeric
integer(kind=sp_kind),dimension(2) :: symbolic,numeric 399 400 integer(kind=sp_kind),dimension(2) :: symbolic,numeric
! As symbolic and numeric are used to store a C pointer, it is necessary to 400 401 ! 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 401 402 ! still use an integer(kind=dp_kind) for 64bits machines
! An other possibility : integer(kind=sp_kind),dimension(2) :: symbolic,numeric 402 403 ! An other possibility : integer(kind=sp_kind),dimension(2) :: symbolic,numeric
real(kind=dp_kind),dimension(90) :: info 403 404 real(kind=dp_kind),dimension(90) :: info
real(kind=dp_kind),dimension(20) :: control 404 405 real(kind=dp_kind),dimension(20) :: control
integer(kind=sp_kind) :: sys 405 406 integer(kind=sp_kind) :: sys
406 407
status=0 407 408 status=0
! we use a working copy of Ti and Tj to perform 1-based to 0-based translation 408 409 ! we use a working copy of Ti and Tj to perform 1-based to 0-based translation
allocate(wTi(nz),wTj(nz)) 409 410 allocate(wTi(nz),wTj(nz))
wTi=Ti-1 410 411 wTi=Ti-1
wTj=Tj-1 411 412 wTj=Tj-1
allocate(A(nz)) 412 413 allocate(A(nz))
allocate(Ap(n+1),Ai(nz)) 413 414 allocate(Ap(n+1),Ai(nz))
414 415
! perform the triplet to compressed column form -> Ap,Ai,Ax,Az 415 416 ! 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) 416 417 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 417 418 ! if status is not zero a problem has occured
if (status /= 0) then 418 419 if (status /= 0) then
write(*,*) "Problem during umfpack_di_triplet_to_col" 419 420 write(*,*) "Problem during umfpack_di_triplet_to_col"
endif 420 421 endif
421 422
! Define defaults control values 422 423 ! Define defaults control values
call umfpack_di_defaults(control) 423 424 call umfpack_di_defaults(control)
424 425
! Symbolic analysis 425 426 ! Symbolic analysis
call umfpack_di_symbolic(n,n,Ap,Ai,A,symbolic, control, info) 426 427 call umfpack_di_symbolic(n,n,Ap,Ai,A,symbolic, control, info)
! info(1) should be zero 427 428 ! info(1) should be zero
if (info(1) /= 0) then 428 429 if (info(1) /= 0) then
write(*,*) "Problem during symbolic analysis" 429 430 write(*,*) "Problem during symbolic analysis : ",trim(umfpack_return_code(int(info(1))))
status=info(1) 430 431 status=info(1)
endif 431 432 endif
432 433
! Numerical factorization 433 434 ! Numerical factorization
call umfpack_di_numeric (Ap, Ai, A, symbolic, numeric, control, info) 434 435 call umfpack_di_numeric (Ap, Ai, A, symbolic, numeric, control, info)
! info(1) should be zero 435 436 ! info(1) should be zero
if (info(1) /= 0) then 436 437 if (info(1) /= 0) then
write(*,*) "Problem during numerical factorization" 437 438 write(*,*) "Problem during numerical factorization : ",trim(umfpack_return_code(int(info(1))))
status=info(1) 438 439 status=info(1)
endif 439 440 endif
440 441
! free the C symbolic pointer 441 442 ! free the C symbolic pointer
call umfpack_di_free_symbolic (symbolic) 442 443 call umfpack_di_free_symbolic (symbolic)
443 444
! if parameter det is present, the determinant of the matrix is calculated 444 445 ! if parameter det is present, the determinant of the matrix is calculated
if (present(det) ) then 445 446 if (present(det) ) then
call umfpack_di_get_determinant(det(1),det(2),numeric,info,status) 446 447 call umfpack_di_get_determinant(det(1),det(2),numeric,info,status)
! info(1) should be zero 447 448 ! info(1) should be zero
if (info(1) /= 0) then 448 449 if (info(1) /= 0) then
if ( (info(1) < 1) .or. (info(1) >3) ) then ! not a warning 449 450 if ( (info(1) < 1) .or. (info(1) >3) ) then ! not a warning
write(*,*) "Problem during sparse determinant, returned code : ",info(1) 450 451 write(*,*) "Problem during sparse determinant : ",trim(umfpack_return_code(int(info(1))))
endif 451 452 endif
status=info(1) 452 453 status=info(1)
endif 453 454 endif
endif 454 455 endif
455 456
sys=0 456 457 sys=0
! sys may be used to define type of solving -> see umfpack.h 457 458 ! sys may be used to define type of solving -> see umfpack.h
! Solving 458 459 ! Solving
call umfpack_di_solve (sys, Ap, Ai, A, x, B, numeric, control, info) 459 460 call umfpack_di_solve (sys, Ap, Ai, A, x, B, numeric, control, info)
! info(1) should be zero 460 461 ! info(1) should be zero
if (info(1) /= 0) then 461 462 if (info(1) /= 0) then
write(*,*) "Problem during solving" 462 463 write(*,*) "Problem during solving : ",trim(umfpack_return_code(int(info(1))))
status=info(1) 463 464 status=info(1)
endif 464 465 endif
465 466
! free the C numeric pointer 466 467 ! free the C numeric pointer
call umfpack_di_free_numeric (numeric) 467 468 call umfpack_di_free_numeric (numeric)
468 469
deallocate(A) 469 470 deallocate(A)
deallocate(wTi,wTj) 470 471 deallocate(wTi,wTj)
end subroutine 471 472 end subroutine
472 473
473 474
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 474 475 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! 475 476 !
! SPARSE DETERMINANT 476 477 ! SPARSE DETERMINANT
! 477 478 !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 478 479 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine fvn_zl_sparse_det(n,nz,T,Ti,Tj,det,status) 479 480 subroutine fvn_zl_sparse_det(n,nz,T,Ti,Tj,det,status)
implicit none 480 481 implicit none
integer(kind=dp_kind), intent(in) :: n,nz 481 482 integer(kind=dp_kind), intent(in) :: n,nz
complex(kind=dp_kind),dimension(nz),intent(in) :: T 482 483 complex(kind=dp_kind),dimension(nz),intent(in) :: T
integer(kind=dp_kind),dimension(nz),intent(in) :: Ti,Tj 483 484 integer(kind=dp_kind),dimension(nz),intent(in) :: Ti,Tj
real(kind=dp_kind), dimension(3), intent(out) :: det 484 485 real(kind=dp_kind), dimension(3), intent(out) :: det
integer(kind=dp_kind), intent(out) :: status 485 486 integer(kind=dp_kind), intent(out) :: status
486 487
integer(kind=dp_kind),dimension(:),allocatable :: wTi,wTj 487 488 integer(kind=dp_kind),dimension(:),allocatable :: wTi,wTj
real(kind=dp_kind),dimension(:),allocatable :: Tx,Tz 488 489 real(kind=dp_kind),dimension(:),allocatable :: Tx,Tz
real(kind=dp_kind),dimension(:),allocatable :: Ax,Az 489 490 real(kind=dp_kind),dimension(:),allocatable :: Ax,Az
integer(kind=dp_kind),dimension(:),allocatable :: Ap,Ai 490 491 integer(kind=dp_kind),dimension(:),allocatable :: Ap,Ai
integer(kind=dp_kind) :: symbolic,numeric 491 492 integer(kind=dp_kind) :: symbolic,numeric
real(kind=dp_kind),dimension(90) :: info 492 493 real(kind=dp_kind),dimension(90) :: info
real(kind=dp_kind),dimension(20) :: control 493 494 real(kind=dp_kind),dimension(20) :: control
real(kind=dp_kind) :: Mx,Mz,Ex 494 495 real(kind=dp_kind) :: Mx,Mz,Ex
495 496
496 497
status=0 497 498 status=0
498 499
! we use a working copy of Ti and Tj to perform 1-based to 0-based translation 499 500 ! 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 500 501 ! Tx and Tz are the real and imaginary parts of T
allocate(wTi(nz),wTj(nz)) 501 502 allocate(wTi(nz),wTj(nz))
allocate(Tx(nz),Tz(nz)) 502 503 allocate(Tx(nz),Tz(nz))
Tx=dble(T) 503 504 Tx=dble(T)
Tz=aimag(T) 504 505 Tz=aimag(T)
wTi=Ti-1 505 506 wTi=Ti-1
wTj=Tj-1 506 507 wTj=Tj-1
allocate(Ax(nz),Az(nz)) 507 508 allocate(Ax(nz),Az(nz))
allocate(Ap(n+1),Ai(nz)) 508 509 allocate(Ap(n+1),Ai(nz))
509 510
! perform the triplet to compressed column form -> Ap,Ai,Ax,Az 510 511 ! 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) 511 512 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 512 513 ! if status is not zero a problem has occured
if (status /= 0) then 513 514 if (status /= 0) then
write(*,*) "Problem during umfpack_zl_triplet_to_col" 514 515 write(*,*) "Problem during umfpack_zl_triplet_to_col"
endif 515 516 endif
516 517
! Define defaults control values 517 518 ! Define defaults control values
call umfpack_zl_defaults(control) 518 519 call umfpack_zl_defaults(control)
519 520
! Symbolic analysis 520 521 ! Symbolic analysis
call umfpack_zl_symbolic(n,n,Ap,Ai,Ax,Az,symbolic, control, info) 521 522 call umfpack_zl_symbolic(n,n,Ap,Ai,Ax,Az,symbolic, control, info)
! info(1) should be zero 522 523 ! info(1) should be zero
if (info(1) /= 0) then 523 524 if (info(1) /= 0) then
write(*,*) "Problem during symbolic analysis" 524 525 write(*,*) "Problem during symbolic analysis : ",trim(umfpack_return_code(int(info(1))))
status=info(1) 525 526 status=info(1)
endif 526 527 endif
527 528
! Numerical factorization 528 529 ! Numerical factorization
call umfpack_zl_numeric (Ap, Ai, Ax, Az, symbolic, numeric, control, info) 529 530 call umfpack_zl_numeric (Ap, Ai, Ax, Az, symbolic, numeric, control, info)
! info(1) should be zero 530 531 ! info(1) should be zero
if (info(1) /= 0) then 531 532 if (info(1) /= 0) then
write(*,*) "Problem during numerical factorization" 532 533 write(*,*) "Problem during numerical factorization : ",trim(umfpack_return_code(int(info(1))))
status=info(1) 533 534 status=info(1)
endif 534 535 endif
535 536
! free the C symbolic pointer 536 537 ! free the C symbolic pointer
call umfpack_zl_free_symbolic (symbolic) 537 538 call umfpack_zl_free_symbolic (symbolic)
538 539
call umfpack_zl_get_determinant(det(1),det(2),det(3),numeric,info,status) 539 540 call umfpack_zl_get_determinant(det(1),det(2),det(3),numeric,info,status)
! info(1) should be zero 540 541 ! info(1) should be zero
if (info(1) /= 0) then 541 542 if (info(1) /= 0) then
if ( (info(1) < 1) .or. (info(1) >3) ) then ! not a warning 542 543 if ( (info(1) < 1) .or. (info(1) >3) ) then ! not a warning
write(*,*) "Problem during sparse determinant, returned code : ",info(1) 543 544 write(*,*) "Problem during sparse determinant : ",trim(umfpack_return_code(int(info(1))))
endif 544 545 endif
status=info(1) 545 546 status=info(1)
endif 546 547 endif
547 548
! free the C numeric pointer 548 549 ! free the C numeric pointer
call umfpack_zl_free_numeric (numeric) 549 550 call umfpack_zl_free_numeric (numeric)
550 551
deallocate(Ax,Az) 551 552 deallocate(Ax,Az)
deallocate(Tx,Tz) 552 553 deallocate(Tx,Tz)
deallocate(wTi,wTj) 553 554 deallocate(wTi,wTj)
end subroutine 554 555 end subroutine
555 556
556 557
subroutine fvn_zi_sparse_det(n,nz,T,Ti,Tj,det,status) 557 558 subroutine fvn_zi_sparse_det(n,nz,T,Ti,Tj,det,status)
implicit none 558 559 implicit none
integer(kind=sp_kind), intent(in) :: n,nz 559 560 integer(kind=sp_kind), intent(in) :: n,nz
complex(kind=dp_kind),dimension(nz),intent(in) :: T 560 561 complex(kind=dp_kind),dimension(nz),intent(in) :: T
integer(kind=sp_kind),dimension(nz),intent(in) :: Ti,Tj 561 562 integer(kind=sp_kind),dimension(nz),intent(in) :: Ti,Tj
integer(kind=sp_kind), intent(out) :: status 562 563 integer(kind=sp_kind), intent(out) :: status
real(kind=dp_kind), dimension(3), intent(out) :: det 563 564 real(kind=dp_kind), dimension(3), intent(out) :: det
564 565
integer(kind=sp_kind),dimension(:),allocatable :: wTi,wTj 565 566 integer(kind=sp_kind),dimension(:),allocatable :: wTi,wTj
real(kind=dp_kind),dimension(:),allocatable :: Tx,Tz 566 567 real(kind=dp_kind),dimension(:),allocatable :: Tx,Tz
real(kind=dp_kind),dimension(:),allocatable :: Ax,Az 567 568 real(kind=dp_kind),dimension(:),allocatable :: Ax,Az
integer(kind=sp_kind),dimension(:),allocatable :: Ap,Ai 568 569 integer(kind=sp_kind),dimension(:),allocatable :: Ap,Ai
!integer(kind=dp_kind) :: symbolic,numeric 569 570 !integer(kind=dp_kind) :: symbolic,numeric
integer(kind=sp_kind),dimension(2) :: symbolic,numeric 570 571 integer(kind=sp_kind),dimension(2) :: symbolic,numeric
! As symbolic and numeric are used to store a C pointer, it is necessary to 571 572 ! 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 572 573 ! still use an integer(kind=dp_kind) for 64bits machines
! An other possibility : integer(kind=sp_kind),dimension(2) :: symbolic,numeric 573 574 ! An other possibility : integer(kind=sp_kind),dimension(2) :: symbolic,numeric
real(kind=dp_kind),dimension(90) :: info 574 575 real(kind=dp_kind),dimension(90) :: info
real(kind=dp_kind),dimension(20) :: control 575 576 real(kind=dp_kind),dimension(20) :: control
576 577
status=0 577 578 status=0
! we use a working copy of Ti and Tj to perform 1-based to 0-based translation 578 579 ! 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 579 580 ! Tx and Tz are the real and imaginary parts of T
allocate(wTi(nz),wTj(nz)) 580 581 allocate(wTi(nz),wTj(nz))
allocate(Tx(nz),Tz(nz)) 581 582 allocate(Tx(nz),Tz(nz))
Tx=dble(T) 582 583 Tx=dble(T)
Tz=aimag(T) 583 584 Tz=aimag(T)
wTi=Ti-1 584 585 wTi=Ti-1
wTj=Tj-1 585 586 wTj=Tj-1
allocate(Ax(nz),Az(nz)) 586 587 allocate(Ax(nz),Az(nz))
allocate(Ap(n+1),Ai(nz)) 587 588 allocate(Ap(n+1),Ai(nz))
588 589
! perform the triplet to compressed column form -> Ap,Ai,Ax,Az 589 590 ! 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) 590 591 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 591 592 ! if status is not zero a problem has occured
if (status /= 0) then 592 593 if (status /= 0) then
write(*,*) "Problem during umfpack_zl_triplet_to_col" 593 594 write(*,*) "Problem during umfpack_zl_triplet_to_col"
endif 594 595 endif
595 596
! Define defaults control values 596 597 ! Define defaults control values
call umfpack_zi_defaults(control) 597 598 call umfpack_zi_defaults(control)
598 599
! Symbolic analysis 599 600 ! Symbolic analysis
call umfpack_zi_symbolic(n,n,Ap,Ai,Ax,Az,symbolic, control, info) 600 601 call umfpack_zi_symbolic(n,n,Ap,Ai,Ax,Az,symbolic, control, info)
! info(1) should be zero 601 602 ! info(1) should be zero
if (info(1) /= 0) then 602 603 if (info(1) /= 0) then
write(*,*) "Problem during symbolic analysis" 603 604 write(*,*) "Problem during symbolic analysis : ",trim(umfpack_return_code(int(info(1))))
status=info(1) 604 605 status=info(1)
endif 605 606 endif
606 607
! Numerical factorization 607 608 ! Numerical factorization
call umfpack_zi_numeric (Ap, Ai, Ax, Az, symbolic, numeric, control, info) 608 609 call umfpack_zi_numeric (Ap, Ai, Ax, Az, symbolic, numeric, control, info)
! info(1) should be zero 609 610 ! info(1) should be zero
if (info(1) /= 0) then 610 611 if (info(1) /= 0) then
write(*,*) "Problem during numerical factorization" 611 612 write(*,*) "Problem during numerical factorization : ",trim(umfpack_return_code(int(info(1))))
status=info(1) 612 613 status=info(1)
endif 613 614 endif
614 615
! free the C symbolic pointer 615 616 ! free the C symbolic pointer
call umfpack_zi_free_symbolic (symbolic) 616 617 call umfpack_zi_free_symbolic (symbolic)
617 618
call umfpack_zi_get_determinant(det(1),det(2),det(3),numeric,info,status) 618 619 call umfpack_zi_get_determinant(det(1),det(2),det(3),numeric,info,status)
! info(1) should be zero 619 620 ! info(1) should be zero
if (info(1) /= 0) then 620 621 if (info(1) /= 0) then
if ( (info(1) < 1) .or. (info(1) >3) ) then ! not a warning 621 622 if ( (info(1) < 1) .or. (info(1) >3) ) then ! not a warning
write(*,*) "Problem during sparse determinant, returned code : ",info(1) 622 623 write(*,*) "Problem during sparse determinant : ",trim(umfpack_return_code(int(info(1))))
endif 623 624 endif
status=info(1) 624 625 status=info(1)
endif 625 626 endif
626 627
! free the C numeric pointer 627 628 ! free the C numeric pointer
call umfpack_zi_free_numeric (numeric) 628 629 call umfpack_zi_free_numeric (numeric)
629 630
deallocate(Ax,Az) 630 631 deallocate(Ax,Az)
deallocate(Tx,Tz) 631 632 deallocate(Tx,Tz)
deallocate(wTi,wTj) 632 633 deallocate(wTi,wTj)
end subroutine 633 634 end subroutine
634 635
subroutine fvn_dl_sparse_det(n,nz,T,Ti,Tj,det,status) 635 636 subroutine fvn_dl_sparse_det(n,nz,T,Ti,Tj,det,status)
implicit none 636 637 implicit none
integer(kind=dp_kind), intent(in) :: n,nz 637 638 integer(kind=dp_kind), intent(in) :: n,nz
real(kind=dp_kind),dimension(nz),intent(in) :: T 638 639 real(kind=dp_kind),dimension(nz),intent(in) :: T
integer(kind=dp_kind),dimension(nz),intent(in) :: Ti,Tj 639 640 integer(kind=dp_kind),dimension(nz),intent(in) :: Ti,Tj
integer(kind=dp_kind), intent(out) :: status 640 641 integer(kind=dp_kind), intent(out) :: status
real(kind=dp_kind), dimension(2), intent(out) :: det 641 642 real(kind=dp_kind), dimension(2), intent(out) :: det
642 643
integer(kind=dp_kind),dimension(:),allocatable :: wTi,wTj 643 644 integer(kind=dp_kind),dimension(:),allocatable :: wTi,wTj
real(kind=dp_kind),dimension(:),allocatable :: A 644 645 real(kind=dp_kind),dimension(:),allocatable :: A
integer(kind=dp_kind),dimension(:),allocatable :: Ap,Ai 645 646 integer(kind=dp_kind),dimension(:),allocatable :: Ap,Ai
!integer(kind=dp_kind) :: symbolic,numeric 646 647 !integer(kind=dp_kind) :: symbolic,numeric
integer(kind=dp_kind) :: symbolic,numeric 647 648 integer(kind=dp_kind) :: symbolic,numeric
real(kind=dp_kind),dimension(90) :: info 648 649 real(kind=dp_kind),dimension(90) :: info
real(kind=dp_kind),dimension(20) :: control 649 650 real(kind=dp_kind),dimension(20) :: control
650 651
status=0 651 652 status=0
! we use a working copy of Ti and Tj to perform 1-based to 0-based translation 652 653 ! we use a working copy of Ti and Tj to perform 1-based to 0-based translation
allocate(wTi(nz),wTj(nz)) 653 654 allocate(wTi(nz),wTj(nz))
wTi=Ti-1 654 655 wTi=Ti-1
wTj=Tj-1 655 656 wTj=Tj-1
allocate(A(nz)) 656 657 allocate(A(nz))
allocate(Ap(n+1),Ai(nz)) 657 658 allocate(Ap(n+1),Ai(nz))
658 659
! perform the triplet to compressed column form -> Ap,Ai,Ax,Az 659 660 ! 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) 660 661 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 661 662 ! if status is not zero a problem has occured
if (status /= 0) then 662 663 if (status /= 0) then
write(*,*) "Problem during umfpack_dl_triplet_to_col" 663 664 write(*,*) "Problem during umfpack_dl_triplet_to_col"
endif 664 665 endif
665 666
! Define defaults control values 666 667 ! Define defaults control values
call umfpack_dl_defaults(control) 667 668 call umfpack_dl_defaults(control)
668 669
! Symbolic analysis 669 670 ! Symbolic analysis
call umfpack_dl_symbolic(n,n,Ap,Ai,A,symbolic, control, info) 670 671 call umfpack_dl_symbolic(n,n,Ap,Ai,A,symbolic, control, info)
! info(1) should be zero 671 672 ! info(1) should be zero
if (info(1) /= 0) then 672 673 if (info(1) /= 0) then
write(*,*) "Problem during symbolic analysis" 673 674 write(*,*) "Problem during symbolic analysis : ",trim(umfpack_return_code(int(info(1))))
status=info(1) 674 675 status=info(1)
endif 675 676 endif
676 677
! Numerical factorization 677 678 ! Numerical factorization
call umfpack_dl_numeric (Ap, Ai, A, symbolic, numeric, control, info) 678 679 call umfpack_dl_numeric (Ap, Ai, A, symbolic, numeric, control, info)
! info(1) should be zero 679 680 ! info(1) should be zero