OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
zsol_aux.F
Go to the documentation of this file.
1C
2C This file is part of MUMPS 5.5.1, released
3C on Tue Jul 12 13:17:24 UTC 2022
4C
5C
6C Copyright 1991-2022 CERFACS, CNRS, ENS Lyon, INP Toulouse, Inria,
7C Mumps Technologies, University of Bordeaux.
8C
9C This version of MUMPS is provided to you free of charge. It is
10C released under the CeCILL-C license
11C (see doc/CeCILL-C_V1-en.txt, doc/CeCILL-C_V1-fr.txt, and
12C https://cecill.info/licences/Licence_CeCILL-C_V1-en.html)
13C
14 SUBROUTINE zmumps_freetopso( N, KEEP28, IWCB, LIWW,
15 & W, LWC,
16 & POSWCB,IWPOSCB,PTRICB,PTRACB)
17 IMPLICIT NONE
18 INTEGER(8), INTENT(IN) :: LWC
19 INTEGER(8), INTENT(INOUT) :: POSWCB
20 INTEGER N,LIWW,IWPOSCB, KEEP28
21 INTEGER IWCB(LIWW),PTRICB(KEEP28)
22 INTEGER(8) :: PTRACB(KEEP28)
23 COMPLEX(kind=8) W(LWC)
24 INTEGER SIZFI, SIZFR
25 IF ( iwposcb .eq. liww ) RETURN
26 DO WHILE ( iwcb( iwposcb + 2 ) .eq. 0 )
27 sizfr = iwcb( iwposcb + 1 )
28 sizfi = 2
29 iwposcb = iwposcb + sizfi
30 poswcb = poswcb + sizfr
31 IF ( iwposcb .eq. liww ) RETURN
32 END DO
33 RETURN
34 END SUBROUTINE zmumps_freetopso
35 SUBROUTINE zmumps_compso(N,KEEP28,IWCB,LIWW,W,LWC,
36 & POSWCB,IWPOSCB,PTRICB,PTRACB)
37 IMPLICIT NONE
38 INTEGER(8), INTENT(IN) :: LWC
39 INTEGER(8), INTENT(INOUT) :: POSWCB
40 INTEGER N,LIWW,IWPOSCB,KEEP28
41 INTEGER IWCB(LIWW),PTRICB(KEEP28)
42 INTEGER(8) :: PTRACB(KEEP28)
43 COMPLEX(kind=8) W(LWC)
44 INTEGER IPTIW,SIZFI,LONGI
45 INTEGER(8) :: IPTA, LONGR, SIZFR, I8
46 INTEGER :: I
47 iptiw = iwposcb
48 ipta = poswcb
49 longi = 0
50 longr = 0_8
51 IF ( iptiw .EQ. liww ) RETURN
5210 CONTINUE
53 IF (iwcb(iptiw+2).EQ.0) THEN
54 sizfr = int(iwcb(iptiw+1),8)
55 sizfi = 2
56 IF (longi.NE.0) THEN
57 DO 20 i=0,longi-1
58 iwcb(iptiw + sizfi - i) = iwcb(iptiw - i)
59 20 CONTINUE
60 DO 30 i8=0,longr-1
61 w(ipta + sizfr - i8) = w(ipta - i8)
62 30 CONTINUE
63 ENDIF
64 DO 40 i=1,keep28
65 IF ((ptricb(i).LE.(iptiw+1)).AND.
66 & (ptricb(i).GT.iwposcb) ) THEN
67 ptricb(i) = ptricb(i) + sizfi
68 ptracb(i) = ptracb(i) + sizfr
69 ENDIF
7040 CONTINUE
71 iwposcb = iwposcb + sizfi
72 iptiw = iptiw + sizfi
73 poswcb = poswcb + sizfr
74 ipta = ipta + sizfr
75 ELSE
76 sizfr = int(iwcb(iptiw+1),8)
77 sizfi = 2
78 iptiw = iptiw + sizfi
79 longi = longi + sizfi
80 ipta = ipta + sizfr
81 longr = longr + sizfr
82 ENDIF
83 IF (iptiw.NE.liww) GOTO 10
84 RETURN
85 END SUBROUTINE zmumps_compso
86 SUBROUTINE zmumps_sol_x(A, NZ8, N, IRN, ICN, Z, KEEP,KEEP8,
87 & EFF_SIZE_SCHUR, SYM_PERM )
88 INTEGER, INTENT(IN) :: N, KEEP(500)
89 INTEGER(8), INTENT(IN) :: NZ8
90 INTEGER(8), INTENT(IN) :: KEEP8(150)
91 INTEGER, INTENT(IN) :: IRN(NZ8), ICN(NZ8)
92 COMPLEX(kind=8), INTENT(IN) :: A(NZ8)
93 DOUBLE PRECISION, INTENT(OUT) :: Z(N)
94 INTEGER, INTENT(IN) :: EFF_SIZE_SCHUR, SYM_PERM(N)
95 INTEGER :: I, J
96 LOGICAL :: SKIP_COLinSchur
97 DOUBLE PRECISION, PARAMETER :: ZERO = 0.0d0
98 INTEGER(8) :: K
99 INTRINSIC abs
100 DO 10 i = 1, n
101 z(i) = zero
102 10 CONTINUE
103 skip_colinschur = (eff_size_schur.GT.0)
104 IF (keep(264).EQ.0) THEN
105 IF (keep(50) .EQ.0) THEN
106 DO k = 1_8, nz8
107 i = irn(k)
108 j = icn(k)
109 IF ((i .LT. 1) .OR. (i .GT. n)) cycle
110 IF ((j .LT. 1) .OR. (j .GT. n)) cycle
111 IF ( skip_colinschur.AND.
112 & (sym_perm(j).GT.n-eff_size_schur)) cycle
113 IF ( skip_colinschur.AND.
114 & (sym_perm(i).GT.n-eff_size_schur)) cycle
115 z(i) = z(i) + abs(a(k))
116 ENDDO
117 ELSE
118 DO k = 1_8, nz8
119 i = irn(k)
120 j = icn(k)
121 IF ((i .LT. 1) .OR. (i .GT. n)) cycle
122 IF ((j .LT. 1) .OR. (j .GT. n)) cycle
123 IF ( skip_colinschur.AND.
124 & ( (sym_perm(i).GT.n-eff_size_schur)
125 & .OR.
126 & (sym_perm(j).GT.n-eff_size_schur)
127 & )
128 & ) cycle
129 z(i) = z(i) + abs(a(k))
130 IF (j.NE.i) THEN
131 z(j) = z(j) + abs(a(k))
132 ENDIF
133 ENDDO
134 ENDIF
135 ELSE
136 IF (keep(50) .EQ.0) THEN
137 IF (skip_colinschur) THEN
138 DO k = 1_8, nz8
139 j = icn(k)
140 IF ( sym_perm(j).GT.n-eff_size_schur ) cycle
141 i = irn(k)
142 IF ( sym_perm(i).GT.n-eff_size_schur ) cycle
143 z(i) = z(i) + abs(a(k))
144 ENDDO
145 ELSE
146 DO k = 1_8, nz8
147 i = irn(k)
148 j = icn(k)
149 z(i) = z(i) + abs(a(k))
150 ENDDO
151 ENDIF
152 ELSE
153 DO k = 1_8, nz8
154 i = irn(k)
155 j = icn(k)
156 IF ( skip_colinschur.AND.
157 & ( (sym_perm(i).GT.n-eff_size_schur)
158 & .OR.
159 & (sym_perm(j).GT.n-eff_size_schur)
160 & )
161 & ) cycle
162 z(i) = z(i) + abs(a(k))
163 IF (j.NE.i) THEN
164 z(j) = z(j) + abs(a(k))
165 ENDIF
166 ENDDO
167 ENDIF
168 ENDIF
169 RETURN
170 END SUBROUTINE zmumps_sol_x
171 SUBROUTINE zmumps_scal_x(A, NZ8, N, IRN, ICN, Z,
172 & KEEP, KEEP8, COLSCA,
173 & EFF_SIZE_SCHUR, SYM_PERM )
174 INTEGER, INTENT(IN) :: N, KEEP(500)
175 INTEGER(8), INTENT(IN) :: NZ8
176 INTEGER(8), INTENT(IN) :: KEEP8(150)
177 INTEGER, INTENT(IN) :: IRN(NZ8), ICN(NZ8)
178 COMPLEX(kind=8), INTENT(IN) :: A(NZ8)
179 DOUBLE PRECISION, INTENT(IN) :: COLSCA(N)
180 DOUBLE PRECISION, INTENT(OUT) :: Z(N)
181 INTEGER, INTENT(IN) :: EFF_SIZE_SCHUR, SYM_PERM(N)
182 DOUBLE PRECISION, PARAMETER :: ZERO = 0.0d0
183 INTEGER :: I, J
184 INTEGER(8) :: K
185 LOGICAL :: SKIP_COLinSchur
186 DO 10 i = 1, n
187 z(i) = zero
188 10 CONTINUE
189 skip_colinschur = (eff_size_schur.GT.0)
190 IF (keep(50) .EQ.0) THEN
191 DO k = 1_8, nz8
192 i = irn(k)
193 j = icn(k)
194 IF ((i .LT. 1) .OR. (i .GT. n)) cycle
195 IF ((j .LT. 1) .OR. (j .GT. n)) cycle
196 IF ( skip_colinschur.AND.
197 & (sym_perm(j).GT.n-eff_size_schur)) cycle
198 IF ( skip_colinschur.AND.
199 & (sym_perm(i).GT.n-eff_size_schur)) cycle
200 z(i) = z(i) + abs(a(k)*colsca(j))
201 ENDDO
202 ELSE
203 DO k = 1, nz8
204 i = irn(k)
205 j = icn(k)
206 IF ((i .LT. 1) .OR. (i .GT. n)) cycle
207 IF ((j .LT. 1) .OR. (j .GT. n)) cycle
208 IF ( skip_colinschur.AND.
209 & ( (sym_perm(i).GT.n-eff_size_schur)
210 & .OR.
211 & (sym_perm(j).GT.n-eff_size_schur)
212 & )
213 & ) cycle
214 z(i) = z(i) + abs(a(k)*colsca(j))
215 IF (j.NE.i) THEN
216 z(j) = z(j) + abs(a(k)*colsca(i))
217 ENDIF
218 ENDDO
219 ENDIF
220 RETURN
221 END SUBROUTINE zmumps_scal_x
222 SUBROUTINE zmumps_sol_y(A, NZ8, N, IRN, ICN, RHS, X, R, W,
223 & KEEP,KEEP8)
224 IMPLICIT NONE
225 INTEGER, INTENT(IN) :: N, KEEP(500)
226 INTEGER(8), INTENT(IN) :: NZ8
227 INTEGER(8), INTENT(IN) :: KEEP8(150)
228 INTEGER, INTENT(IN) :: IRN(NZ8), ICN(NZ8)
229 COMPLEX(kind=8), INTENT(IN) :: A(NZ8), RHS(N), X(N)
230 DOUBLE PRECISION, INTENT(OUT) :: W(N)
231 COMPLEX(kind=8), INTENT(OUT) :: R(N)
232 INTEGER I, J
233 INTEGER(8) :: K8
234 DOUBLE PRECISION, PARAMETER :: ZERO = 0.0d0
235 COMPLEX(kind=8) D
236 DO i = 1, n
237 r(i) = rhs(i)
238 w(i) = zero
239 ENDDO
240 IF (keep(264).EQ.0) THEN
241 IF (keep(50) .EQ.0) THEN
242 DO k8 = 1_8, nz8
243 i = irn(k8)
244 j = icn(k8)
245 IF ((i .GT. n) .OR. (j .GT. n) .OR. (i .LT. 1) .OR.
246 & (j .LT. 1)) cycle
247 d = a(k8) * x(j)
248 r(i) = r(i) - d
249 w(i) = w(i) + abs(d)
250 ENDDO
251 ELSE
252 DO k8 = 1_8, nz8
253 i = irn(k8)
254 j = icn(k8)
255 IF ((i .GT. n) .OR. (j .GT. n) .OR. (i .LT. 1) .OR.
256 & (j .LT. 1)) cycle
257 d = a(k8) * x(j)
258 r(i) = r(i) - d
259 w(i) = w(i) + abs(d)
260 IF (i.NE.j) THEN
261 d = a(k8) * x(i)
262 r(j) = r(j) - d
263 w(j) = w(j) + abs(d)
264 ENDIF
265 ENDDO
266 ENDIF
267 ELSE
268 IF (keep(50) .EQ.0) THEN
269 DO k8 = 1_8, nz8
270 i = irn(k8)
271 j = icn(k8)
272 d = a(k8) * x(j)
273 r(i) = r(i) - d
274 w(i) = w(i) + abs(d)
275 ENDDO
276 ELSE
277 DO k8 = 1_8, nz8
278 i = irn(k8)
279 j = icn(k8)
280 d = a(k8) * x(j)
281 r(i) = r(i) - d
282 w(i) = w(i) + abs(d)
283 IF (i.NE.j) THEN
284 d = a(k8) * x(i)
285 r(j) = r(j) - d
286 w(j) = w(j) + abs(d)
287 ENDIF
288 ENDDO
289 ENDIF
290 ENDIF
291 RETURN
292 END SUBROUTINE zmumps_sol_y
293 SUBROUTINE zmumps_sol_mulr(N, R, W)
294 INTEGER, intent(in) :: N
295 DOUBLE PRECISION, intent(in) :: W(N)
296 COMPLEX(kind=8), intent(inout) :: R(N)
297 INTEGER I
298 DO 10 i = 1, n
299 r(i) = r(i) * w(i)
300 10 CONTINUE
301 RETURN
302 END SUBROUTINE zmumps_sol_mulr
303 SUBROUTINE zmumps_sol_b(N, KASE, X, EST, W, IW, GRAIN)
304 INTEGER, intent(in) :: N
305 INTEGER, intent(inout) :: KASE
306 INTEGER IW(N)
307 COMPLEX(kind=8) W(N), X(N)
308 DOUBLE PRECISION, intent(inout) :: EST
309 INTEGER, intent(in) :: GRAIN
310 INTRINSIC abs, nint, real, sign
311 INTEGER ZMUMPS_IXAMAX
312 EXTERNAL zmumps_ixamax
313 INTEGER ITMAX
314 parameter(itmax = 5)
315 INTEGER I, ITER, J, JLAST, JUMP
316 DOUBLE PRECISION ALTSGN
317 DOUBLE PRECISION TEMP
318 SAVE iter, j, jlast, jump
319 COMPLEX(kind=8) ZERO, ONE
320 parameter( zero = (0.0d0,0.0d0) )
321 parameter( one = (1.0d0,0.0d0) )
322 DOUBLE PRECISION, PARAMETER :: RZERO = 0.0d0
323 DOUBLE PRECISION, PARAMETER :: RONE = 1.0d0
324 IF (kase .EQ. 0) THEN
325 DO 10 i = 1, n
326 x(i) = one / dble(n)
327 10 CONTINUE
328 kase = 1
329 jump = 1
330 RETURN
331 ENDIF
332 SELECT CASE (jump)
333 CASE (1)
334 GOTO 20
335 CASE(2)
336 GOTO 40
337 CASE(3)
338 GOTO 70
339 CASE(4)
340 GOTO 120
341 CASE(5)
342 GOTO 160
343 CASE DEFAULT
344 END SELECT
345 20 CONTINUE
346 IF (n .EQ. 1) THEN
347 w(1) = x(1)
348 est = abs(w(1))
349 GOTO 190
350 ENDIF
351 DO 30 i = 1, n
352 x(i) = cmplx( sign(rone,dble(x(i))), kind=kind(x))
353 iw(i) = nint(dble(x(i)))
354 30 CONTINUE
355 kase = 2
356 jump = 2
357 RETURN
358 40 CONTINUE
359 j = zmumps_ixamax(n, x, 1, grain)
360 iter = 2
361 50 CONTINUE
362 DO 60 i = 1, n
363 x(i) = zero
364 60 CONTINUE
365 x(j) = one
366 kase = 1
367 jump = 3
368 RETURN
369 70 CONTINUE
370 DO 80 i = 1, n
371 w(i) = x(i)
372 80 CONTINUE
373 DO 90 i = 1, n
374 IF (nint(sign(rone, dble(x(i)))) .NE. iw(i)) GOTO 100
375 90 CONTINUE
376 GOTO 130
377 100 CONTINUE
378 DO 110 i = 1, n
379 x(i) = cmplx( sign(rone, dble(x(i))), kind=kind(x) )
380 iw(i) = nint(dble(x(i)))
381 110 CONTINUE
382 kase = 2
383 jump = 4
384 RETURN
385 120 CONTINUE
386 jlast = j
387 j = zmumps_ixamax(n, x, 1, grain)
388 IF ((abs(x(jlast)) .NE. abs(x(j))) .AND. (iter .LT. itmax)) THEN
389 iter = iter + 1
390 GOTO 50
391 ENDIF
392 130 CONTINUE
393 est = rzero
394 DO 140 i = 1, n
395 est = est + abs(w(i))
396 140 CONTINUE
397 altsgn = rone
398 DO 150 i = 1, n
399 x(i) = cmplx(altsgn * (rone + dble(i - 1) / dble(n - 1)),
400 & kind=kind(x))
401 altsgn = -altsgn
402 150 CONTINUE
403 kase = 1
404 jump = 5
405 RETURN
406 160 CONTINUE
407 temp = rzero
408 DO 170 i = 1, n
409 temp = temp + abs(x(i))
410 170 CONTINUE
411 temp = 2.0d0 * temp / dble(3 * n)
412 IF (temp .GT. est) THEN
413 DO 180 i = 1, n
414 w(i) = x(i)
415 180 CONTINUE
416 est = temp
417 ENDIF
418 190 kase = 0
419 RETURN
420 END SUBROUTINE zmumps_sol_b
421 SUBROUTINE zmumps_qd2( MTYPE, N, NZ8, ASPK, IRN, ICN,
422 & LHS, WRHS, W, RHS, KEEP,KEEP8)
423 IMPLICIT NONE
424 INTEGER MTYPE, N
425 INTEGER(8), INTENT(IN) :: NZ8
426 INTEGER, INTENT(IN) :: IRN( NZ8 ), ICN( NZ8 )
427 INTEGER KEEP(500)
428 INTEGER(8) KEEP8(150)
429 COMPLEX(kind=8), INTENT(IN) :: ASPK( NZ8 )
430 COMPLEX(kind=8), INTENT(IN) :: LHS( N ), WRHS( N )
431 COMPLEX(kind=8), INTENT(OUT):: RHS( N )
432 DOUBLE PRECISION, INTENT(OUT):: W( N )
433 INTEGER I, J
434 INTEGER(8) :: K8
435 DOUBLE PRECISION, PARAMETER :: DZERO = 0.0d0
436 DO i = 1, n
437 w(i) = dzero
438 rhs(i) = wrhs(i)
439 ENDDO
440 IF ( keep(50) .EQ. 0 ) THEN
441 IF (mtype .EQ. 1) THEN
442 IF (keep(264).EQ.0) THEN
443 DO k8 = 1_8, nz8
444 i = irn(k8)
445 j = icn(k8)
446 IF ((i .LE. 0) .OR. (i .GT. n) .OR. (j .LE. 0) .OR.
447 & (j .GT. n)) cycle
448 rhs(i) = rhs(i) - aspk(k8) * lhs(j)
449 w(i) = w(i) + abs(aspk(k8))
450 ENDDO
451 ELSE
452 DO k8 = 1_8, nz8
453 i = irn(k8)
454 j = icn(k8)
455 rhs(i) = rhs(i) - aspk(k8) * lhs(j)
456 w(i) = w(i) + abs(aspk(k8))
457 ENDDO
458 ENDIF
459 ELSE
460 IF (keep(264).EQ.0) THEN
461 DO k8 = 1_8, nz8
462 i = irn(k8)
463 j = icn(k8)
464 IF ((i .LE. 0) .OR. (i .GT. n) .OR. (j .LE. 0) .OR.
465 & (j .GT. n)) cycle
466 rhs(j) = rhs(j) - aspk(k8) * lhs(i)
467 w(j) = w(j) + abs(aspk(k8))
468 ENDDO
469 ELSE
470 DO k8 = 1_8, nz8
471 i = irn(k8)
472 j = icn(k8)
473 rhs(j) = rhs(j) - aspk(k8) * lhs(i)
474 w(j) = w(j) + abs(aspk(k8))
475 ENDDO
476 ENDIF
477 ENDIF
478 ELSE
479 IF (keep(264).EQ.0) THEN
480 DO k8 = 1_8, nz8
481 i = irn(k8)
482 j = icn(k8)
483 IF ((i .LE. 0) .OR. (i .GT. n) .OR. (j .LE. 0) .OR.
484 & (j .GT. n)) cycle
485 rhs(i) = rhs(i) - aspk(k8) * lhs(j)
486 w(i) = w(i) + abs(aspk(k8))
487 IF (j.NE.i) THEN
488 rhs(j) = rhs(j) - aspk(k8) * lhs(i)
489 w(j) = w(j) + abs(aspk(k8))
490 ENDIF
491 ENDDO
492 ELSE
493 DO k8 = 1_8, nz8
494 i = irn(k8)
495 j = icn(k8)
496 rhs(i) = rhs(i) - aspk(k8) * lhs(j)
497 w(i) = w(i) + abs(aspk(k8))
498 IF (j.NE.i) THEN
499 rhs(j) = rhs(j) - aspk(k8) * lhs(i)
500 w(j) = w(j) + abs(aspk(k8))
501 ENDIF
502 ENDDO
503 ENDIF
504 ENDIF
505 RETURN
506 END SUBROUTINE zmumps_qd2
507 SUBROUTINE zmumps_eltqd2( MTYPE, N,
508 & NELT, ELTPTR, LELTVAR, ELTVAR, NA_ELT8, A_ELT,
509 & LHS, WRHS, W, RHS, KEEP,KEEP8 )
510 IMPLICIT NONE
511 INTEGER MTYPE, N, NELT, LELTVAR
512 INTEGER(8), INTENT(IN) :: NA_ELT8
513 INTEGER ELTPTR(NELT+1), ELTVAR(LELTVAR)
514 INTEGER KEEP(500)
515 INTEGER(8) KEEP8(150)
516 COMPLEX(kind=8) A_ELT(NA_ELT8)
517 COMPLEX(kind=8) LHS( N ), WRHS( N ), RHS( N )
518 DOUBLE PRECISION W(N)
519 CALL ZMUMPS_MV_ELT(N, NELT, ELTPTR, ELTVAR, A_ELT,
520 & lhs, rhs, keep(50), mtype )
521 rhs = wrhs - rhs
522 CALL zmumps_sol_x_elt( mtype, n,
523 & nelt, eltptr, leltvar, eltvar, na_elt8, a_elt,
524 & w, keep,keep8 )
525 RETURN
526 END SUBROUTINE zmumps_eltqd2
527 SUBROUTINE zmumps_sol_x_elt( MTYPE, N,
528 & NELT, ELTPTR, LELTVAR, ELTVAR, NA_ELT8, A_ELT,
529 & W, KEEP,KEEP8 )
530 IMPLICIT NONE
531 INTEGER MTYPE, N, NELT, LELTVAR
532 INTEGER(8), INTENT(IN) :: NA_ELT8
533 INTEGER ELTPTR(NELT+1), ELTVAR(LELTVAR)
534 INTEGER KEEP(500)
535 INTEGER(8) KEEP8(150)
536 COMPLEX(kind=8) A_ELT(NA_ELT8)
537 DOUBLE PRECISION TEMP
538 DOUBLE PRECISION W(N)
539 INTEGER I, J, IEL, SIZEI, IELPTR
540 INTEGER(8) :: K8
541 DOUBLE PRECISION DZERO
542 parameter(dzero = 0.0d0)
543 w = dzero
544 k8 = 1_8
545 DO iel = 1, nelt
546 sizei = eltptr( iel + 1 ) - eltptr( iel )
547 ielptr = eltptr( iel ) - 1
548 IF ( keep(50).EQ.0 ) THEN
549 IF (mtype.EQ.1) THEN
550 DO j = 1, sizei
551 DO i = 1, sizei
552 w( eltvar( ielptr + i) ) =
553 & w( eltvar( ielptr + i) )
554 & + abs(a_elt( k8 ))
555 k8 = k8 + 1_8
556 END DO
557 END DO
558 ELSE
559 DO j = 1, sizei
560 temp = w( eltvar( ielptr + j ) )
561 DO i = 1, sizei
562 temp = temp + abs( a_elt(k8))
563 k8 = k8 + 1_8
564 END DO
565 w(eltvar( ielptr + j )) =
566 & w(eltvar( ielptr + j )) + temp
567 END DO
568 ENDIF
569 ELSE
570 DO j = 1, sizei
571 w(eltvar( ielptr + j )) =
572 & w(eltvar( ielptr + j )) + abs(a_elt( k8 ))
573 k8 = k8 + 1_8
574 DO i = j+1, sizei
575 w(eltvar( ielptr + j )) =
576 & w(eltvar( ielptr + j )) + abs(a_elt( k8 ))
577 w(eltvar( ielptr + i ) ) =
578 & w(eltvar( ielptr + i )) + abs(a_elt( k8 ))
579 k8 = k8 + 1_8
580 END DO
581 ENDDO
582 ENDIF
583 ENDDO
584 RETURN
585 END SUBROUTINE zmumps_sol_x_elt
586 SUBROUTINE zmumps_sol_scalx_elt(MTYPE, N,
587 & NELT, ELTPTR, LELTVAR, ELTVAR, NA_ELT8, A_ELT,
588 & W, KEEP,KEEP8, COLSCA )
589 IMPLICIT NONE
590 INTEGER MTYPE, N, NELT, LELTVAR
591 INTEGER(8), INTENT(IN) :: NA_ELT8
592 INTEGER ELTPTR(NELT+1), ELTVAR(LELTVAR)
593 INTEGER KEEP(500)
594 INTEGER(8) KEEP8(150)
595 DOUBLE PRECISION COLSCA(N)
596 COMPLEX(kind=8) A_ELT(NA_ELT8)
597 DOUBLE PRECISION W(N)
598 DOUBLE PRECISION TEMP, TEMP2
599 INTEGER I, J, IEL, SIZEI, IELPTR
600 INTEGER(8) :: K8
601 DOUBLE PRECISION DZERO
602 PARAMETER(DZERO = 0.0d0)
603 w = dzero
604 k8 = 1_8
605 DO iel = 1, nelt
606 sizei = eltptr( iel + 1 ) - eltptr( iel )
607 ielptr = eltptr( iel ) - 1
608 IF ( keep(50).EQ.0 ) THEN
609 IF (mtype.EQ.1) THEN
610 DO j = 1, sizei
611 temp2 = abs(colsca(eltvar( ielptr + j) ))
612 DO i = 1, sizei
613 w( eltvar( ielptr + i) ) =
614 & w( eltvar( ielptr + i) )
615 & + abs(a_elt( k8 )) * temp2
616 k8 = k8 + 1_8
617 END DO
618 END DO
619 ELSE
620 DO j = 1, sizei
621 temp = w( eltvar( ielptr + j ) )
622 temp2= abs(colsca(eltvar( ielptr + j) ))
623 DO i = 1, sizei
624 temp = temp + abs(a_elt( k8 )) * temp2
625 k8 = k8 + 1_8
626 END DO
627 w(eltvar( ielptr + j )) =
628 & w(eltvar( ielptr + j )) + temp
629 END DO
630 ENDIF
631 ELSE
632 DO j = 1, sizei
633 w(eltvar( ielptr + j )) =
634 & w(eltvar( ielptr + j )) +
635 & abs( a_elt( k8 )*colsca(eltvar( ielptr + j)) )
636 k8 = k8 + 1_8
637 DO i = j+1, sizei
638 w(eltvar( ielptr + j )) =
639 & w(eltvar( ielptr + j )) +
640 & abs(a_elt( k8 )*colsca(eltvar( ielptr + j)))
641 w(eltvar( ielptr + i ) ) =
642 & w(eltvar( ielptr + i )) +
643 & abs(a_elt( k8 )*colsca(eltvar( ielptr + i)))
644 k8 = k8 + 1_8
645 END DO
646 ENDDO
647 ENDIF
648 ENDDO
649 RETURN
650 END SUBROUTINE zmumps_sol_scalx_elt
651 SUBROUTINE zmumps_eltyd( MTYPE, N, NELT, ELTPTR,
652 & LELTVAR, ELTVAR, NA_ELT8, A_ELT,
653 & SAVERHS, X, Y, W, K50 )
654 IMPLICIT NONE
655 INTEGER N, NELT, K50, MTYPE, LELTVAR
656 INTEGER(8) :: NA_ELT8
657 INTEGER ELTPTR( NELT + 1 ), ELTVAR( LELTVAR )
658 COMPLEX(kind=8) A_ELT( NA_ELT8 ), X( N ), Y( N ),
659 & SAVERHS(N)
660 DOUBLE PRECISION W(N)
661 INTEGER IEL, I , J, K, SIZEI, IELPTR
662 DOUBLE PRECISION ZERO
663 COMPLEX(kind=8) TEMP
664 DOUBLE PRECISION TEMP2
665 parameter( zero = 0.0d0 )
666 y = saverhs
667 w = zero
668 k = 1
669 DO iel = 1, nelt
670 sizei = eltptr( iel + 1 ) - eltptr( iel )
671 ielptr = eltptr( iel ) - 1
672 IF ( k50 .eq. 0 ) THEN
673 IF ( mtype .eq. 1 ) THEN
674 DO j = 1, sizei
675 temp = x( eltvar( ielptr + j ) )
676 DO i = 1, sizei
677 y( eltvar( ielptr + i ) ) =
678 & y( eltvar( ielptr + i ) ) -
679 & a_elt( k ) * temp
680 w( eltvar( ielptr + i ) ) =
681 & w( eltvar( ielptr + i ) ) +
682 & abs( a_elt( k ) * temp )
683 k = k + 1
684 END DO
685 END DO
686 ELSE
687 DO j = 1, sizei
688 temp = y( eltvar( ielptr + j ) )
689 temp2 = w( eltvar( ielptr + j ) )
690 DO i = 1, sizei
691 temp = temp -
692 & a_elt( k ) * x( eltvar( ielptr + i ) )
693 temp2 = temp2 + abs(
694 & a_elt( k ) * x( eltvar( ielptr + i ) ) )
695 k = k + 1
696 END DO
697 y( eltvar( ielptr + j ) ) = temp
698 w( eltvar( ielptr + j ) ) = temp2
699 END DO
700 END IF
701 ELSE
702 DO j = 1, sizei
703 y( eltvar( ielptr + j ) ) =
704 & y( eltvar( ielptr + j ) ) -
705 & a_elt( k ) * x( eltvar( ielptr + j ) )
706 w( eltvar( ielptr + j ) ) =
707 & w( eltvar( ielptr + j ) ) + abs(
708 & a_elt( k ) * x( eltvar( ielptr + j ) ) )
709 k = k + 1
710 DO i = j+1, sizei
711 y( eltvar( ielptr + i ) ) =
712 & y( eltvar( ielptr + i ) ) -
713 & a_elt( k ) * x( eltvar( ielptr + j ) )
714 y( eltvar( ielptr + j ) ) =
715 & y( eltvar( ielptr + j ) ) -
716 & a_elt( k ) * x( eltvar( ielptr + i ) )
717 w( eltvar( ielptr + i ) ) =
718 & w( eltvar( ielptr + i ) ) + abs(
719 & a_elt( k ) * x( eltvar( ielptr + j ) ) )
720 w( eltvar( ielptr + j ) ) =
721 & w( eltvar( ielptr + j ) ) + abs(
722 & a_elt( k ) * x( eltvar( ielptr + i ) ) )
723 k = k + 1
724 END DO
725 END DO
726 END IF
727 END DO
728 RETURN
729 END SUBROUTINE zmumps_eltyd
730 SUBROUTINE zmumps_solve_get_ooc_node(
731 & INODE,PTRFAC,KEEP,A,LA,STEP,
732 & KEEP8,N,MUST_BE_PERMUTED,IERR)
733 USE zmumps_ooc
734 IMPLICIT NONE
735 INTEGER INODE,KEEP(500),N
736 INTEGER(8) KEEP8(150)
737 INTEGER(8) :: LA
738 INTEGER(8) :: PTRFAC(KEEP(28))
739 INTEGER STEP(N)
740 INTEGER IERR
741 COMPLEX(kind=8) A(LA)
742 INTEGER RETURN_VALUE
743 LOGICAL MUST_BE_PERMUTED
744 return_value=zmumps_solve_is_inode_in_mem(inode,ptrfac,
745 & keep(28),a,la,ierr)
746 IF(return_value.EQ.ooc_node_not_in_mem)THEN
747 IF(ierr.LT.0)THEN
748 RETURN
749 ENDIF
750 CALL zmumps_solve_alloc_factor_space(inode,ptrfac,
751 & keep,keep8,a,ierr)
752 IF(ierr.LT.0)THEN
753 RETURN
754 ENDIF
755 CALL zmumps_read_ooc(
756 & a(ptrfac(step(inode))),
757 & inode,ierr
758 & )
759 IF(ierr.LT.0)THEN
760 RETURN
761 ENDIF
762 ELSE
763 IF(ierr.LT.0)THEN
764 RETURN
765 ENDIF
766 ENDIF
767 IF(return_value.NE.ooc_node_permuted)THEN
768 must_be_permuted=.true.
769 CALL zmumps_solve_modify_state_node(inode)
770 ELSE
771 must_be_permuted=.false.
772 ENDIF
773 RETURN
774 END SUBROUTINE zmumps_solve_get_ooc_node
775 SUBROUTINE zmumps_build_mapping_info(id)
776 USE zmumps_struc_def
777 IMPLICIT NONE
778 include 'mpif.h'
779 TYPE(zmumps_struc), TARGET :: id
780 INTEGER, ALLOCATABLE, DIMENSION(:) :: LOCAL_LIST
781 INTEGER :: I,IERR,TMP,NSTEPS,N_LOCAL_LIST
782 INTEGER :: MASTER,TAG_SIZE,TAG_LIST
783 INTEGER :: STATUS(MPI_STATUS_SIZE)
784 LOGICAL :: I_AM_SLAVE
785 parameter(master=0, tag_size=85,tag_list=86)
786 i_am_slave = (id%MYID .NE. master
787 & .OR. ((id%MYID.EQ.master).AND.(id%KEEP(46).EQ.1)))
788 nsteps = id%KEEP(28)
789 ALLOCATE(local_list(nsteps),stat=ierr)
790 IF(ierr.GT.0) THEN
791 WRITE(*,*)'Problem in solve: error allocating LOCAL_LIST'
792 CALL mumps_abort()
793 END IF
794 n_local_list = 0
795 IF(i_am_slave) THEN
796 DO i=1,nsteps
797 IF(id%PTLUST_S(i).NE.0) THEN
798 n_local_list = n_local_list + 1
799 local_list(n_local_list) = i
800 END IF
801 END DO
802 IF(id%MYID.NE.master) THEN
803 CALL mpi_send(n_local_list, 1,
804 & mpi_integer, master, tag_size, id%COMM,ierr)
805 CALL mpi_send(local_list, n_local_list,
806 & mpi_integer, master, tag_list, id%COMM,ierr)
807 DEALLOCATE(local_list)
808 ALLOCATE(id%IPTR_WORKING(1),
809 & id%WORKING(1),
810 & stat=ierr)
811 IF(ierr.GT.0) THEN
812 WRITE(*,*)'Problem in solve: error allocating ',
813 & 'IPTR_WORKING and WORKING'
814 CALL mumps_abort()
815 END IF
816 END IF
817 END IF
818 IF(id%MYID.EQ.master) THEN
819 ALLOCATE(id%IPTR_WORKING(id%NPROCS+1), stat=ierr)
820 IF(ierr.GT.0) THEN
821 WRITE(*,*)'Problem in solve: error allocating IPTR_WORKING'
822 CALL mumps_abort()
823 END IF
824 id%IPTR_WORKING = 0
825 id%IPTR_WORKING(1) = 1
826 id%IPTR_WORKING(master+2) = n_local_list
827 DO i=1, id%NPROCS-1
828 CALL mpi_recv(tmp, 1, mpi_integer, mpi_any_source,
829 & tag_size, id%COMM, status, ierr)
830 id%IPTR_WORKING(status(mpi_source)+2) = tmp
831 END DO
832 DO i=2, id%NPROCS+1
833 id%IPTR_WORKING(i) = id%IPTR_WORKING(i)
834 & + id%IPTR_WORKING(i-1)
835 END DO
836 ALLOCATE(id%WORKING(id%IPTR_WORKING(id%NPROCS+1)-1),stat=ierr)
837 IF(ierr.GT.0) THEN
838 WRITE(*,*)'Problem in solve: error allocating LOCAL_LIST'
839 CALL mumps_abort()
840 END IF
841 tmp = master + 1
842 IF (i_am_slave) THEN
843 id%WORKING(id%IPTR_WORKING(tmp):id%IPTR_WORKING(tmp+1)-1)
844 & = local_list(1:id%IPTR_WORKING(tmp+1)
845 & -id%IPTR_WORKING(tmp))
846 ENDIF
847 DO i=1,id%NPROCS-1
848 CALL mpi_recv(local_list, nsteps, mpi_integer,
849 & mpi_any_source, tag_list, id%COMM, status, ierr)
850 tmp = status(mpi_source)+1
851 id%WORKING(id%IPTR_WORKING(tmp):id%IPTR_WORKING(tmp+1)-1)
852 & = local_list(1:id%IPTR_WORKING(tmp+1)-
853 & id%IPTR_WORKING(tmp))
854 END DO
855 DEALLOCATE(local_list)
856 END IF
857 END SUBROUTINE zmumps_build_mapping_info
858 SUBROUTINE zmumps_sol_omega(N, RHS,
859 & X, Y, R_W, C_W, IW, IFLAG,
860 & OMEGA, NOITER, TESTConv,
861 & LP, ARRET, GRAIN )
862 IMPLICIT NONE
863 INTEGER N, IFLAG
864 INTEGER IW(N,2)
865 COMPLEX(kind=8) RHS(N)
866 COMPLEX(kind=8) X(N), Y(N)
867 DOUBLE PRECISION R_W(N,2)
868 COMPLEX(kind=8) C_W(N)
869 INTEGER LP, NOITER
870 LOGICAL TESTConv
871 DOUBLE PRECISION OMEGA(2)
872 DOUBLE PRECISION ARRET
873 INTEGER, intent(in) :: GRAIN
874 DOUBLE PRECISION, PARAMETER :: CGCE=0.2d0
875 DOUBLE PRECISION, PARAMETER :: CTAU=1.0d3
876 INTEGER I, IMAX
877 DOUBLE PRECISION OM1, OM2, DXMAX
878 DOUBLE PRECISION TAU, DD
879 DOUBLE PRECISION OLDOMG(2)
880 DOUBLE PRECISION, PARAMETER :: ZERO=0.0d0
881 DOUBLE PRECISION, PARAMETER :: ONE=1.0d0
882 INTEGER ZMUMPS_IXAMAX
883 INTRINSIC abs, max
884 SAVE om1, oldomg
885 imax = zmumps_ixamax(n, x, 1, grain)
886 dxmax = abs(x(imax))
887 omega(1) = zero
888 omega(2) = zero
889 DO i = 1, n
890 tau = (r_w(i, 2) * dxmax + abs(rhs(i))) * dble(n) * ctau
891 dd = r_w(i, 1) + abs(rhs(i))
892 IF (dd .GT. tau * epsilon(ctau)) THEN
893 omega(1) = max(omega(1), abs(y(i)) / dd)
894 iw(i, 1) = 1
895 ELSE
896 IF (tau .GT. zero) THEN
897 omega(2) = max(omega(2),
898 & abs(y(i)) / (dd + r_w(i, 2) * dxmax))
899 ENDIF
900 iw(i, 1) = 2
901 ENDIF
902 ENDDO
903 IF (testconv) THEN
904 om2 = omega(1) + omega(2)
905 IF (om2 .LT. arret ) THEN
906 iflag = 1
907 GOTO 70
908 ENDIF
909 IF (noiter .GE. 1) THEN
910 IF (om2 .GT. om1 * cgce) THEN
911 IF (om2 .GT. om1) THEN
912 omega(1) = oldomg(1)
913 omega(2) = oldomg(2)
914 DO i = 1, n
915 x(i) = c_w(i)
916 ENDDO
917 iflag = 2
918 GOTO 70
919 ENDIF
920 iflag = 3
921 GOTO 70
922 ENDIF
923 ENDIF
924 DO i = 1, n
925 c_w(i) = x(i)
926 ENDDO
927 oldomg(1) = omega(1)
928 oldomg(2) = omega(2)
929 om1 = om2
930 ENDIF
931 iflag = 0
932 RETURN
933 70 CONTINUE
934 RETURN
935 END SUBROUTINE zmumps_sol_omega
936 SUBROUTINE zmumps_sol_lcond(N, RHS,
937 & X, Y, D, R_W, C_W, IW, KASE,
938 & OMEGA, ERX, COND,
939 & LP, KEEP,KEEP8 )
940 IMPLICIT NONE
941 INTEGER N, KASE, KEEP(500)
942 INTEGER(8) KEEP8(150)
943 INTEGER IW(N,2)
944 COMPLEX(kind=8) RHS(N)
945 COMPLEX(kind=8) X(N), Y(N)
946 DOUBLE PRECISION D(N)
947 DOUBLE PRECISION R_W(N,2)
948 COMPLEX(kind=8) C_W(N)
949 INTEGER LP
950 DOUBLE PRECISION COND(2),OMEGA(2)
951 LOGICAL LCOND1, LCOND2
952 INTEGER JUMP, I, IMAX
953 DOUBLE PRECISION ERX, DXMAX
954 DOUBLE PRECISION DXIMAX
955 DOUBLE PRECISION, PARAMETER :: ZERO = 0.0d0
956 DOUBLE PRECISION, PARAMETER :: ONE = 1.0d0
957 INTEGER ZMUMPS_IXAMAX
958 INTRINSIC abs, max
959 SAVE lcond1, lcond2, jump, dximax, dxmax
960 IF (kase .EQ. 0) THEN
961 lcond1 = .false.
962 lcond2 = .false.
963 cond(1) = one
964 cond(2) = one
965 erx = zero
966 jump = 1
967 ENDIF
968 SELECT CASE (jump)
969 CASE (1)
970 GOTO 30
971 CASE(2)
972 GOTO 10
973 CASE(3)
974 GOTO 110
975 CASE(4)
976 GOTO 150
977 CASE(5)
978 GOTO 35
979 CASE DEFAULT
980 END SELECT
981 10 CONTINUE
982 30 CONTINUE
983 35 CONTINUE
984 imax = zmumps_ixamax(n, x, 1, keep(361))
985 dxmax = abs(x(imax))
986 DO i = 1, n
987 IF (iw(i, 1) .EQ. 1) THEN
988 r_w(i, 1) = r_w(i, 1) + abs(rhs(i))
989 r_w(i, 2) = zero
990 lcond1 = .true.
991 ELSE
992 r_w(i, 2) = r_w(i, 2) * dxmax + r_w(i, 1)
993 r_w(i, 1) = zero
994 lcond2 = .true.
995 ENDIF
996 ENDDO
997 DO i = 1, n
998 c_w(i) = x(i) * d(i)
999 ENDDO
1000 imax = zmumps_ixamax(n, c_w(1), 1, keep(361))
1001 dximax = abs(c_w(imax))
1002 IF (.NOT.lcond1) GOTO 130
1003 100 CONTINUE
1004 CALL zmumps_sol_b(n, kase, y, cond(1), c_w, iw(1, 2), keep(361))
1005 IF (kase .EQ. 0) GOTO 120
1006 IF (kase .EQ. 1) CALL zmumps_sol_mulr(n, y, d)
1007 IF (kase .EQ. 2) CALL zmumps_sol_mulr(n, y, r_w)
1008 jump = 3
1009 RETURN
1010 110 CONTINUE
1011 IF (kase .EQ. 1) CALL zmumps_sol_mulr(n, y, r_w)
1012 IF (kase .EQ. 2) CALL zmumps_sol_mulr(n, y, d)
1013 GOTO 100
1014 120 CONTINUE
1015 IF (dximax .GT. zero) cond(1) = cond(1) / dximax
1016 erx = omega(1) * cond(1)
1017 130 CONTINUE
1018 IF (.NOT.lcond2) GOTO 170
1019 kase = 0
1020 140 CONTINUE
1021 CALL zmumps_sol_b(n, kase, y, cond(2), c_w, iw(1, 2), keep(361))
1022 IF (kase .EQ. 0) GOTO 160
1023 IF (kase .EQ. 1) CALL zmumps_sol_mulr(n, y, d)
1024 IF (kase .EQ. 2) CALL zmumps_sol_mulr(n, y, r_w(1, 2))
1025 jump = 4
1026 RETURN
1027 150 CONTINUE
1028 IF (kase .EQ. 1) CALL zmumps_sol_mulr(n, y, r_w(1, 2))
1029 IF (kase .EQ. 2) CALL zmumps_sol_mulr(n, y, d)
1030 GOTO 140
1031 160 IF (dximax .GT. zero) THEN
1032 cond(2) = cond(2) / dximax
1033 ENDIF
1034 erx = erx + omega(2) * cond(2)
1035 170 CONTINUE
1036 RETURN
1037 END SUBROUTINE zmumps_sol_lcond
1038 SUBROUTINE zmumps_sol_cpy_fs2rhscomp( JBDEB, JBFIN, NBROWS,
1039 & KEEP, RHSCOMP, NRHS, LRHSCOMP, FIRST_ROW_RHSCOMP, W, LD_W,
1040 & FIRST_ROW_W )
1041 INTEGER :: JBDEB, JBFIN, NBROWS
1042 INTEGER :: NRHS, LRHSCOMP
1043 INTEGER :: FIRST_ROW_RHSCOMP
1044 INTEGER, INTENT(IN) :: KEEP(500)
1045 COMPLEX(kind=8), INTENT(INOUT) :: RHSCOMP(LRHSCOMP,NRHS)
1046 INTEGER :: LD_W, FIRST_ROW_W
1047 COMPLEX(kind=8) :: W(LD_W*(JBFIN-JBDEB+1))
1048 INTEGER :: JJ, K, ISHIFT
1049!$OMP PARALLEL DO PRIVATE(ISHIFT, JJ), IF
1050!$OMP& (JBFIN-JBDEB+1 > 2*KEEP(362) .AND.
1051!$OMP& NBROWS * (JBFIN-JBDEB+1) > 2*KEEP(363))
1052 DO k = jbdeb, jbfin
1053 ishift = first_row_w + ld_w * (k-jbdeb)
1054 DO jj = 0, nbrows-1
1055 rhscomp(first_row_rhscomp+jj,k) = w(ishift+jj)
1056 END DO
1057 END DO
1058!$OMP END PARALLEL DO
1059 RETURN
1060 END SUBROUTINE zmumps_sol_cpy_fs2rhscomp
1061 SUBROUTINE zmumps_sol_bwd_gthr( JBDEB, JBFIN, J1, J2,
1062 & RHSCOMP, NRHS, LRHSCOMP, W, LD_W, FIRST_ROW_W,
1063 & IW, LIW, KEEP, N, POSINRHSCOMP_BWD )
1064 INTEGER, INTENT(IN) :: JBDEB, JBFIN, J1, J2
1065 INTEGER, INTENT(IN) :: NRHS, LRHSCOMP
1066 INTEGER, INTENT(IN) :: FIRST_ROW_W, LD_W, LIW
1067 INTEGER, INTENT(IN) :: IW(LIW)
1068 INTEGER, INTENT(IN) :: KEEP(500)
1069 COMPLEX(kind=8), INTENT(INOUT) :: RHSCOMP(LRHSCOMP,NRHS)
1070 COMPLEX(kind=8) :: W(LD_W*(JBFIN-JBDEB+1))
1071 INTEGER, INTENT(IN) :: N
1072 INTEGER, INTENT(IN) :: POSINRHSCOMP_BWD(N)
1073 INTEGER :: ISHIFT, JJ, K, IPOSINRHSCOMP
1074!$OMP PARALLEL DO PRIVATE(JJ,ISHIFT,IPOSINRHSCOMP), IF
1075!$OMP& ((JBFIN-JBDEB+1 > 2*KEEP(362) .AND.
1076!$OMP& (JBFIN-JBDEB+1)*(J2-KEEP(253)-J1+1)>2*KEEP(363)))
1077 DO k=jbdeb, jbfin
1078 ishift = first_row_w+(k-jbdeb)*ld_w
1079 DO jj = j1, j2-keep(253)
1080 iposinrhscomp = abs(posinrhscomp_bwd(iw(jj)))
1081 w(ishift+jj-j1)= rhscomp(iposinrhscomp,k)
1082 ENDDO
1083 ENDDO
1084!$OMP END PARALLEL DO
1085 RETURN
1086 END SUBROUTINE zmumps_sol_bwd_gthr
1087 SUBROUTINE zmumps_sol_q(MTYPE, IFLAG, N,
1088 & LHS, WRHS, W, RES, GIVNORM, ANORM, XNORM, SCLNRM,
1089 & MPRINT, ICNTL, KEEP,KEEP8)
1090 INTEGER MTYPE,N,IFLAG,ICNTL(60), KEEP(500)
1091 INTEGER(8) KEEP8(150)
1092 COMPLEX(kind=8) RES(N),LHS(N)
1093 COMPLEX(kind=8) WRHS(N)
1094 DOUBLE PRECISION W(N)
1095 DOUBLE PRECISION RESMAX,RESL2,XNORM, SCLNRM
1096 DOUBLE PRECISION ANORM,DZERO
1097 LOGICAL GIVNORM,PROK
1098 INTEGER MPRINT, MP
1099 INTEGER K
1100 INTRINSIC abs, max, sqrt
1101 mp = icntl(2)
1102 prok = (mprint .GT. 0)
1103 dzero = 0.0d0
1104 IF (.NOT.givnorm) anorm = dzero
1105 resmax = dzero
1106 resl2 = dzero
1107 DO 40 k = 1, n
1108 resmax = max(resmax, abs(res(k)))
1109 resl2 = resl2 + abs(res(k)) * abs(res(k))
1110 IF (.NOT.givnorm) anorm = max(anorm, w(k))
1111 40 CONTINUE
1112 xnorm = dzero
1113 DO 50 k = 1, n
1114 xnorm = max(xnorm, abs(lhs(k)))
1115 50 CONTINUE
1116 IF ( xnorm .EQ. dzero .OR. (exponent(xnorm) .LT.
1117 & minexponent(xnorm) + keep(122) )
1118 & .OR.
1119 & ( exponent(anorm)+exponent(xnorm) .LT.
1120 & minexponent(xnorm) + keep(122) )
1121 & .OR.
1122 & ( exponent(anorm) + exponent(xnorm) -exponent(resmax)
1123 & .LT. minexponent(xnorm) + keep(122) )
1124 & ) THEN
1125 IF (mod(iflag/2,2) .EQ. 0) THEN
1126 iflag = iflag + 2
1127 ENDIF
1128 IF ((mp .GT. 0) .AND. (icntl(4) .GE. 2)) WRITE( mp, * )
1129 & ' max-NORM of computed solut. is zero or close to zero. '
1130 ENDIF
1131 IF (resmax .EQ. dzero) THEN
1132 sclnrm = dzero
1133 ELSE
1134 sclnrm = resmax / (anorm * xnorm)
1135 ENDIF
1136 resl2 = sqrt(resl2)
1137 IF (prok) WRITE( mprint, 90 ) resmax, resl2, anorm, xnorm,
1138 & sclnrm
1139 90 FORMAT (/' RESIDUAL IS ............ (MAX-NORM) =',1pd9.2/
1140 & ' .. (2-NORM) =',1pd9.2/
1141 & ' RINFOG(4):NORM OF input Matrix (MAX-NORM)=',1pd9.2/
1142 & ' RINFOG(5):NORM OF Computed SOLUT (MAX-NORM)=',1pd9.2/
1143 & ' RINFOG(6):SCALED RESIDUAL ...... (MAX-NORM)=',1pd9.2)
1144 RETURN
1145 END SUBROUTINE zmumps_sol_q
1146 SUBROUTINE zmumps_solve_fwd_trsolve (A, LA, APOS, NPIV, LDADIAG,
1147 & NRHS_B, WCB, LWCB, LDA_WCB, PPIV_COURANT, MTYPE, KEEP)
1148 INTEGER, INTENT(IN) :: MTYPE, LDADIAG, NPIV, KEEP(500)
1149 INTEGER, INTENT(IN) :: NRHS_B, LDA_WCB
1150 INTEGER(8), INTENT(IN) :: LA, APOS, LWCB, PPIV_COURANT
1151 COMPLEX(kind=8), INTENT(IN) :: A(LA)
1152 COMPLEX(kind=8), INTENT(INOUT) :: WCB(LWCB)
1153 COMPLEX(kind=8) ONE
1154 PARAMETER ( ONE=(1.0d0,0.0d0) )
1155 IF (keep(50).NE.0 .OR. mtype .eq. 1 ) THEN
1156#if defined(MUMPS_USE_BLAS2)
1157 IF ( nrhs_b == 1 ) THEN
1158 CALL ztrsv( 'U', 'T', 'U', npiv, a(apos), ldadiag,
1159 & wcb(ppiv_courant), 1 )
1160 ELSE
1161#endif
1162 CALL ztrsm( 'L','U','T','U', npiv, nrhs_b, one,
1163 & a(apos), ldadiag, wcb(ppiv_courant),
1164 & lda_wcb )
1165#if defined(MUMPS_USE_BLAS2)
1166 ENDIF
1167#endif
1168 ELSE
1169#if defined(MUMPS_USE_BLAS2)
1170 IF ( nrhs_b == 1 ) THEN
1171 CALL ztrsv( 'L', 'N', 'N', npiv, a(apos), ldadiag,
1172 & wcb(ppiv_courant), 1 )
1173 ELSE
1174#endif
1175 CALL ztrsm( 'L','L','N','N', npiv, nrhs_b, one,
1176 & a(apos), ldadiag, wcb(ppiv_courant),
1177 & lda_wcb )
1178#if defined(MUMPS_USE_BLAS2)
1179 ENDIF
1180#endif
1181 ENDIF
1182 RETURN
1183 END SUBROUTINE zmumps_solve_fwd_trsolve
1184 SUBROUTINE zmumps_solve_bwd_trsolve (A, LA, APOS, NPIV, LDADIAG,
1185 & NRHS_B, WCB, LWCB, LDA_WCB, PPIV_COURANT, MTYPE, KEEP)
1186 INTEGER, INTENT(IN) :: MTYPE, LDADIAG, NPIV, KEEP(500)
1187 INTEGER, INTENT(IN) :: NRHS_B, LDA_WCB
1188 INTEGER(8), INTENT(IN) :: LA, APOS, LWCB, PPIV_COURANT
1189 COMPLEX(kind=8), INTENT(IN) :: A(LA)
1190 COMPLEX(kind=8), INTENT(INOUT) :: WCB(LWCB)
1191 COMPLEX(kind=8) ONE
1192 PARAMETER ( ONE=(1.0d0,0.0d0) )
1193 IF (mtype .eq. 1 ) THEN
1194#if defined(MUMPS_USE_BLAS2)
1195 IF ( nrhs_b == 1 ) THEN
1196 CALL ztrsv( 'L', 'T', 'N', npiv, a(apos), ldadiag,
1197 & wcb(ppiv_courant), 1 )
1198 ELSE
1199#endif
1200 CALL ztrsm( 'L','L','T','N', npiv, nrhs_b, one,
1201 & a(apos), ldadiag, wcb(ppiv_courant),
1202 & lda_wcb )
1203#if defined(MUMPS_USE_BLAS2)
1204 ENDIF
1205#endif
1206 ELSE
1207#if defined(MUMPS_USE_BLAS2)
1208 IF ( nrhs_b == 1 ) THEN
1209 CALL ztrsv( 'U', 'N', 'U', npiv, a(apos), ldadiag,
1210 & wcb(ppiv_courant), 1 )
1211 ELSE
1212#endif
1213 CALL ztrsm( 'L','U','N','U', npiv, nrhs_b, one,
1214 & a(apos), ldadiag, wcb(ppiv_courant),
1215 & lda_wcb )
1216#if defined(MUMPS_USE_BLAS2)
1217 ENDIF
1218#endif
1219 ENDIF
1220 RETURN
1221 END SUBROUTINE zmumps_solve_bwd_trsolve
1222 SUBROUTINE zmumps_solve_fwd_panels(
1223 & A, LA, APOS, NPIV, IW,
1224 & NRHS_B, WCB, LWCB, LDA_WCB,
1225 & PPIV_COURANT, MTYPE, KEEP)
1226 INTEGER, INTENT(IN) :: MTYPE, NPIV, KEEP(500)
1227 INTEGER, INTENT(IN) :: IW(NPIV)
1228 INTEGER, INTENT(IN) :: NRHS_B, LDA_WCB
1229 INTEGER(8), INTENT(IN) :: LA, APOS, LWCB, PPIV_COURANT
1230 COMPLEX(kind=8), INTENT(IN) :: A(LA)
1231 COMPLEX(kind=8), INTENT(INOUT) :: WCB(LWCB)
1232 INTEGER :: NB_TARGET
1233 INTEGER :: NBPANELS
1234 INTEGER :: NBROWS_PANEL, NBCOLS_PANEL, ICOL_BEG, ICOL_END
1235 INTEGER(8) :: PANEL_APOS, PPIV_PANEL
1236 COMPLEX(kind=8), PARAMETER :: ONE=(1.0d0,0.0d0)
1237 IF (keep(459) .LE. 1) THEN
1238 WRITE(*,*) " Internal error in ZMUMPS_SOLVE_FWD_PANELS"
1239 CALL mumps_abort()
1240 ENDIF
1241 CALL mumps_ldltpanel_nbtarget( npiv, nb_target, keep )
1242 panel_apos = apos
1243 nbpanels = 0
1244 icol_beg = 1
1245 nbrows_panel = npiv
1246 ppiv_panel = ppiv_courant
1247 DO WHILE ( icol_beg .LE. npiv )
1248 nbpanels = nbpanels + 1
1249 icol_end = min(nb_target * nbpanels, npiv)
1250 IF ( iw(icol_end) .LT. 0 ) icol_end=icol_end+1
1251 nbcols_panel = icol_end - icol_beg + 1
1252 CALL zmumps_solve_fwd_trsolve (a, la, panel_apos,
1253 & nbcols_panel, nbcols_panel,
1254 & nrhs_b, wcb, lwcb, lda_wcb, ppiv_panel, mtype, keep)
1255 IF ( nbrows_panel .GT. nbcols_panel ) THEN
1256 CALL zmumps_solve_gemm_update( a, la,
1257 & panel_apos + int(nbcols_panel,8) * int(nbcols_panel,8),
1258 & nbcols_panel, nbcols_panel, nbrows_panel-nbcols_panel,
1259 & nrhs_b, wcb, lwcb, ppiv_panel, lda_wcb,
1260 & ppiv_panel+nbcols_panel, lda_wcb,
1261 & mtype, keep, one )
1262 ENDIF
1263 icol_beg = icol_end + 1
1264 panel_apos = panel_apos + int(nbcols_panel,8) *
1265 & int(nbrows_panel,8)
1266 nbrows_panel = nbrows_panel - nbcols_panel
1267 ppiv_panel = ppiv_panel + nbcols_panel
1268 ENDDO
1269 RETURN
1270 END SUBROUTINE zmumps_solve_fwd_panels
1271 SUBROUTINE zmumps_solve_bwd_panels(
1272 & A, LA, APOS, NPIV, IW,
1273 & NRHS_B, WCB, LWCB, LDA_WCB,
1274 & PPIV_COURANT, MTYPE, KEEP)
1275 INTEGER, INTENT(IN) :: MTYPE, NPIV, KEEP(500)
1276 INTEGER, INTENT(IN) :: IW(NPIV)
1277 INTEGER, INTENT(IN) :: NRHS_B, LDA_WCB
1278 INTEGER(8), INTENT(IN) :: LA, APOS, LWCB, PPIV_COURANT
1279 COMPLEX(kind=8), INTENT(IN) :: A(LA)
1280 COMPLEX(kind=8), INTENT(INOUT) :: WCB(LWCB)
1281 INTEGER, PARAMETER :: PANEL_TABSIZE = 20
1282 INTEGER(8) :: PANEL_POS(PANEL_TABSIZE)
1283 INTEGER :: PANEL_COL(PANEL_TABSIZE)
1284 INTEGER :: IPANEL, NBPANELS, NB_TARGET
1285 INTEGER :: NBROWS_PANEL, NBCOLS_PANEL
1286 INTEGER(8) :: PPIV_PANEL
1287 INTEGER :: MTYPE_TEMP
1288 COMPLEX(kind=8), PARAMETER :: ONE=(1.0d0,0.0d0)
1289 IF (keep(459) .LE. 1) THEN
1290 WRITE(*,*) " Internal error 1 in ZMUMPS_SOLVE_BWD_PANELS"
1291 CALL mumps_abort()
1292 ENDIF
1293 IF ( keep(459)+1 .GT. panel_tabsize ) THEN
1294 WRITE(*,*) " Internal error 2 in ZMUMPS_SOLVE_BWD_PANELS"
1295 CALL mumps_abort()
1296 ENDIF
1297 CALL mumps_ldltpanel_panelinfos( npiv, keep, iw,
1298 &nb_target, nbpanels, panel_col, panel_pos, panel_tabsize,
1299 &.false. )
1300 DO ipanel = nbpanels, 1, -1
1301 nbcols_panel = panel_col( ipanel+1 ) - panel_col( ipanel )
1302 nbrows_panel = npiv - panel_col( ipanel ) + 1
1303 ppiv_panel = ppiv_courant + panel_col( ipanel ) - 1
1304 IF ( nbrows_panel .GT. nbcols_panel ) THEN
1305 mtype_temp = 0
1306 CALL zmumps_solve_gemm_update( a, la,
1307 & apos-1_8+panel_pos(ipanel)+
1308 & int(nbcols_panel,8)*int(nbcols_panel,8),
1309 & nbrows_panel-nbcols_panel, nbcols_panel,
1310 & nbcols_panel,
1311 & nrhs_b, wcb, lwcb, ppiv_panel+nbcols_panel, lda_wcb,
1312 & ppiv_panel, lda_wcb,
1313 & mtype_temp, keep, one )
1314 ENDIF
1315 CALL zmumps_solve_bwd_trsolve (a, la,
1316 & apos+panel_pos(ipanel)-1_8,
1317 & nbcols_panel, nbcols_panel,
1318 & nrhs_b, wcb, lwcb, lda_wcb, ppiv_panel, mtype, keep)
1319 ENDDO
1320 RETURN
1321 END SUBROUTINE zmumps_solve_bwd_panels
1322 SUBROUTINE zmumps_solve_gemm_update
1323 & (a, la, apos1, nx, lda, ny,
1324 & nrhs_b, wcb, lwcb, ptrx, ldx,
1325 & ptry, ldy,
1326 & mtype, keep, coef_y )
1327 INTEGER, INTENT(IN) :: MTYPE, NY, NX, KEEP(500)
1328 INTEGER, INTENT(IN) :: NRHS_B, LDY, LDA, LDX
1329 INTEGER(8), INTENT(IN) :: LA, APOS1, LWCB, PTRX,
1330 & ptry
1331 COMPLEX(kind=8), INTENT(IN) :: A(LA)
1332 COMPLEX(kind=8), INTENT(INOUT) :: WCB(LWCB)
1333 COMPLEX(kind=8), INTENT(IN) :: COEF_Y
1334 COMPLEX(kind=8) ALPHA, ZERO, ONE
1335 PARAMETER (ZERO=(0.0d0,0.0d0), one=(1.0d0,0.0d0),
1336 & alpha=(-1.0d0,0.0d0))
1337 IF ( nx .NE. 0 .AND. ny.NE.0 ) THEN
1338 IF ( mtype .eq. 1 ) THEN
1339#if defined(MUMPS_USE_BLAS2)
1340 IF ( nrhs_b == 1 ) THEN
1341 CALL zgemv('T', nx, ny, alpha, a(apos1),
1342 & lda, wcb(ptrx), 1, coef_y,
1343 & wcb(ptry), 1)
1344 ELSE
1345#endif
1346 CALL zgemm('T', 'N', ny, nrhs_b, nx, alpha,
1347 & a(apos1), lda, wcb(ptrx), ldx, coef_y,
1348 & wcb(ptry), ldy)
1349#if defined(MUMPS_USE_BLAS2)
1350 END IF
1351#endif
1352 ELSE
1353#if defined(mumps_use_blas2)
1354 IF ( nrhs_b == 1 ) THEN
1355 CALL zgemv('N',ny, nx, alpha, a(apos1),
1356 & lda, wcb(ptrx), 1,
1357 & coef_y, wcb(ptry), 1 )
1358 ELSE
1359#endif
1360 CALL zgemm('N', 'N', ny, nrhs_b, nx, alpha,
1361 & a(apos1), lda, wcb(ptrx), ldx,
1362 & coef_y, wcb(ptry), ldy)
1363#if defined(MUMPS_USE_BLAS2)
1364 END IF
1365#endif
1366 END IF
1367 END IF
1368 RETURN
1369 END SUBROUTINE zmumps_solve_gemm_update
1370 SUBROUTINE zmumps_sol_ld_and_reload_panel (
1371 & INODE, N, NPIV, LIELL, NELIM, NSLAVES,
1372 & PPIV_COURANT,
1373 & IW, IPOS, LIW,
1374 & A, LA, APOS,
1375 & WCB, LWCB, LD_WCBPIV,
1376 & RHSCOMP, LRHSCOMP, NRHS,
1377 & POSINRHSCOMP_FWD, JBDEB, JBFIN,
1378 & MTYPE, KEEP, OOCWRITE_COMPATIBLE_WITH_BLR,
1379 & IGNORE_K459
1380 & )
1381 USE zmumps_ooc
1382 IMPLICIT NONE
1383 INTEGER, INTENT(IN) :: MTYPE, INODE, N, NPIV, LIELL,
1384 & nelim, nslaves
1385 INTEGER, INTENT(IN) :: LRHSCOMP, NRHS, LIW, JBDEB, JBFIN
1386 INTEGER, INTENT(IN) :: IW(LIW), IPOS, POSINRHSCOMP_FWD(N)
1387 INTEGER(8), INTENT(IN) :: LWCB, APOS, LA, PPIV_COURANT
1388 INTEGER, INTENT(IN) :: LD_WCBPIV
1389 INTEGER, INTENT(IN) :: KEEP(500)
1390 COMPLEX(kind=8), INTENT(IN) :: WCB( LWCB )
1391 COMPLEX(kind=8), INTENT(IN) :: A( LA )
1392 COMPLEX(kind=8), INTENT(INOUT) :: RHSCOMP(LRHSCOMP, NRHS)
1393 LOGICAL, INTENT(IN) :: OOCWRITE_COMPATIBLE_WITH_BLR
1394 LOGICAL, INTENT(IN) :: IGNORE_K459
1395 INTEGER :: J1, J3
1396 INTEGER :: IPOSINRHSCOMP, JJ, K, NBK,
1397 & ldaj, nrhs_b
1398 INTEGER(8) :: IFR8 , APOS1, APOS2, APOSOFF, IFR_ini8,
1399 & poswcb1, poswcb2
1400 COMPLEX(kind=8) :: VALPIV, A11, A22, A12, DETPIV
1401 INTEGER, PARAMETER :: PANEL_TABSIZE = 20
1402 INTEGER(8) :: PANEL_POS(PANEL_TABSIZE)
1403 INTEGER :: PANEL_COL(PANEL_TABSIZE)
1404 INTEGER :: IPANEL, ICOL, NBPANELS, NB_TARGET
1405 LOGICAL :: SKIP_IT
1406 LOGICAL :: OMP_FLAG
1407 COMPLEX(kind=8) ONE
1408 PARAMETER ( ONE=(1.0d0,0.0d0) )
1409 IF ( npiv.EQ. 0 ) RETURN
1410 nrhs_b = jbfin-jbdeb+1
1411 IF ( mtype .EQ. 1 .OR. keep(50) .NE. 0 ) THEN
1412 j1 = ipos + 1
1413 j3 = ipos + npiv
1414 ELSE
1415 j1 = ipos + liell + 1
1416 j3 = ipos + liell + npiv
1417 END IF
1418 iposinrhscomp = posinrhscomp_fwd(iw(j1))
1419 IF ( keep(50) .eq. 0 ) THEN
1420 omp_flag = .false.
1421!$ OMP_FLAG=(int(NRHS_B,8)*int(NPIV,8).GE.int(KEEP(363),8))
1422 IF (omp_flag) THEN
1423!$OMP PARALLEL DO PRIVATE(IFR8) COLLAPSE(2)
1424 DO k = jbdeb, jbfin
1425 DO ifr8 = 0_8, int(npiv-1,8)
1426 rhscomp(iposinrhscomp+ifr8, k) =
1427 & wcb(ppiv_courant+(k-jbdeb)*ld_wcbpiv+ifr8)
1428 ENDDO
1429 ENDDO
1430!$OMP END PARALLEL DO
1431 ELSE
1432 DO k = jbdeb, jbfin
1433 DO ifr8 = 0_8, int(npiv-1,8)
1434 rhscomp(iposinrhscomp+ifr8, k) =
1435 & wcb(ppiv_courant+(k-jbdeb)*ld_wcbpiv+ifr8)
1436 ENDDO
1437 ENDDO
1438 ENDIF
1439 ELSE
1440 CALL mumps_ldltpanel_panelinfos( npiv, keep, iw(ipos+liell+1),
1441 & nb_target, nbpanels, panel_col, panel_pos, panel_tabsize,
1442 & ignore_k459 )
1443 ifr_ini8 = ppiv_courant
1444!$ OMP_FLAG = ( JBFIN-JBDEB+1.GE.KEEP(362) .AND.
1445!$ & ((J3-J1+1)*(JBFIN-JBDEB+1) .GE. KEEP(363)))
1446!$OMP PARALLEL DO PRIVATE(JJ,IFR8,NBK,APOS1,APOS2,APOSOFF,VALPIV,
1447!$OMP& IPANEL,ICOL,
1448!$OMP& POSWCB1,POSWCB2,A11,A22,A12,DETPIV,LDAJ,SKIP_IT)
1449!$OMP& IF(OMP_FLAG)
1450 DO k = jbdeb, jbfin
1451 DO jj = j1, j3
1452 ipanel = (jj-j1)/nb_target + 1
1453 IF ( jj-j1+1 .LT. panel_col(ipanel) ) ipanel = ipanel -1
1454 icol = jj-j1+1 - panel_col(ipanel) + 1
1455 ldaj = panel_col(ipanel+1) - panel_col(ipanel)
1456 apos1 = apos-1_8+panel_pos( ipanel ) + int(icol-1,8) *
1457 & int(ldaj+1,8)
1458 ifr8 = ifr_ini8 + int(k-jbdeb,8)*int(ld_wcbpiv,8) +
1459 & int(jj-j1,8)
1460 IF ( jj .NE. j1 ) THEN
1461 IF ( iw(liell+jj-1) .LT. 0 ) THEN
1462 skip_it = .true.
1463 ELSE
1464 skip_it = .false.
1465 ENDIF
1466 ELSE
1467 skip_it = .false.
1468 ENDIF
1469 IF (skip_it) THEN
1470 ELSE IF ( iw(jj+liell) .GT. 0 ) THEN
1471 valpiv = one/a( apos1 )
1472 rhscomp(iposinrhscomp+jj-j1 , k ) =
1473 & wcb( ifr8 ) * valpiv
1474 apos1 = apos1 + int(ldaj + 1,8)
1475 ELSE
1476 apos2 = apos1+int(ldaj+1,8)
1477 aposoff=apos1+1_8
1478 a11 = a(apos1)
1479 a22 = a(apos2)
1480 a12 = a(aposoff)
1481 detpiv = a11*a22 - a12**2
1482 a22 = a11/detpiv
1483 a11 = a(apos2)/detpiv
1484 a12 = -a12/detpiv
1485 poswcb1 = ifr8
1486 poswcb2 = poswcb1+1_8
1487 rhscomp(iposinrhscomp+jj-j1,k) =
1488 & wcb(poswcb1)*a11
1489 & + wcb(poswcb2)*a12
1490 rhscomp(iposinrhscomp+jj-j1+1,k) =
1491 & wcb(poswcb1)*a12
1492 & + wcb(poswcb2)*a22
1493 ENDIF
1494 ENDDO
1495 ENDDO
1496!$OMP END PARALLEL DO
1497 END IF
1498 RETURN
1499 END SUBROUTINE zmumps_sol_ld_and_reload_panel
1500 SUBROUTINE zmumps_sol_ld_and_reload (
1501 & INODE, N, NPIV, LIELL, NELIM, NSLAVES,
1502 & PPIV_COURANT,
1503 & IW, IPOS, LIW,
1504 & A, LA, APOS,
1505 & WCB, LWCB, LD_WCBPIV,
1506 & RHSCOMP, LRHSCOMP, NRHS,
1507 & POSINRHSCOMP_FWD, JBDEB, JBFIN,
1508 & MTYPE, KEEP, OOCWRITE_COMPATIBLE_WITH_BLR,
1509 & IGNORE_K459
1510 & )
1511 USE zmumps_ooc
1512 INTEGER, INTENT(IN) :: MTYPE, INODE, N, NPIV, LIELL,
1513 & nelim, nslaves
1514 INTEGER, INTENT(IN) :: LRHSCOMP, NRHS, LIW, JBDEB, JBFIN
1515 INTEGER, INTENT(IN) :: IW(LIW), IPOS, POSINRHSCOMP_FWD(N)
1516 INTEGER(8), INTENT(IN) :: LWCB, APOS, LA, PPIV_COURANT
1517 INTEGER, INTENT(IN) :: LD_WCBPIV
1518 INTEGER, INTENT(IN) :: KEEP(500)
1519 COMPLEX(kind=8), INTENT(IN) :: WCB( LWCB )
1520 COMPLEX(kind=8), INTENT(IN) :: A( LA )
1521 COMPLEX(kind=8), INTENT(INOUT) :: RHSCOMP(LRHSCOMP, NRHS)
1522 LOGICAL, INTENT(IN) :: OOCWRITE_COMPATIBLE_WITH_BLR
1523 LOGICAL, INTENT(IN) :: IGNORE_K459
1524 INTEGER :: TempNROW, J1, J3, PANEL_SIZE, TYPEF
1525 INTEGER :: IPOSINRHSCOMP, JJ, K, NBK, LDAJ,
1526 & ldaj_ini, nbk_ini, ldaj_first_panel, nrhs_b
1527 INTEGER(8) :: IFR8 , APOS1, APOS2, APOSOFF, IFR_ini8,
1528 & poswcb1, poswcb2
1529 COMPLEX(kind=8) :: VALPIV, A11, A22, A12, DETPIV
1530!$ LOGICAL :: OMP_FLAG
1531 COMPLEX(kind=8) ONE
1532 parameter( one=(1.0d0,0.0d0) )
1533 nrhs_b = jbfin-jbdeb+1
1534 IF ( mtype .EQ. 1 .OR. keep(50) .NE. 0 ) THEN
1535 j1 = ipos + 1
1536 j3 = ipos + npiv
1537 ELSE
1538 j1 = ipos + liell + 1
1539 j3 = ipos + liell + npiv
1540 END IF
1541 iposinrhscomp = posinrhscomp_fwd(iw(j1))
1542 IF ( keep(50) .eq. 0 ) THEN
1543!$ OMP_FLAG=(NRHS_B.GE.KEEP(362).AND.NRHS_B*NPIV.GE.KEEP(363))
1544!$OMP PARALLEL DO PRIVATE(IFR8) IF (OMP_FLAG)
1545 DO k=jbdeb,jbfin
1546 ifr8 = ppiv_courant + (k-jbdeb)*ld_wcbpiv
1547 rhscomp(iposinrhscomp:iposinrhscomp+npiv-1, k) =
1548 & wcb(ifr8:ifr8+int(npiv-1,8))
1549 ENDDO
1550!$OMP END PARALLEL DO
1551 ELSE
1552 ifr8 = ppiv_courant - 1_8
1553 IF (keep(201).EQ.1.AND.oocwrite_compatible_with_blr) THEN
1554 IF (mtype.EQ.1) THEN
1555 IF ((mtype.EQ.1).AND.nslaves.NE.0) THEN
1556 tempnrow= npiv+nelim
1557 ldaj_first_panel=tempnrow
1558 ELSE
1559 tempnrow= liell
1560 ldaj_first_panel=tempnrow
1561 ENDIF
1562 typef=typef_l
1563 ELSE
1564 tempnrow= npiv
1565 ldaj_first_panel=liell
1566 typef= typef_u
1567 ENDIF
1568 panel_size = zmumps_ooc_panel_size( ldaj_first_panel )
1569 ldaj = tempnrow
1570 ELSE
1571 IF ( keep(459) .GT. 1 .AND. keep(50) .NE. 0
1572 & .AND. .NOT. ignore_k459 ) THEN
1573 CALL mumps_ldltpanel_nbtarget( npiv, panel_size, keep )
1574 ldaj = panel_size
1575 ELSE
1576 panel_size = -1
1577 ldaj = npiv
1578 ENDIF
1579 ENDIF
1580 IF (keep(201).EQ.1.AND.oocwrite_compatible_with_blr) THEN
1581 nbk = 0
1582 ENDIF
1583 ifr_ini8 = ppiv_courant - 1_8
1584 ldaj_ini = ldaj
1585 IF (keep(201).EQ.1.AND.oocwrite_compatible_with_blr)
1586 & nbk_ini = nbk
1587!$ OMP_FLAG = ( JBFIN-JBDEB+1.GE.KEEP(362) .AND.
1588!$ & ((J3-J1+1)*(JBFIN-JBDEB+1) .GE. KEEP(363)))
1589!$OMP PARALLEL DO PRIVATE(JJ,IFR8,NBK,APOS1,APOS2,APOSOFF,VALPIV,
1590!$OMP& POSWCB1, POSWCB2,A11,A22,A12,DETPIV,LDAJ) IF(OMP_FLAG)
1591 DO k = jbdeb, jbfin
1592 ifr8 = ifr_ini8 + int(k-jbdeb,8)*int(ld_wcbpiv,8)
1593 nbk = nbk_ini
1594 apos1 = apos
1595 ldaj = ldaj_ini
1596 jj = j1
1597 DO
1598 IF (jj .GT. j3) EXIT
1599 ifr8 = ifr8 + 1_8
1600 IF (iw(jj+liell) .GT. 0) THEN
1601 valpiv = one/a( apos1 )
1602 rhscomp(iposinrhscomp+jj-j1 , k ) =
1603 & wcb( ifr8 ) * valpiv
1604 IF (keep(201).EQ.1.AND.oocwrite_compatible_with_blr)
1605 & THEN
1606 nbk = nbk+1
1607 IF (nbk.EQ.panel_size) THEN
1608 nbk = 0
1609 ldaj = ldaj - panel_size
1610 ENDIF
1611 ENDIF
1612 apos1 = apos1 + int(ldaj + 1,8)
1613 jj = jj+1
1614 ELSE
1615 IF (keep(201).EQ.1.AND.oocwrite_compatible_with_blr)
1616 & THEN
1617 nbk = nbk+1
1618 ENDIF
1619 apos2 = apos1+int(ldaj+1,8)
1620 IF (keep(201).EQ.1.AND.oocwrite_compatible_with_blr)
1621 & THEN
1622 aposoff = apos1+int(ldaj,8)
1623 ELSE
1624 aposoff=apos1+1_8
1625 ENDIF
1626 a11 = a(apos1)
1627 a22 = a(apos2)
1628 a12 = a(aposoff)
1629 detpiv = a11*a22 - a12**2
1630 a22 = a11/detpiv
1631 a11 = a(apos2)/detpiv
1632 a12 = -a12/detpiv
1633 poswcb1 = ifr8
1634 poswcb2 = poswcb1+1_8
1635 rhscomp(iposinrhscomp+jj-j1,k) =
1636 & wcb(poswcb1)*a11
1637 & + wcb(poswcb2)*a12
1638 rhscomp(iposinrhscomp+jj-j1+1,k) =
1639 & wcb(poswcb1)*a12
1640 & + wcb(poswcb2)*a22
1641 IF (keep(201).EQ.1.AND.oocwrite_compatible_with_blr)
1642 & THEN
1643 nbk = nbk+1
1644 IF (nbk.GE.panel_size) THEN
1645 ldaj = ldaj - nbk
1646 nbk = 0
1647 ENDIF
1648 ENDIF
1649 apos1 = apos2 + int(ldaj + 1,8)
1650 jj = jj+2
1651 ifr8 = ifr8+1_8
1652 ENDIF
1653 ENDDO
1654 ENDDO
1655!$OMP END PARALLEL DO
1656 END IF
1657 RETURN
1658 END SUBROUTINE zmumps_sol_ld_and_reload
1659 SUBROUTINE zmumps_set_scaling_loc( scaling_data, N, ILOC, LILOC,
1660 & COMM, MYID, I_AM_SLAVE, MASTER, NB_BYTES, NB_BYTES_MAX,
1661 & K16_8, LP, LPOK, ICNTL, INFO )
1662 IMPLICIT NONE
1663 type scaling_data_t
1664 sequence
1665 DOUBLE PRECISION, dimension(:), pointer :: SCALING
1666 DOUBLE PRECISION, dimension(:), pointer :: SCALING_LOC
1667 end type scaling_data_t
1668 type (scaling_data_t), INTENT(INOUT) :: scaling_data
1669 INTEGER, INTENT(IN) :: N, LILOC, COMM, MYID, MASTER, LP
1670 INTEGER, INTENT(IN) :: ILOC(LILOC)
1671 INTEGER(8), INTENT(INOUT) :: NB_BYTES, NB_BYTES_MAX
1672 INTEGER(8), INTENT(IN) :: K16_8
1673 LOGICAL, INTENT(IN) :: I_AM_SLAVE, LPOK
1674 INTEGER, INTENT(INOUT) :: INFO(80)
1675 INTEGER, INTENT(IN) :: ICNTL(60)
1676 DOUBLE PRECISION, POINTER, DIMENSION(:) :: SCALING
1677 INTEGER :: I, IERR_MPI, allocok
1678 include 'mpif.h'
1679 NULLIFY(scaling_data%SCALING_LOC)
1680 IF (i_am_slave) THEN
1681 ALLOCATE(scaling_data%SCALING_LOC(max(1,liloc)),
1682 & stat=allocok)
1683 IF (allocok > 0) THEN
1684 info(1)=-13
1685 info(2)=max(1,liloc)
1686 GOTO 35
1687 ENDIF
1688 nb_bytes = nb_bytes + int(max(1,liloc),8)*k16_8
1689 nb_bytes_max = max(nb_bytes_max,nb_bytes)
1690 ENDIF
1691 IF (myid .NE. master) THEN
1692 ALLOCATE(scaling(n), stat=allocok)
1693 IF (allocok > 0) THEN
1694 IF (lpok) THEN
1695 WRITE(lp,*) 'Error allocating temporary scaling array'
1696 ENDIF
1697 info(1)=-13
1698 info(2)=n
1699 GOTO 35
1700 ENDIF
1701 nb_bytes = nb_bytes + int(n,8)*k16_8
1702 nb_bytes_max = max(nb_bytes_max,nb_bytes)
1703 ELSE
1704 scaling => scaling_data%SCALING
1705 ENDIF
1706 35 CONTINUE
1707 CALL mumps_propinfo( icntl(1), info(1),
1708 & comm, myid )
1709 IF (info(1) .LT. 0) GOTO 90
1710 CALL mpi_bcast( scaling(1), n, mpi_double_precision,
1711 & master, comm, ierr_mpi)
1712 IF ( i_am_slave ) THEN
1713 DO i = 1, liloc
1714 IF (iloc(i) .GE. 1 .AND. iloc(i) .LE. n) THEN
1715 scaling_data%SCALING_LOC(i) = scaling(iloc(i))
1716 ENDIF
1717 ENDDO
1718 ENDIF
1719 90 CONTINUE
1720 IF (myid.NE. master) THEN
1721 IF (associated(scaling)) THEN
1722 DEALLOCATE(scaling)
1723 nb_bytes = nb_bytes - int(n,8)*k16_8
1724 ENDIF
1725 ENDIF
1726 NULLIFY(scaling)
1727 IF (info(1) .LT. 0) THEN
1728 IF (associated(scaling_data%SCALING_LOC)) THEN
1729 DEALLOCATE(scaling_data%SCALING_LOC)
1730 NULLIFY(scaling_data%SCALING_LOC)
1731 ENDIF
1732 ENDIF
1733 RETURN
1734 END SUBROUTINE zmumps_set_scaling_loc
cmplx
float cmplx[2]
Definition pblas.h:136
mumps_abort
#define mumps_abort
Definition VE_Metis.h:25
mumps_propinfo
subroutine mumps_propinfo(icntl, info, comm, id)
Definition bcast_errors.F:15
zgemv
subroutine zgemv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy)
ZGEMV
Definition zgemv.f:158
ztrsv
subroutine ztrsv(uplo, trans, diag, n, a, lda, x, incx)
ZTRSV
Definition ztrsv.f:149
zgemm
subroutine zgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
ZGEMM
Definition zgemm.f:187
ztrsm
subroutine ztrsm(side, uplo, transa, diag, m, n, alpha, a, lda, b, ldb)
ZTRSM
Definition ztrsm.f:180
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
mpi_recv
subroutine mpi_recv(buf, cnt, datatype, source, tag, comm, status, ierr)
Definition mpi.f:461
mpi_send
subroutine mpi_send(buf, cnt, datatype, dest, tag, comm, ierr)
Definition mpi.f:480
mpi_bcast
subroutine mpi_bcast(buffer, cnt, datatype, root, comm, ierr)
Definition mpi.f:205
zmumps_ooc
Definition zmumps_ooc.F:14
zmumps_ooc::zmumps_solve_is_inode_in_mem
integer function zmumps_solve_is_inode_in_mem(inode, ptrfac, nsteps, a, la, ierr)
Definition zmumps_ooc.F:1326
zmumps_ooc::zmumps_solve_modify_state_node
subroutine zmumps_solve_modify_state_node(inode)
Definition zmumps_ooc.F:1388
zmumps_ooc::zmumps_read_ooc
subroutine, public zmumps_read_ooc(dest, inode, ierr)
Definition zmumps_ooc.F:394
zmumps_ooc::zmumps_ooc_panel_size
integer function, public zmumps_ooc_panel_size(nnmax)
Definition zmumps_ooc.F:3404
zmumps_ooc::ooc_node_permuted
integer ooc_node_permuted
Definition zmumps_ooc.F:24
zmumps_ooc::ooc_node_not_in_mem
integer ooc_node_not_in_mem
Definition zmumps_ooc.F:24
zmumps_ooc::zmumps_solve_alloc_factor_space
subroutine, public zmumps_solve_alloc_factor_space(inode, ptrfac, keep, keep8, a, ierr)
Definition zmumps_ooc.F:1475
zmumps_struc_def
Definition zmumps_struc_def.F:14
mumps_ldltpanel_panelinfos
subroutine mumps_ldltpanel_panelinfos(npiv, keep, iw, nb_target, nbpanels, panel_col, panel_pos, panel_tabsize, ignore_k459)
Definition tools_common.F:1361
mumps_ldltpanel_nbtarget
subroutine mumps_ldltpanel_nbtarget(npiv, nb_target, keep)
Definition tools_common.F:1314
zmumps_eltyd
subroutine zmumps_eltyd(mtype, n, nelt, eltptr, leltvar, eltvar, na_elt8, a_elt, saverhs, x, y, w, k50)
Definition zsol_aux.F:654
zmumps_solve_gemm_update
subroutine zmumps_solve_gemm_update(a, la, apos1, nx, lda, ny, nrhs_b, wcb, lwcb, ptrx, ldx, ptry, ldy, mtype, keep, coef_y)
Definition zsol_aux.F:1327
zmumps_solve_fwd_panels
subroutine zmumps_solve_fwd_panels(a, la, apos, npiv, iw, nrhs_b, wcb, lwcb, lda_wcb, ppiv_courant, mtype, keep)
Definition zsol_aux.F:1226
zmumps_solve_get_ooc_node
subroutine zmumps_solve_get_ooc_node(inode, ptrfac, keep, a, la, step, keep8, n, must_be_permuted, ierr)
Definition zsol_aux.F:733
zmumps_scal_x
subroutine zmumps_scal_x(a, nz8, n, irn, icn, z, keep, keep8, colsca, eff_size_schur, sym_perm)
Definition zsol_aux.F:174
zmumps_solve_bwd_trsolve
subroutine zmumps_solve_bwd_trsolve(a, la, apos, npiv, ldadiag, nrhs_b, wcb, lwcb, lda_wcb, ppiv_courant, mtype, keep)
Definition zsol_aux.F:1186
zmumps_sol_ld_and_reload
subroutine zmumps_sol_ld_and_reload(inode, n, npiv, liell, nelim, nslaves, ppiv_courant, iw, ipos, liw, a, la, apos, wcb, lwcb, ld_wcbpiv, rhscomp, lrhscomp, nrhs, posinrhscomp_fwd, jbdeb, jbfin, mtype, keep, oocwrite_compatible_with_blr, ignore_k459)
Definition zsol_aux.F:1511
zmumps_eltqd2
subroutine zmumps_eltqd2(mtype, n, nelt, eltptr, leltvar, eltvar, na_elt8, a_elt, lhs, wrhs, w, rhs, keep, keep8)
Definition zsol_aux.F:510
zmumps_sol_b
subroutine zmumps_sol_b(n, kase, x, est, w, iw, grain)
Definition zsol_aux.F:304
zmumps_solve_fwd_trsolve
subroutine zmumps_solve_fwd_trsolve(a, la, apos, npiv, ldadiag, nrhs_b, wcb, lwcb, lda_wcb, ppiv_courant, mtype, keep)
Definition zsol_aux.F:1148
zmumps_sol_x
subroutine zmumps_sol_x(a, nz8, n, irn, icn, z, keep, keep8, eff_size_schur, sym_perm)
Definition zsol_aux.F:88
zmumps_sol_y
subroutine zmumps_sol_y(a, nz8, n, irn, icn, rhs, x, r, w, keep, keep8)
Definition zsol_aux.F:224
zmumps_sol_bwd_gthr
subroutine zmumps_sol_bwd_gthr(jbdeb, jbfin, j1, j2, rhscomp, nrhs, lrhscomp, w, ld_w, first_row_w, iw, liw, keep, n, posinrhscomp_bwd)
Definition zsol_aux.F:1064
zmumps_sol_x_elt
subroutine zmumps_sol_x_elt(mtype, n, nelt, eltptr, leltvar, eltvar, na_elt8, a_elt, w, keep, keep8)
Definition zsol_aux.F:530
zmumps_qd2
subroutine zmumps_qd2(mtype, n, nz8, aspk, irn, icn, lhs, wrhs, w, rhs, keep, keep8)
Definition zsol_aux.F:423
zmumps_sol_cpy_fs2rhscomp
subroutine zmumps_sol_cpy_fs2rhscomp(jbdeb, jbfin, nbrows, keep, rhscomp, nrhs, lrhscomp, first_row_rhscomp, w, ld_w, first_row_w)
Definition zsol_aux.F:1041
zmumps_sol_lcond
subroutine zmumps_sol_lcond(n, rhs, x, y, d, r_w, c_w, iw, kase, omega, erx, cond, lp, keep, keep8)
Definition zsol_aux.F:940
zmumps_build_mapping_info
subroutine zmumps_build_mapping_info(id)
Definition zsol_aux.F:776
zmumps_sol_scalx_elt
subroutine zmumps_sol_scalx_elt(mtype, n, nelt, eltptr, leltvar, eltvar, na_elt8, a_elt, w, keep, keep8, colsca)
Definition zsol_aux.F:589
zmumps_sol_omega
subroutine zmumps_sol_omega(n, rhs, x, y, r_w, c_w, iw, iflag, omega, noiter, testconv, lp, arret, grain)
Definition zsol_aux.F:862
zmumps_sol_q
subroutine zmumps_sol_q(mtype, iflag, n, lhs, wrhs, w, res, givnorm, anorm, xnorm, sclnrm, mprint, icntl, keep, keep8)
Definition zsol_aux.F:1090
zmumps_sol_mulr
subroutine zmumps_sol_mulr(n, r, w)
Definition zsol_aux.F:294
zmumps_sol_ld_and_reload_panel
subroutine zmumps_sol_ld_and_reload_panel(inode, n, npiv, liell, nelim, nslaves, ppiv_courant, iw, ipos, liw, a, la, apos, wcb, lwcb, ld_wcbpiv, rhscomp, lrhscomp, nrhs, posinrhscomp_fwd, jbdeb, jbfin, mtype, keep, oocwrite_compatible_with_blr, ignore_k459)
Definition zsol_aux.F:1381
zmumps_solve_bwd_panels
subroutine zmumps_solve_bwd_panels(a, la, apos, npiv, iw, nrhs_b, wcb, lwcb, lda_wcb, ppiv_courant, mtype, keep)
Definition zsol_aux.F:1275
zmumps_compso
subroutine zmumps_compso(n, keep28, iwcb, liww, w, lwc, poswcb, iwposcb, ptricb, ptracb)
Definition zsol_aux.F:37
zmumps_set_scaling_loc
subroutine zmumps_set_scaling_loc(scaling_data, n, iloc, liloc, comm, myid, i_am_slave, master, nb_bytes, nb_bytes_max, k16_8, lp, lpok, icntl, info)
Definition zsol_aux.F:1662
zmumps_freetopso
subroutine zmumps_freetopso(n, keep28, iwcb, liww, w, lwc, poswcb, iwposcb, ptricb, ptracb)
Definition zsol_aux.F:17