lxsd / fvn

Download as
Browse Code ยป

Commit 4bd2122857ff6d0fa38dcdb157799912272ea76c

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

Changed sparse routine interface to avoid duplicating datas in the routines, 

which causes problems with larges systems.
    1) Complex inputs must be splitted into 2 arrays : real and imaginary parts
    2) Ti and Tj must now be  0-based

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

Showing 7 changed files with 2189 additions and 2218 deletions Side-by-side Diff

No preview for this file type
... ... @@ -283,17 +283,20 @@
283 283 \subsubsection{Sparse solving}
284 284 The provided routines solves the equation $Ax=B$ where A is sparse and given in its triplet form.
285 285  
  286 +Note that interface has changed from previous versions, splitting complex inputs into real and imaginary part and using 0-based indices, to fit umfpack routines inputs and avoid duplicating datas inside the routines.
  287 +
286 288 \begin{verbatim}
287 289 Module : fvn_sparse
288 290 call fvn_sparse_solve(n,nz,T,Ti,Tj,B,x,status,det)
  291 +call fvn_sparse_solve(n,nz,Tx,Tz,Ti,Tj,Bx,Bz,x,status,det)
289 292 \end{verbatim}
290 293 \begin{itemize}
291 294 \item For this family of subroutines the two letters (zl,zi,dl,di) of the specific interface name decribe the arguments's type. z is for complex(8), d for real(8), l for integer(8) and i for integer(4)
292 295 \item \texttt{n} (in) is an integer equal to the matrix rank
293 296 \item \texttt{nz} (in) is an integer equal to the number of non-zero elements
294   - \item \texttt{T(nz)} (in) is a complex/real array containing the non-zero elements
295   - \item \textt{Ti(nz)},\texttt{Tj(nz)} (in) are the indexes of the corresponding element of \texttt{T} in the original matrix.
296   - \item \texttt{B(n)} (in) is a complex/real array containing the second member of the equation.
  297 + \item \texttt{T(nz) or Tx(nz),Tz(nz)} (in) is a real array (or two real arrays for real/imaginary in case of a complex system) containing the non-zero elements
  298 + \item \textt{Ti(nz)},\texttt{Tj(nz)} (in) are the indexes of the corresponding element of \texttt{T} in the original matrix, this has to be 0-based as in C.
  299 + \item \texttt{B(n) or Bx(n),Bz(n)} (in) is a real array or two real arrays for real/imaginary in case of a complex system) containing the second member of the equation.
297 300 \item \texttt{x(n)} (out) is a complex/real array containing the solution
298 301 \item \texttt{status} (out) is an integer which contain non-zero is something went wrong
299 302 \item \texttt{det} (out), is an optional real(8) array of dimension 2 for dl and di specific interface (real systems) and dimension 3 for zl and zi interface (complex systems)
300 303  
... ... @@ -310,8 +313,10 @@
310 313 integer(8), parameter :: nz=12
311 314 integer(8), parameter :: n=5
312 315 complex(8),dimension(nz) :: A
  316 + real(8), dimension(nz) :: Ax,Az
313 317 integer(8),dimension(nz) :: Ti,Tj
314 318 complex(8),dimension(n) :: B,x
  319 + real(8), dimension(n) :: Bx,Bz
315 320 integer(8) :: status
316 321  
317 322 A = (/ (2.,0.),(3.,0.),(3.,0.),(-1.,0.),(4.,0.),(4.,0.),(-3.,0.),&
318 323  
319 324  
... ... @@ -319,10 +324,18 @@
319 324 B = (/ (8.,0.), (45.,0.), (-3.,0.), (3.,0.), (19.,0.)/)
320 325 Ti = (/ 1,2,1,3,5,2,3,4,5,3,2,5 /)
321 326 Tj = (/ 1,1,2,2,2,3,3,3,3,4,5,5 /)
  327 + Ax=real(A)
  328 + Az=aimag(A)
  329 + Bx=real(B)
  330 + Bz=aimag(B)
322 331  
  332 + ! 1-based to 0-based translation
  333 + Ti=Ti-1
  334 + Tj=Tj-1
  335 +
323 336 !specific routine that will be used here
324   - !call fvn_zl_sparse_solve(n,nz,A,Ti,Tj,B,x,status)
325   - call fvn_sparse_solve(n,nz,A,Ti,Tj,B,x,status)
  337 + !call fvn_zl_sparse_solve(n,nz,Ax,Az,Ti,Tj,Bx,Bz,x,status)
  338 + call fvn_sparse_solve(n,nz,Ax,Az,Ti,Tj,Bx,Bz,x,status)
326 339 write(*,*) x
327 340  
328 341 end program
... ... @@ -345,6 +358,10 @@
345 358 Ti = (/ 1,2,1,3,5,2,3,4,5,3,2,5 /)
346 359 Tj = (/ 1,1,2,2,2,3,3,3,3,4,5,5 /)
347 360  
  361 +! 1-based to 0-based translation
  362 +Ti=Ti-1
  363 +Tj=Tj-1
  364 +
348 365 !specific routine that will be used here
349 366 !call fvn_di_sparse_solve(n,nz,A,Ti,Tj,B,x,status)
350 367 call fvn_sparse_solve(n,nz,A,Ti,Tj,B,x,status)
351 368  
... ... @@ -362,14 +379,15 @@
362 379 \begin{verbatim}
363 380 Module : fvn_sparse
364 381 call fvn_sparse_det(n,nz,T,Ti,Tj,det,status)
  382 +call fvn_sparse_det(n,nz,Tx,Tz,Ti,Tj,det,status)
365 383 \end{verbatim}
366 384  
367 385 \begin{itemize}
368 386 \item For this family of subroutines the two letters (zl,zi,dl,di) of the specific interface name decribe the arguments's type. z is for complex(8), d for real(8), l for integer(8) and i for integer(4)
369 387 \item \texttt{n} (in) is an integer equal to the matrix rank
370 388 \item \texttt{nz} (in) is an integer equal to the number of non-zero elements
371   - \item \texttt{T(nz)} (in) is a complex/real array containing the non-zero elements
372   - \item \textt{Ti(nz)},\texttt{Tj(nz)} (in) are the indexes of the corresponding element of \texttt{T} in the original matrix.
  389 + \item \texttt{T(nz) or Tx(nz),Tz(nz)} (in) is a real array (or two real arrays for real/imaginary in case of a complex system) containing the non-zero elements
  390 + \item \textt{Ti(nz)},\texttt{Tj(nz)} (in) are the indexes of the corresponding element of \texttt{T} in the original matrix, it has to be 0-based as in C.
373 391 \item \texttt{det} (out), a real(8) array of dimension 2 for dl and di specific interface (real systems) and dimension 3 for zl and zi interface (complex systems)
374 392 \item \texttt{status} (out) is an integer which contain non-zero is something went wrong
375 393 \end{itemize}
376 394  
... ... @@ -395,9 +413,11 @@
395 413 integer(kind=sp_kind), parameter :: nz=12
396 414 integer(kind=sp_kind), parameter :: n=5
397 415 complex(kind=dp_kind),dimension(nz) :: A
  416 +real(kind=dp_kind), dimension(nz) :: Ax,Az
398 417 complex(kind=dp_kind),dimension(n,n) :: As
399 418 integer(kind=sp_kind),dimension(nz) :: Ti,Tj
400 419 complex(kind=dp_kind),dimension(n) :: B,x
  420 +real(kind=dp_kind), dimension(n) :: Bx,Bz
401 421 integer(kind=sp_kind) :: status,i
402 422 real(kind=dp_kind),dimension(3) :: det
403 423 character(len=80) :: fmcmplx
... ... @@ -410,6 +430,11 @@
410 430 Ti = (/ 1,2,1,3,5,2,3,4,5,3,2,5 /)
411 431 Tj = (/ 1,1,2,2,2,3,3,3,3,4,5,5 /)
412 432  
  433 +Ax=real(A)
  434 +Az=aimag(A)
  435 +Bx=real(B)
  436 +Bz=aimag(B)
  437 +
413 438 ! Reconstruction of the matrix in standard form
414 439 As=0.
415 440 do i=1,nz
416 441  
... ... @@ -428,13 +453,13 @@
428 453  
429 454 ! can use either specific interface, fvn_zi_sparse_det
430 455 ! either generic one fvn_sparse_det
431   -call fvn_zi_sparse_det(n,nz,A,Ti,Tj,det,status)
  456 +call fvn_zi_sparse_det(n,nz,Ax,Az,Ti,Tj,det,status)
432 457 write(*,*)
433 458 write(*,*) "Sparse Det = ",cmplx(det(1),det(2),kind=dp_kind)*10**det(3)
434 459 ! can use either specific interface fvn_zi_sparse_solve
435 460 ! either generic one fvn_sparse_solve
436 461 ! parameter det is optional
437   -call fvn_zi_sparse_solve(n,nz,A,Ti,Tj,B,x,status,det)
  462 +call fvn_zi_sparse_solve(n,nz,Ax,Az,Ti,Tj,Bx,Bz,x,status,det)
438 463 write(*,*)
439 464 write(*,*) "Sparse Det as solve option= ",cmplx(det(1),det(2),kind=dp_kind)*10**det(3)
440 465 write(*,*)
fvn_sparse/fvn_sparse.f90
Diff comments   View file @ 4bd2122
... ... @@ -65,8 +65,8 @@
65 65 ! Solve Ax=B using UMFPACK
66 66 !
67 67 ! Where A is a sparse matrix given in its triplet form
68   -! T -> non zero elements
69   -! Ti,Tj -> row and column index (1-based) of the given elt
  68 +! T -> non zero elements (Tx and Tz for real and imaginary part if complex)
  69 +! Ti,Tj -> row and column index (0-based) of the given elt
70 70 ! n : rank of matrix A
71 71 ! nz : number of non zero elts
72 72 !
73 73  
74 74  
75 75  
76 76  
... ... @@ -74,22 +74,20 @@
74 74 ! * = zl : double complex + integer(kind=dp_kind)
75 75 ! * = zi : double complex + integer(kind=sp_kind)
76 76 !
77   -subroutine fvn_zl_sparse_solve(n,nz,T,Ti,Tj,B,x,status,det)
  77 +subroutine fvn_zl_sparse_solve(n,nz,Tx,Tz,Ti,Tj,Bx,Bz,x,status,det)
78 78 implicit none
79 79 integer(kind=dp_kind), intent(in) :: n,nz
80   -complex(kind=dp_kind),dimension(nz),intent(in) :: T
  80 +real(kind=dp_kind),dimension(nz),intent(in) :: Tx,Tz
81 81 integer(kind=dp_kind),dimension(nz),intent(in) :: Ti,Tj
82   -complex(kind=dp_kind),dimension(n),intent(in) :: B
  82 +real(kind=dp_kind),dimension(n),intent(in) :: Bx,Bz
83 83 complex(kind=dp_kind),dimension(n),intent(out) :: x
84 84 integer(kind=dp_kind), intent(out) :: status
85 85 real(kind=dp_kind), dimension(3), optional, intent(out) :: det
86 86  
87   -integer(kind=dp_kind),dimension(:),allocatable :: wTi,wTj
88   -real(kind=dp_kind),dimension(:),allocatable :: Tx,Tz
89 87 real(kind=dp_kind),dimension(:),allocatable :: Ax,Az
90 88 integer(kind=dp_kind),dimension(:),allocatable :: Ap,Ai
91 89 integer(kind=dp_kind) :: symbolic,numeric
92   -real(kind=dp_kind),dimension(:),allocatable :: xx,xz,bx,bz
  90 +real(kind=dp_kind),dimension(:),allocatable :: xx,xz
93 91 real(kind=dp_kind),dimension(90) :: info
94 92 real(kind=dp_kind),dimension(20) :: control
95 93 integer(kind=dp_kind) :: sys
96 94  
... ... @@ -97,19 +95,11 @@
97 95  
98 96 status=0
99 97  
100   -! we use a working copy of Ti and Tj to perform 1-based to 0-based translation
101   -! Tx and Tz are the real and imaginary parts of T
102   -allocate(wTi(nz),wTj(nz))
103   -allocate(Tx(nz),Tz(nz))
104   -Tx=dble(T)
105   -Tz=aimag(T)
106   -wTi=Ti-1
107   -wTj=Tj-1
108 98 allocate(Ax(nz),Az(nz))
109 99 allocate(Ap(n+1),Ai(nz))
110 100  
111 101 ! perform the triplet to compressed column form -> Ap,Ai,Ax,Az
112   -call umfpack_zl_triplet_to_col(n,n,nz,wTi,wTj,Tx,Tz,Ap,Ai,Ax,Az,status)
  102 +call umfpack_zl_triplet_to_col(n,n,nz,Ti,Tj,Tx,Tz,Ap,Ai,Ax,Az,status)
113 103 ! if status is not zero a problem has occured
114 104 if (status /= 0) then
115 105 write(*,*) "Problem during umfpack_zl_triplet_to_col : ",trim(umfpack_return_code(int(status,kind=sp_kind)))
... ... @@ -151,9 +141,7 @@
151 141  
152 142  
153 143  
154   -allocate(bx(n),bz(n),xx(n),xz(n))
155   -bx=dble(B)
156   -bz=aimag(B)
  144 +allocate(xx(n),xz(n))
157 145 sys=0
158 146 ! sys may be used to define type of solving -> see umfpack.h
159 147  
160 148  
161 149  
162 150  
163 151  
164 152  
... ... @@ -171,28 +159,24 @@
171 159  
172 160 x=cmplx(xx,xz,dp_kind)
173 161  
174   -deallocate(bx,bz,xx,xz)
  162 +deallocate(xx,xz)
175 163 deallocate(Ax,Az)
176   -deallocate(Tx,Tz)
177   -deallocate(wTi,wTj)
178 164 end subroutine
179 165  
180 166  
181 167  
182 168  
183 169  
184   -subroutine fvn_zi_sparse_solve(n,nz,T,Ti,Tj,B,x,status,det)
  170 +subroutine fvn_zi_sparse_solve(n,nz,Tx,Tz,Ti,Tj,Bx,Bz,x,status,det)
185 171 implicit none
186 172 integer(kind=sp_kind), intent(in) :: n,nz
187   -complex(kind=dp_kind),dimension(nz),intent(in) :: T
  173 +real(kind=dp_kind),dimension(nz),intent(in) :: Tx,Tz
188 174 integer(kind=sp_kind),dimension(nz),intent(in) :: Ti,Tj
189   -complex(kind=dp_kind),dimension(n),intent(in) :: B
  175 +real(kind=dp_kind),dimension(n),intent(in) :: Bx,Bz
190 176 complex(kind=dp_kind),dimension(n),intent(out) :: x
191 177 integer(kind=sp_kind), intent(out) :: status
192 178 real(kind=dp_kind), dimension(3), optional, intent(out) :: det
193 179  
194   -integer(kind=sp_kind),dimension(:),allocatable :: wTi,wTj
195   -real(kind=dp_kind),dimension(:),allocatable :: Tx,Tz
196 180 real(kind=dp_kind),dimension(:),allocatable :: Ax,Az
197 181 integer(kind=sp_kind),dimension(:),allocatable :: Ap,Ai
198 182 !integer(kind=dp_kind) :: symbolic,numeric
199 183  
200 184  
... ... @@ -200,25 +184,17 @@
200 184 ! As symbolic and numeric are used to store a C pointer, it is necessary to
201 185 ! still use an integer(kind=dp_kind) for 64bits machines
202 186 ! An other possibility : integer(kind=sp_kind),dimension(2) :: symbolic,numeric
203   -real(kind=dp_kind),dimension(:),allocatable :: xx,xz,bx,bz
  187 +real(kind=dp_kind),dimension(:),allocatable :: xx,xz
204 188 real(kind=dp_kind),dimension(90) :: info
205 189 real(kind=dp_kind),dimension(20) :: control
206 190 integer(kind=sp_kind) :: sys
207 191  
208 192 status=0
209   -! we use a working copy of Ti and Tj to perform 1-based to 0-based translation
210   -! Tx and Tz are the real and imaginary parts of T
211   -allocate(wTi(nz),wTj(nz))
212   -allocate(Tx(nz),Tz(nz))
213   -Tx=dble(T)
214   -Tz=aimag(T)
215   -wTi=Ti-1
216   -wTj=Tj-1
217 193 allocate(Ax(nz),Az(nz))
218 194 allocate(Ap(n+1),Ai(nz))
219 195  
220 196 ! perform the triplet to compressed column form -> Ap,Ai,Ax,Az
221   -call umfpack_zi_triplet_to_col(n,n,nz,wTi,wTj,Tx,Tz,Ap,Ai,Ax,Az,status)
  197 +call umfpack_zi_triplet_to_col(n,n,nz,Ti,Tj,Tx,Tz,Ap,Ai,Ax,Az,status)
222 198 ! if status is not zero a problem has occured
223 199 if (status /= 0) then
224 200 write(*,*) "Problem during umfpack_zl_triplet_to_col : ",trim(umfpack_return_code(status))
... ... @@ -261,9 +237,7 @@
261 237  
262 238  
263 239  
264   -allocate(bx(n),bz(n),xx(n),xz(n))
265   -bx=dble(B)
266   -bz=aimag(B)
  240 +allocate(xx(n),xz(n))
267 241 sys=0
268 242 ! sys may be used to define type of solving -> see umfpack.h
269 243  
270 244  
... ... @@ -280,10 +254,8 @@
280 254  
281 255 x=cmplx(xx,xz,dp_kind)
282 256  
283   -deallocate(bx,bz,xx,xz)
  257 +deallocate(xx,xz)
284 258 deallocate(Ax,Az)
285   -deallocate(Tx,Tz)
286   -deallocate(wTi,wTj)
287 259 end subroutine
288 260  
289 261  
... ... @@ -301,7 +273,6 @@
301 273 integer(kind=dp_kind), intent(out) :: status
302 274 real(kind=dp_kind), dimension(2), optional, intent(out) :: det
303 275  
304   -integer(kind=dp_kind),dimension(:),allocatable :: wTi,wTj
305 276 real(kind=dp_kind),dimension(:),allocatable :: A
306 277 integer(kind=dp_kind),dimension(:),allocatable :: Ap,Ai
307 278 !integer(kind=dp_kind) :: symbolic,numeric
308 279  
... ... @@ -311,15 +282,11 @@
311 282 integer(kind=dp_kind) :: sys
312 283  
313 284 status=0
314   -! we use a working copy of Ti and Tj to perform 1-based to 0-based translation
315   -allocate(wTi(nz),wTj(nz))
316   -wTi=Ti-1
317   -wTj=Tj-1
318 285 allocate(A(nz))
319 286 allocate(Ap(n+1),Ai(nz))
320 287  
321 288 ! perform the triplet to compressed column form -> Ap,Ai,Ax,Az
322   -call umfpack_dl_triplet_to_col(n,n,nz,wTi,wTj,T,Ap,Ai,A,status)
  289 +call umfpack_dl_triplet_to_col(n,n,nz,Ti,Tj,T,Ap,Ai,A,status)
323 290 ! if status is not zero a problem has occured
324 291 if (status /= 0) then
325 292 write(*,*) "Problem during umfpack_dl_triplet_to_col : ",trim(umfpack_return_code(int(status,kind=sp_kind)))
... ... @@ -375,7 +342,6 @@
375 342 call umfpack_dl_free_numeric (numeric)
376 343  
377 344 deallocate(A)
378   -deallocate(wTi,wTj)
379 345 end subroutine
380 346  
381 347  
... ... @@ -393,7 +359,6 @@
393 359 integer(kind=sp_kind), intent(out) :: status
394 360 real(kind=dp_kind), dimension(2), optional, intent(out) :: det
395 361  
396   -integer(kind=sp_kind),dimension(:),allocatable :: wTi,wTj
397 362 real(kind=dp_kind),dimension(:),allocatable :: A
398 363 integer(kind=sp_kind),dimension(:),allocatable :: Ap,Ai
399 364 !integer(kind=dp_kind) :: symbolic,numeric
400 365  
... ... @@ -406,15 +371,11 @@
406 371 integer(kind=sp_kind) :: sys
407 372  
408 373 status=0
409   -! we use a working copy of Ti and Tj to perform 1-based to 0-based translation
410   -allocate(wTi(nz),wTj(nz))
411   -wTi=Ti-1
412   -wTj=Tj-1
413 374 allocate(A(nz))
414 375 allocate(Ap(n+1),Ai(nz))
415 376  
416 377 ! perform the triplet to compressed column form -> Ap,Ai,Ax,Az
417   -call umfpack_di_triplet_to_col(n,n,nz,wTi,wTj,T,Ap,Ai,A,status)
  378 +call umfpack_di_triplet_to_col(n,n,nz,Ti,Tj,T,Ap,Ai,A,status)
418 379 ! if status is not zero a problem has occured
419 380 if (status /= 0) then
420 381 write(*,*) "Problem during umfpack_di_triplet_to_col : ",trim(umfpack_return_code(status))
... ... @@ -468,7 +429,6 @@
468 429 call umfpack_di_free_numeric (numeric)
469 430  
470 431 deallocate(A)
471   -deallocate(wTi,wTj)
472 432 end subroutine
473 433  
474 434  
475 435  
476 436  
... ... @@ -477,16 +437,14 @@
477 437 ! SPARSE DETERMINANT
478 438 !
479 439 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
480   -subroutine fvn_zl_sparse_det(n,nz,T,Ti,Tj,det,status)
  440 +subroutine fvn_zl_sparse_det(n,nz,Tx,Tz,Ti,Tj,det,status)
481 441 implicit none
482 442 integer(kind=dp_kind), intent(in) :: n,nz
483   -complex(kind=dp_kind),dimension(nz),intent(in) :: T
  443 +real(kind=dp_kind),dimension(nz),intent(in) :: Tx,Tz
484 444 integer(kind=dp_kind),dimension(nz),intent(in) :: Ti,Tj
485 445 real(kind=dp_kind), dimension(3), intent(out) :: det
486 446 integer(kind=dp_kind), intent(out) :: status
487 447  
488   -integer(kind=dp_kind),dimension(:),allocatable :: wTi,wTj
489   -real(kind=dp_kind),dimension(:),allocatable :: Tx,Tz
490 448 real(kind=dp_kind),dimension(:),allocatable :: Ax,Az
491 449 integer(kind=dp_kind),dimension(:),allocatable :: Ap,Ai
492 450 integer(kind=dp_kind) :: symbolic,numeric
493 451  
... ... @@ -497,19 +455,11 @@
497 455  
498 456 status=0
499 457  
500   -! we use a working copy of Ti and Tj to perform 1-based to 0-based translation
501   -! Tx and Tz are the real and imaginary parts of T
502   -allocate(wTi(nz),wTj(nz))
503   -allocate(Tx(nz),Tz(nz))
504   -Tx=dble(T)
505   -Tz=aimag(T)
506   -wTi=Ti-1
507   -wTj=Tj-1
508 458 allocate(Ax(nz),Az(nz))
509 459 allocate(Ap(n+1),Ai(nz))
510 460  
511 461 ! perform the triplet to compressed column form -> Ap,Ai,Ax,Az
512   -call umfpack_zl_triplet_to_col(n,n,nz,wTi,wTj,Tx,Tz,Ap,Ai,Ax,Az,status)
  462 +call umfpack_zl_triplet_to_col(n,n,nz,Ti,Tj,Tx,Tz,Ap,Ai,Ax,Az,status)
513 463 ! if status is not zero a problem has occured
514 464 if (status /= 0) then
515 465 write(*,*) "Problem during umfpack_zl_triplet_to_col ",trim(umfpack_return_code(int(status,kind=sp_kind)))
516 466  
517 467  
518 468  
... ... @@ -550,21 +500,17 @@
550 500 call umfpack_zl_free_numeric (numeric)
551 501  
552 502 deallocate(Ax,Az)
553   -deallocate(Tx,Tz)
554   -deallocate(wTi,wTj)
555 503 end subroutine
556 504  
557 505  
558   -subroutine fvn_zi_sparse_det(n,nz,T,Ti,Tj,det,status)
  506 +subroutine fvn_zi_sparse_det(n,nz,Tx,Tz,Ti,Tj,det,status)
559 507 implicit none
560 508 integer(kind=sp_kind), intent(in) :: n,nz
561   -complex(kind=dp_kind),dimension(nz),intent(in) :: T
  509 +real(kind=dp_kind),dimension(nz),intent(in) :: Tx,Tz
562 510 integer(kind=sp_kind),dimension(nz),intent(in) :: Ti,Tj
563 511 integer(kind=sp_kind), intent(out) :: status
564 512 real(kind=dp_kind), dimension(3), intent(out) :: det
565 513  
566   -integer(kind=sp_kind),dimension(:),allocatable :: wTi,wTj
567   -real(kind=dp_kind),dimension(:),allocatable :: Tx,Tz
568 514 real(kind=dp_kind),dimension(:),allocatable :: Ax,Az
569 515 integer(kind=sp_kind),dimension(:),allocatable :: Ap,Ai
570 516 !integer(kind=dp_kind) :: symbolic,numeric
571 517  
... ... @@ -576,19 +522,11 @@
576 522 real(kind=dp_kind),dimension(20) :: control
577 523  
578 524 status=0
579   -! we use a working copy of Ti and Tj to perform 1-based to 0-based translation
580   -! Tx and Tz are the real and imaginary parts of T
581   -allocate(wTi(nz),wTj(nz))
582   -allocate(Tx(nz),Tz(nz))
583   -Tx=dble(T)
584   -Tz=aimag(T)
585   -wTi=Ti-1
586   -wTj=Tj-1
587 525 allocate(Ax(nz),Az(nz))
588 526 allocate(Ap(n+1),Ai(nz))
589 527  
590 528 ! perform the triplet to compressed column form -> Ap,Ai,Ax,Az
591   -call umfpack_zi_triplet_to_col(n,n,nz,wTi,wTj,Tx,Tz,Ap,Ai,Ax,Az,status)
  529 +call umfpack_zi_triplet_to_col(n,n,nz,Ti,Tj,Tx,Tz,Ap,Ai,Ax,Az,status)
592 530 ! if status is not zero a problem has occured
593 531 if (status /= 0) then
594 532 write(*,*) "Problem during umfpack_zl_triplet_to_col : ",trim(umfpack_return_code(status))
... ... @@ -629,8 +567,6 @@
629 567 call umfpack_zi_free_numeric (numeric)
630 568  
631 569 deallocate(Ax,Az)
632   -deallocate(Tx,Tz)
633   -deallocate(wTi,wTj)
634 570 end subroutine
635 571  
636 572 subroutine fvn_dl_sparse_det(n,nz,T,Ti,Tj,det,status)
... ... @@ -641,7 +577,6 @@
641 577 integer(kind=dp_kind), intent(out) :: status
642 578 real(kind=dp_kind), dimension(2), intent(out) :: det
643 579  
644   -integer(kind=dp_kind),dimension(:),allocatable :: wTi,wTj
645 580 real(kind=dp_kind),dimension(:),allocatable :: A
646 581 integer(kind=dp_kind),dimension(:),allocatable :: Ap,Ai
647 582 !integer(kind=dp_kind) :: symbolic,numeric
648 583  
... ... @@ -650,15 +585,11 @@
650 585 real(kind=dp_kind),dimension(20) :: control
651 586  
652 587 status=0
653   -! we use a working copy of Ti and Tj to perform 1-based to 0-based translation
654   -allocate(wTi(nz),wTj(nz))
655   -wTi=Ti-1
656   -wTj=Tj-1
657 588 allocate(A(nz))
658 589 allocate(Ap(n+1),Ai(nz))
659 590  
660 591 ! perform the triplet to compressed column form -> Ap,Ai,Ax,Az
661   -call umfpack_dl_triplet_to_col(n,n,nz,wTi,wTj,T,Ap,Ai,A,status)
  592 +call umfpack_dl_triplet_to_col(n,n,nz,Ti,Tj,T,Ap,Ai,A,status)
662 593 ! if status is not zero a problem has occured
663 594 if (status /= 0) then
664 595 write(*,*) "Problem during umfpack_dl_triplet_to_col : ",trim(umfpack_return_code(int(status,kind=sp_kind)))
... ... @@ -699,7 +630,6 @@
699 630 call umfpack_dl_free_numeric (numeric)
700 631  
701 632 deallocate(A)
702   -deallocate(wTi,wTj)
703 633 end subroutine
704 634  
705 635  
... ... @@ -711,7 +641,6 @@
711 641 integer(kind=sp_kind), intent(out) :: status
712 642 real(kind=dp_kind), dimension(2), intent(out) :: det
713 643  
714   -integer(kind=sp_kind),dimension(:),allocatable :: wTi,wTj
715 644 real(kind=dp_kind),dimension(:),allocatable :: A
716 645 integer(kind=sp_kind),dimension(:),allocatable :: Ap,Ai
717 646 !integer(kind=dp_kind) :: symbolic,numeric
718 647  
... ... @@ -725,15 +654,11 @@
725 654  
726 655  
727 656 status=0
728   -! we use a working copy of Ti and Tj to perform 1-based to 0-based translation
729   -allocate(wTi(nz),wTj(nz))
730   -wTi=Ti-1
731   -wTj=Tj-1
732 657 allocate(A(nz))
733 658 allocate(Ap(n+1),Ai(nz))
734 659  
735 660 ! perform the triplet to compressed column form -> Ap,Ai,Ax,Az
736   -call umfpack_di_triplet_to_col(n,n,nz,wTi,wTj,T,Ap,Ai,A,status)
  661 +call umfpack_di_triplet_to_col(n,n,nz,Ti,Tj,T,Ap,Ai,A,status)
737 662 ! if status is not zero a problem has occured
738 663 if (status /= 0) then
739 664 write(*,*) "Problem during umfpack_di_triplet_to_col : ",trim(umfpack_return_code(status))
... ... @@ -775,7 +700,6 @@
775 700 call umfpack_di_free_numeric (numeric)
776 701  
777 702 deallocate(A)
778   -deallocate(wTi,wTj)
779 703 end subroutine
780 704  
781 705  
fvn_test/test_sparse_di.f90
Diff comments   View file @ 4bd2122
... ... @@ -22,6 +22,11 @@
22 22 do i=1,nz
23 23 As(Ti(i),Tj(i))=A(i)
24 24 end do
  25 +
  26 +! sparse routines must be fed up with 0-based indices
  27 +Ti=Ti-1
  28 +Tj=Tj-1
  29 +
25 30 write(*,*) "Matrix in standard representation :"
26 31 do i=1,5
27 32 write(*,'(5f8.4)') As(i,:)
fvn_test/test_sparse_dl.f90
Diff comments   View file @ 4bd2122
... ... @@ -22,6 +22,11 @@
22 22 do i=1,nz
23 23 As(Ti(i),Tj(i))=A(i)
24 24 end do
  25 +
  26 +! sparse routines must be fed up with 0-based indices
  27 +Ti=Ti-1
  28 +Tj=Tj-1
  29 +
25 30 write(*,*) "Matrix in standard representation :"
26 31 do i=1,5
27 32 write(*,'(5f8.4)') As(i,:)
fvn_test/test_sparse_zi.f90
Diff comments   View file @ 4bd2122
... ... @@ -4,9 +4,11 @@
4 4 integer(kind=sp_kind), parameter :: nz=12
5 5 integer(kind=sp_kind), parameter :: n=5
6 6 complex(kind=dp_kind),dimension(nz) :: A
  7 +real(kind=dp_kind),dimension(nz) :: Ax,Az
7 8 complex(kind=dp_kind),dimension(n,n) :: As
8 9 integer(kind=sp_kind),dimension(nz) :: Ti,Tj
9 10 complex(kind=dp_kind),dimension(n) :: B,x
  11 +real(kind=dp_kind),dimension(n) :: Bx,Bz
10 12 integer(kind=sp_kind) :: status,i
11 13 real(kind=dp_kind),dimension(3) :: det
12 14 character(len=80) :: fmcmplx
13 15  
... ... @@ -19,12 +21,21 @@
19 21 Ti = (/ 1,2,1,3,5,2,3,4,5,3,2,5 /)
20 22 Tj = (/ 1,1,2,2,2,3,3,3,3,4,5,5 /)
21 23  
  24 +Ax=real(A)
  25 +Az=aimag(A)
  26 +Bx=real(B)
  27 +Bz=aimag(B)
  28 +
22 29 ! Reconstruction of the matrix in standard form
23 30 As=0.
24 31 do i=1,nz
25 32 As(Ti(i),Tj(i))=A(i)
26 33 end do
27 34  
  35 +! sparse routines must be fed up with 0-based indices
  36 +Ti=Ti-1
  37 +Tj=Tj-1
  38 +
28 39 write(*,*) "Matrix in standard representation :"
29 40 do i=1,5
30 41 write(*,fmcmplx) As(i,:)
31 42  
... ... @@ -37,13 +48,13 @@
37 48  
38 49 ! can use either specific interface, fvn_zi_sparse_det
39 50 ! either generic one fvn_sparse_det
40   -call fvn_zi_sparse_det(n,nz,A,Ti,Tj,det,status)
  51 +call fvn_zi_sparse_det(n,nz,Ax,Az,Ti,Tj,det,status)
41 52 write(*,*)
42 53 write(*,*) "Sparse Det = ",cmplx(det(1),det(2),kind=dp_kind)*10**det(3)
43 54 ! can use either specific interface fvn_zi_sparse_solve
44 55 ! either generic one fvn_sparse_solve
45 56 ! parameter det is optional
46   -call fvn_zi_sparse_solve(n,nz,A,Ti,Tj,B,x,status,det)
  57 +call fvn_zi_sparse_solve(n,nz,Ax,Az,Ti,Tj,Bx,Bz,x,status,det)
47 58 write(*,*)
48 59 write(*,*) "Sparse Det as solve option= ",cmplx(det(1),det(2),kind=dp_kind)*10**det(3)
49 60 write(*,*)
fvn_test/test_sparse_zl.f90
Diff comments   View file @ 4bd2122
... ... @@ -5,9 +5,11 @@
5 5 integer(kind=dp_kind), parameter :: nz=12
6 6 integer(kind=dp_kind), parameter :: n=5
7 7 complex(kind=dp_kind),dimension(nz) :: A
  8 +real(kind=dp_kind),dimension(nz) :: Ax,Az
8 9 complex(kind=dp_kind),dimension(n,n) :: As
9 10 integer(kind=dp_kind),dimension(nz) :: Ti,Tj
10 11 complex(kind=dp_kind),dimension(n) :: B,x
  12 +real(kind=dp_kind),dimension(n) :: Bx,Bz
11 13 integer(kind=dp_kind) :: status,i
12 14 real(kind=dp_kind),dimension(3) :: det
13 15 character(len=80) :: fmcmplx
14 16  
... ... @@ -20,12 +22,21 @@
20 22 Ti = (/ 1,2,1,3,5,2,3,4,5,3,2,5 /)
21 23 Tj = (/ 1,1,2,2,2,3,3,3,3,4,5,5 /)
22 24  
  25 +Ax=real(A)
  26 +Az=aimag(A)
  27 +Bx=real(B)
  28 +Bz=aimag(B)
  29 +
23 30 ! Reconstruction of the matrix in standard form
24 31 As=0.
25 32 do i=1,nz
26 33 As(Ti(i),Tj(i))=A(i)
27 34 end do
28 35  
  36 +! sparse routines must be fed up with 0-based indices
  37 +Ti=Ti-1
  38 +Tj=Tj-1
  39 +
29 40 write(*,*) "Matrix in standard representation :"
30 41 do i=1,5
31 42 write(*,fmcmplx) As(i,:)
32 43  
... ... @@ -37,13 +48,13 @@
37 48 write(*,fmcmplx) B
38 49 ! can use either specific interface, fvn_zl_sparse_det
39 50 ! either generic one fvn_sparse_det
40   -call fvn_zl_sparse_det(n,nz,A,Ti,Tj,det,status)
  51 +call fvn_zl_sparse_det(n,nz,Ax,Az,Ti,Tj,det,status)
41 52 write(*,*)
42 53 write(*,*) "Sparse Det = ",cmplx(det(1),det(2),kind=dp_kind)*10**det(3)
43 54 ! can use either specific interface fvn_zl_sparse_solve
44 55 ! either generic one fvn_sparse_solve
45 56 ! parameter det is optional
46   -call fvn_zl_sparse_solve(n,nz,A,Ti,Tj,B,x,status,det)
  57 +call fvn_zl_sparse_solve(n,nz,Ax,Az,Ti,Tj,Bx,Bz,x,status,det)
47 58 write(*,*)
48 59 write(*,*) "Sparse Det as solve option= ",cmplx(det(1),det(2),kind=dp_kind)*10**det(3)
49 60 write(*,*)