lxsd / fvn

Download as
Browse Code »

Commit db6500dcfb6c82b1630b7a6d35067cd1ee90cf78

Authored by William Daniau
1 parent b5f099f3cf
Exists in multi

First draft for a multi resolution

Showing 1 changed file with 83 additions and 0 deletions Inline Diff

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