OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
dsol_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 dmumps_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 DOUBLE PRECISION 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 dmumps_freetopso
35 SUBROUTINE dmumps_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 DOUBLE PRECISION 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 dmumps_compso
86 SUBROUTINE dmumps_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 DOUBLE PRECISION, 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 dmumps_sol_x
171 SUBROUTINE dmumps_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 DOUBLE PRECISION, 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 dmumps_scal_x
222 SUBROUTINE dmumps_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 DOUBLE PRECISION, INTENT(IN) :: A(NZ8), RHS(N), X(N)
230 DOUBLE PRECISION, INTENT(OUT) :: W(N)
231 DOUBLE PRECISION, INTENT(OUT) :: R(N)
232 INTEGER I, J
233 INTEGER(8) :: K8
234 DOUBLE PRECISION, PARAMETER :: ZERO = 0.0d0
235 DOUBLE PRECISION 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 dmumps_sol_y
293 SUBROUTINE dmumps_sol_mulr(N, R, W)
294 INTEGER, intent(in) :: N
295 DOUBLE PRECISION, intent(in) :: W(N)
296 DOUBLE PRECISION, 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 dmumps_sol_mulr
303 SUBROUTINE dmumps_sol_b(N, KASE, X, EST, W, IW, GRAIN)
304 INTEGER, intent(in) :: N
305 INTEGER, intent(inout) :: KASE
306 INTEGER IW(N)
307 DOUBLE PRECISION W(N), X(N)
308 DOUBLE PRECISION, intent(inout) :: EST
309 INTEGER, intent(in) :: GRAIN
310 INTRINSIC abs, nint, real, sign
311 INTEGER DMUMPS_IXAMAX
312 EXTERNAL dmumps_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 DOUBLE PRECISION ZERO, ONE
320 parameter( zero = 0.0d0 )
321 parameter( one = 1.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) = sign( rone,dble(x(i)) )
353 iw(i) = nint(dble(x(i)))
354 30 CONTINUE
355 kase = 2
356 jump = 2
357 RETURN
358 40 CONTINUE
359 j = dmumps_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) = sign(rone, dble(x(i)))
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 = dmumps_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) = altsgn * (rone + dble(i - 1) / dble(n - 1))
400 altsgn = -altsgn
401 150 CONTINUE
402 kase = 1
403 jump = 5
404 RETURN
405 160 CONTINUE
406 temp = rzero
407 DO 170 i = 1, n
408 temp = temp + abs(x(i))
409 170 CONTINUE
410 temp = 2.0d0 * temp / dble(3 * n)
411 IF (temp .GT. est) THEN
412 DO 180 i = 1, n
413 w(i) = x(i)
414 180 CONTINUE
415 est = temp
416 ENDIF
417 190 kase = 0
418 RETURN
419 END SUBROUTINE dmumps_sol_b
420 SUBROUTINE dmumps_qd2( MTYPE, N, NZ8, ASPK, IRN, ICN,
421 & LHS, WRHS, W, RHS, KEEP,KEEP8)
422 IMPLICIT NONE
423 INTEGER MTYPE, N
424 INTEGER(8), INTENT(IN) :: NZ8
425 INTEGER, INTENT(IN) :: IRN( NZ8 ), ICN( NZ8 )
426 INTEGER KEEP(500)
427 INTEGER(8) KEEP8(150)
428 DOUBLE PRECISION, INTENT(IN) :: ASPK( NZ8 )
429 DOUBLE PRECISION, INTENT(IN) :: LHS( N ), WRHS( N )
430 DOUBLE PRECISION, INTENT(OUT):: RHS( N )
431 DOUBLE PRECISION, INTENT(OUT):: W( N )
432 INTEGER I, J
433 INTEGER(8) :: K8
434 DOUBLE PRECISION, PARAMETER :: DZERO = 0.0d0
435 DO i = 1, n
436 w(i) = dzero
437 rhs(i) = wrhs(i)
438 ENDDO
439 IF ( keep(50) .EQ. 0 ) THEN
440 IF (mtype .EQ. 1) THEN
441 IF (keep(264).EQ.0) THEN
442 DO k8 = 1_8, nz8
443 i = irn(k8)
444 j = icn(k8)
445 IF ((i .LE. 0) .OR. (i .GT. n) .OR. (j .LE. 0) .OR.
446 & (j .GT. n)) cycle
447 rhs(i) = rhs(i) - aspk(k8) * lhs(j)
448 w(i) = w(i) + abs(aspk(k8))
449 ENDDO
450 ELSE
451 DO k8 = 1_8, nz8
452 i = irn(k8)
453 j = icn(k8)
454 rhs(i) = rhs(i) - aspk(k8) * lhs(j)
455 w(i) = w(i) + abs(aspk(k8))
456 ENDDO
457 ENDIF
458 ELSE
459 IF (keep(264).EQ.0) THEN
460 DO k8 = 1_8, nz8
461 i = irn(k8)
462 j = icn(k8)
463 IF ((i .LE. 0) .OR. (i .GT. n) .OR. (j .LE. 0) .OR.
464 & (j .GT. n)) cycle
465 rhs(j) = rhs(j) - aspk(k8) * lhs(i)
466 w(j) = w(j) + abs(aspk(k8))
467 ENDDO
468 ELSE
469 DO k8 = 1_8, nz8
470 i = irn(k8)
471 j = icn(k8)
472 rhs(j) = rhs(j) - aspk(k8) * lhs(i)
473 w(j) = w(j) + abs(aspk(k8))
474 ENDDO
475 ENDIF
476 ENDIF
477 ELSE
478 IF (keep(264).EQ.0) THEN
479 DO k8 = 1_8, nz8
480 i = irn(k8)
481 j = icn(k8)
482 IF ((i .LE. 0) .OR. (i .GT. n) .OR. (j .LE. 0) .OR.
483 & (j .GT. n)) cycle
484 rhs(i) = rhs(i) - aspk(k8) * lhs(j)
485 w(i) = w(i) + abs(aspk(k8))
486 IF (j.NE.i) THEN
487 rhs(j) = rhs(j) - aspk(k8) * lhs(i)
488 w(j) = w(j) + abs(aspk(k8))
489 ENDIF
490 ENDDO
491 ELSE
492 DO k8 = 1_8, nz8
493 i = irn(k8)
494 j = icn(k8)
495 rhs(i) = rhs(i) - aspk(k8) * lhs(j)
496 w(i) = w(i) + abs(aspk(k8))
497 IF (j.NE.i) THEN
498 rhs(j) = rhs(j) - aspk(k8) * lhs(i)
499 w(j) = w(j) + abs(aspk(k8))
500 ENDIF
501 ENDDO
502 ENDIF
503 ENDIF
504 RETURN
505 END SUBROUTINE dmumps_qd2
506 SUBROUTINE dmumps_eltqd2( MTYPE, N,
507 & NELT, ELTPTR, LELTVAR, ELTVAR, NA_ELT8, A_ELT,
508 & LHS, WRHS, W, RHS, KEEP,KEEP8 )
509 IMPLICIT NONE
510 INTEGER MTYPE, N, NELT, LELTVAR
511 INTEGER(8), INTENT(IN) :: NA_ELT8
512 INTEGER ELTPTR(NELT+1), ELTVAR(LELTVAR)
513 INTEGER KEEP(500)
514 INTEGER(8) KEEP8(150)
515 DOUBLE PRECISION A_ELT(NA_ELT8)
516 DOUBLE PRECISION LHS( N ), WRHS( N ), RHS( N )
517 DOUBLE PRECISION W(N)
518 CALL dmumps_mv_elt(n, nelt, eltptr, eltvar, a_elt,
519 & lhs, rhs, keep(50), mtype )
520 rhs = wrhs - rhs
521 CALL dmumps_sol_x_elt( mtype, n,
522 & nelt, eltptr, leltvar, eltvar, na_elt8, a_elt,
523 & w, keep,keep8 )
524 RETURN
525 END SUBROUTINE dmumps_eltqd2
526 SUBROUTINE dmumps_sol_x_elt( MTYPE, N,
527 & NELT, ELTPTR, LELTVAR, ELTVAR, NA_ELT8, A_ELT,
528 & W, KEEP,KEEP8 )
529 IMPLICIT NONE
530 INTEGER MTYPE, N, NELT, LELTVAR
531 INTEGER(8), INTENT(IN) :: NA_ELT8
532 INTEGER ELTPTR(NELT+1), ELTVAR(LELTVAR)
533 INTEGER KEEP(500)
534 INTEGER(8) KEEP8(150)
535 DOUBLE PRECISION A_ELT(NA_ELT8)
536 DOUBLE PRECISION TEMP
537 DOUBLE PRECISION W(N)
538 INTEGER I, J, IEL, SIZEI, IELPTR
539 INTEGER(8) :: K8
540 DOUBLE PRECISION DZERO
541 parameter(dzero = 0.0d0)
542 w = dzero
543 k8 = 1_8
544 DO iel = 1, nelt
545 sizei = eltptr( iel + 1 ) - eltptr( iel )
546 ielptr = eltptr( iel ) - 1
547 IF ( keep(50).EQ.0 ) THEN
548 IF (mtype.EQ.1) THEN
549 DO j = 1, sizei
550 DO i = 1, sizei
551 w( eltvar( ielptr + i) ) =
552 & w( eltvar( ielptr + i) )
553 & + abs(a_elt( k8 ))
554 k8 = k8 + 1_8
555 END DO
556 END DO
557 ELSE
558 DO j = 1, sizei
559 temp = w( eltvar( ielptr + j ) )
560 DO i = 1, sizei
561 temp = temp + abs( a_elt(k8))
562 k8 = k8 + 1_8
563 END DO
564 w(eltvar( ielptr + j )) =
565 & w(eltvar( ielptr + j )) + temp
566 END DO
567 ENDIF
568 ELSE
569 DO j = 1, sizei
570 w(eltvar( ielptr + j )) =
571 & w(eltvar( ielptr + j )) + abs(a_elt( k8 ))
572 k8 = k8 + 1_8
573 DO i = j+1, sizei
574 w(eltvar( ielptr + j )) =
575 & w(eltvar( ielptr + j )) + abs(a_elt( k8 ))
576 w(eltvar( ielptr + i ) ) =
577 & w(eltvar( ielptr + i )) + abs(a_elt( k8 ))
578 k8 = k8 + 1_8
579 END DO
580 ENDDO
581 ENDIF
582 ENDDO
583 RETURN
584 END SUBROUTINE dmumps_sol_x_elt
585 SUBROUTINE dmumps_sol_scalx_elt(MTYPE, N,
586 & NELT, ELTPTR, LELTVAR, ELTVAR, NA_ELT8, A_ELT,
587 & W, KEEP,KEEP8, COLSCA )
588 IMPLICIT NONE
589 INTEGER MTYPE, N, NELT, LELTVAR
590 INTEGER(8), INTENT(IN) :: NA_ELT8
591 INTEGER ELTPTR(NELT+1), ELTVAR(LELTVAR)
592 INTEGER KEEP(500)
593 INTEGER(8) KEEP8(150)
594 DOUBLE PRECISION COLSCA(N)
595 DOUBLE PRECISION A_ELT(NA_ELT8)
596 DOUBLE PRECISION W(N)
597 DOUBLE PRECISION TEMP, TEMP2
598 INTEGER I, J, IEL, SIZEI, IELPTR
599 INTEGER(8) :: K8
600 DOUBLE PRECISION DZERO
601 parameter(dzero = 0.0d0)
602 w = dzero
603 k8 = 1_8
604 DO iel = 1, nelt
605 sizei = eltptr( iel + 1 ) - eltptr( iel )
606 ielptr = eltptr( iel ) - 1
607 IF ( keep(50).EQ.0 ) THEN
608 IF (mtype.EQ.1) THEN
609 DO j = 1, sizei
610 temp2 = abs(colsca(eltvar( ielptr + j) ))
611 DO i = 1, sizei
612 w( eltvar( ielptr + i) ) =
613 & w( eltvar( ielptr + i) )
614 & + abs(a_elt( k8 )) * temp2
615 k8 = k8 + 1_8
616 END DO
617 END DO
618 ELSE
619 DO j = 1, sizei
620 temp = w( eltvar( ielptr + j ) )
621 temp2= abs(colsca(eltvar( ielptr + j) ))
622 DO i = 1, sizei
623 temp = temp + abs(a_elt( k8 )) * temp2
624 k8 = k8 + 1_8
625 END DO
626 w(eltvar( ielptr + j )) =
627 & w(eltvar( ielptr + j )) + temp
628 END DO
629 ENDIF
630 ELSE
631 DO j = 1, sizei
632 w(eltvar( ielptr + j )) =
633 & w(eltvar( ielptr + j )) +
634 & abs( a_elt( k8 )*colsca(eltvar( ielptr + j)) )
635 k8 = k8 + 1_8
636 DO i = j+1, sizei
637 w(eltvar( ielptr + j )) =
638 & w(eltvar( ielptr + j )) +
639 & abs(a_elt( k8 )*colsca(eltvar( ielptr + j)))
640 w(eltvar( ielptr + i ) ) =
641 & w(eltvar( ielptr + i )) +
642 & abs(a_elt( k8 )*colsca(eltvar( ielptr + i)))
643 k8 = k8 + 1_8
644 END DO
645 ENDDO
646 ENDIF
647 ENDDO
648 RETURN
649 END SUBROUTINE dmumps_sol_scalx_elt
650 SUBROUTINE dmumps_eltyd( MTYPE, N, NELT, ELTPTR,
651 & LELTVAR, ELTVAR, NA_ELT8, A_ELT,
652 & SAVERHS, X, Y, W, K50 )
653 IMPLICIT NONE
654 INTEGER N, NELT, K50, MTYPE, LELTVAR
655 INTEGER(8) :: NA_ELT8
656 INTEGER ELTPTR( NELT + 1 ), ELTVAR( LELTVAR )
657 DOUBLE PRECISION A_ELT( NA_ELT8 ), X( N ), Y( N ),
658 & saverhs(n)
659 DOUBLE PRECISION W(N)
660 INTEGER IEL, I , J, K, SIZEI, IELPTR
661 DOUBLE PRECISION ZERO
662 DOUBLE PRECISION TEMP
663 DOUBLE PRECISION TEMP2
664 parameter( zero = 0.0d0 )
665 y = saverhs
666 w = zero
667 k = 1
668 DO iel = 1, nelt
669 sizei = eltptr( iel + 1 ) - eltptr( iel )
670 ielptr = eltptr( iel ) - 1
671 IF ( k50 .eq. 0 ) THEN
672 IF ( mtype .eq. 1 ) THEN
673 DO j = 1, sizei
674 temp = x( eltvar( ielptr + j ) )
675 DO i = 1, sizei
676 y( eltvar( ielptr + i ) ) =
677 & y( eltvar( ielptr + i ) ) -
678 & a_elt( k ) * temp
679 w( eltvar( ielptr + i ) ) =
680 & w( eltvar( ielptr + i ) ) +
681 & abs( a_elt( k ) * temp )
682 k = k + 1
683 END DO
684 END DO
685 ELSE
686 DO j = 1, sizei
687 temp = y( eltvar( ielptr + j ) )
688 temp2 = w( eltvar( ielptr + j ) )
689 DO i = 1, sizei
690 temp = temp -
691 & a_elt( k ) * x( eltvar( ielptr + i ) )
692 temp2 = temp2 + abs(
693 & a_elt( k ) * x( eltvar( ielptr + i ) ) )
694 k = k + 1
695 END DO
696 y( eltvar( ielptr + j ) ) = temp
697 w( eltvar( ielptr + j ) ) = temp2
698 END DO
699 END IF
700 ELSE
701 DO j = 1, sizei
702 y( eltvar( ielptr + j ) ) =
703 & y( eltvar( ielptr + j ) ) -
704 & a_elt( k ) * x( eltvar( ielptr + j ) )
705 w( eltvar( ielptr + j ) ) =
706 & w( eltvar( ielptr + j ) ) + abs(
707 & a_elt( k ) * x( eltvar( ielptr + j ) ) )
708 k = k + 1
709 DO i = j+1, sizei
710 y( eltvar( ielptr + i ) ) =
711 & y( eltvar( ielptr + i ) ) -
712 & a_elt( k ) * x( eltvar( ielptr + j ) )
713 y( eltvar( ielptr + j ) ) =
714 & y( eltvar( ielptr + j ) ) -
715 & a_elt( k ) * x( eltvar( ielptr + i ) )
716 w( eltvar( ielptr + i ) ) =
717 & w( eltvar( ielptr + i ) ) + abs(
718 & a_elt( k ) * x( eltvar( ielptr + j ) ) )
719 w( eltvar( ielptr + j ) ) =
720 & w( eltvar( ielptr + j ) ) + abs(
721 & a_elt( k ) * x( eltvar( ielptr + i ) ) )
722 k = k + 1
723 END DO
724 END DO
725 END IF
726 END DO
727 RETURN
728 END SUBROUTINE dmumps_eltyd
729 SUBROUTINE dmumps_solve_get_ooc_node(
730 & INODE,PTRFAC,KEEP,A,LA,STEP,
731 & KEEP8,N,MUST_BE_PERMUTED,IERR)
732 USE dmumps_ooc
733 IMPLICIT NONE
734 INTEGER INODE,KEEP(500),N
735 INTEGER(8) KEEP8(150)
736 INTEGER(8) :: LA
737 INTEGER(8) :: PTRFAC(KEEP(28))
738 INTEGER STEP(N)
739 INTEGER IERR
740 DOUBLE PRECISION A(LA)
741 INTEGER RETURN_VALUE
742 LOGICAL MUST_BE_PERMUTED
743 return_value=dmumps_solve_is_inode_in_mem(inode,ptrfac,
744 & keep(28),a,la,ierr)
745 IF(return_value.EQ.ooc_node_not_in_mem)THEN
746 IF(ierr.LT.0)THEN
747 RETURN
748 ENDIF
749 CALL dmumps_solve_alloc_factor_space(inode,ptrfac,
750 & keep,keep8,a,ierr)
751 IF(ierr.LT.0)THEN
752 RETURN
753 ENDIF
754 CALL dmumps_read_ooc(
755 & a(ptrfac(step(inode))),
756 & inode,ierr
757 & )
758 IF(ierr.LT.0)THEN
759 RETURN
760 ENDIF
761 ELSE
762 IF(ierr.LT.0)THEN
763 RETURN
764 ENDIF
765 ENDIF
766 IF(return_value.NE.ooc_node_permuted)THEN
767 must_be_permuted=.true.
768 CALL dmumps_solve_modify_state_node(inode)
769 ELSE
770 must_be_permuted=.false.
771 ENDIF
772 RETURN
773 END SUBROUTINE dmumps_solve_get_ooc_node
774 SUBROUTINE dmumps_build_mapping_info(id)
775 USE dmumps_struc_def
776 IMPLICIT NONE
777 include 'mpif.h'
778 TYPE(dmumps_struc), TARGET :: id
779 INTEGER, ALLOCATABLE, DIMENSION(:) :: LOCAL_LIST
780 INTEGER :: I,IERR,TMP,NSTEPS,N_LOCAL_LIST
781 INTEGER :: MASTER,TAG_SIZE,TAG_LIST
782 INTEGER :: STATUS(MPI_STATUS_SIZE)
783 LOGICAL :: I_AM_SLAVE
784 parameter(master=0, tag_size=85,tag_list=86)
785 i_am_slave = (id%MYID .NE. master
786 & .OR. ((id%MYID.EQ.master).AND.(id%KEEP(46).EQ.1)))
787 nsteps = id%KEEP(28)
788 ALLOCATE(local_list(nsteps),stat=ierr)
789 IF(ierr.GT.0) THEN
790 WRITE(*,*)'Problem in solve: error allocating LOCAL_LIST'
791 CALL mumps_abort()
792 END IF
793 n_local_list = 0
794 IF(i_am_slave) THEN
795 DO i=1,nsteps
796 IF(id%PTLUST_S(i).NE.0) THEN
797 n_local_list = n_local_list + 1
798 local_list(n_local_list) = i
799 END IF
800 END DO
801 IF(id%MYID.NE.master) THEN
802 CALL mpi_send(n_local_list, 1,
803 & mpi_integer, master, tag_size, id%COMM,ierr)
804 CALL mpi_send(local_list, n_local_list,
805 & mpi_integer, master, tag_list, id%COMM,ierr)
806 DEALLOCATE(local_list)
807 ALLOCATE(id%IPTR_WORKING(1),
808 & id%WORKING(1),
809 & stat=ierr)
810 IF(ierr.GT.0) THEN
811 WRITE(*,*)'Problem in solve: error allocating ',
812 & 'IPTR_WORKING and WORKING'
813 CALL mumps_abort()
814 END IF
815 END IF
816 END IF
817 IF(id%MYID.EQ.master) THEN
818 ALLOCATE(id%IPTR_WORKING(id%NPROCS+1), stat=ierr)
819 IF(ierr.GT.0) THEN
820 WRITE(*,*)'Problem in solve: error allocating IPTR_WORKING'
821 CALL mumps_abort()
822 END IF
823 id%IPTR_WORKING = 0
824 id%IPTR_WORKING(1) = 1
825 id%IPTR_WORKING(master+2) = n_local_list
826 DO i=1, id%NPROCS-1
827 CALL mpi_recv(tmp, 1, mpi_integer, mpi_any_source,
828 & tag_size, id%COMM, status, ierr)
829 id%IPTR_WORKING(status(mpi_source)+2) = tmp
830 END DO
831 DO i=2, id%NPROCS+1
832 id%IPTR_WORKING(i) = id%IPTR_WORKING(i)
833 & + id%IPTR_WORKING(i-1)
834 END DO
835 ALLOCATE(id%WORKING(id%IPTR_WORKING(id%NPROCS+1)-1),stat=ierr)
836 IF(ierr.GT.0) THEN
837 WRITE(*,*)'Problem in solve: error allocating LOCAL_LIST'
838 CALL mumps_abort()
839 END IF
840 tmp = master + 1
841 IF (i_am_slave) THEN
842 id%WORKING(id%IPTR_WORKING(tmp):id%IPTR_WORKING(tmp+1)-1)
843 & = local_list(1:id%IPTR_WORKING(tmp+1)
844 & -id%IPTR_WORKING(tmp))
845 ENDIF
846 DO i=1,id%NPROCS-1
847 CALL mpi_recv(local_list, nsteps, mpi_integer,
848 & mpi_any_source, tag_list, id%COMM, status, ierr)
849 tmp = status(mpi_source)+1
850 id%WORKING(id%IPTR_WORKING(tmp):id%IPTR_WORKING(tmp+1)-1)
851 & = local_list(1:id%IPTR_WORKING(tmp+1)-
852 & id%IPTR_WORKING(tmp))
853 END DO
854 DEALLOCATE(local_list)
855 END IF
856 END SUBROUTINE dmumps_build_mapping_info
857 SUBROUTINE dmumps_sol_omega(N, RHS,
858 & X, Y, R_W, C_W, IW, IFLAG,
859 & OMEGA, NOITER, TESTConv,
860 & LP, ARRET, GRAIN )
861 IMPLICIT NONE
862 INTEGER N, IFLAG
863 INTEGER IW(N,2)
864 DOUBLE PRECISION RHS(N)
865 DOUBLE PRECISION X(N), Y(N)
866 DOUBLE PRECISION R_W(N,2)
867 DOUBLE PRECISION C_W(N)
868 INTEGER LP, NOITER
869 LOGICAL TESTConv
870 DOUBLE PRECISION OMEGA(2)
871 DOUBLE PRECISION ARRET
872 INTEGER, intent(in) :: GRAIN
873 DOUBLE PRECISION, PARAMETER :: CGCE=0.2d0
874 DOUBLE PRECISION, PARAMETER :: CTAU=1.0d3
875 INTEGER I, IMAX
876 DOUBLE PRECISION OM1, OM2, DXMAX
877 DOUBLE PRECISION TAU, DD
878 DOUBLE PRECISION OLDOMG(2)
879 DOUBLE PRECISION, PARAMETER :: ZERO=0.0d0
880 DOUBLE PRECISION, PARAMETER :: ONE=1.0d0
881 INTEGER DMUMPS_IXAMAX
882 INTRINSIC abs, max
883 SAVE om1, oldomg
884 imax = dmumps_ixamax(n, x, 1, grain)
885 dxmax = abs(x(imax))
886 omega(1) = zero
887 omega(2) = zero
888 DO i = 1, n
889 tau = (r_w(i, 2) * dxmax + abs(rhs(i))) * dble(n) * ctau
890 dd = r_w(i, 1) + abs(rhs(i))
891 IF (dd .GT. tau * epsilon(ctau)) THEN
892 omega(1) = max(omega(1), abs(y(i)) / dd)
893 iw(i, 1) = 1
894 ELSE
895 IF (tau .GT. zero) THEN
896 omega(2) = max(omega(2),
897 & abs(y(i)) / (dd + r_w(i, 2) * dxmax))
898 ENDIF
899 iw(i, 1) = 2
900 ENDIF
901 ENDDO
902 IF (testconv) THEN
903 om2 = omega(1) + omega(2)
904 IF (om2 .LT. arret ) THEN
905 iflag = 1
906 GOTO 70
907 ENDIF
908 IF (noiter .GE. 1) THEN
909 IF (om2 .GT. om1 * cgce) THEN
910 IF (om2 .GT. om1) THEN
911 omega(1) = oldomg(1)
912 omega(2) = oldomg(2)
913 DO i = 1, n
914 x(i) = c_w(i)
915 ENDDO
916 iflag = 2
917 GOTO 70
918 ENDIF
919 iflag = 3
920 GOTO 70
921 ENDIF
922 ENDIF
923 DO i = 1, n
924 c_w(i) = x(i)
925 ENDDO
926 oldomg(1) = omega(1)
927 oldomg(2) = omega(2)
928 om1 = om2
929 ENDIF
930 iflag = 0
931 RETURN
932 70 CONTINUE
933 RETURN
934 END SUBROUTINE dmumps_sol_omega
935 SUBROUTINE dmumps_sol_lcond(N, RHS,
936 & X, Y, D, R_W, C_W, IW, KASE,
937 & OMEGA, ERX, COND,
938 & LP, KEEP,KEEP8 )
939 IMPLICIT NONE
940 INTEGER N, KASE, KEEP(500)
941 INTEGER(8) KEEP8(150)
942 INTEGER IW(N,2)
943 DOUBLE PRECISION RHS(N)
944 DOUBLE PRECISION X(N), Y(N)
945 DOUBLE PRECISION D(N)
946 DOUBLE PRECISION R_W(N,2)
947 DOUBLE PRECISION C_W(N)
948 INTEGER LP
949 DOUBLE PRECISION COND(2),OMEGA(2)
950 LOGICAL LCOND1, LCOND2
951 INTEGER JUMP, I, IMAX
952 DOUBLE PRECISION ERX, DXMAX
953 DOUBLE PRECISION DXIMAX
954 DOUBLE PRECISION, PARAMETER :: ZERO = 0.0d0
955 DOUBLE PRECISION, PARAMETER :: ONE = 1.0d0
956 INTEGER DMUMPS_IXAMAX
957 INTRINSIC abs, max
958 SAVE lcond1, lcond2, jump, dximax, dxmax
959 IF (kase .EQ. 0) THEN
960 lcond1 = .false.
961 lcond2 = .false.
962 cond(1) = one
963 cond(2) = one
964 erx = zero
965 jump = 1
966 ENDIF
967 SELECT CASE (jump)
968 CASE (1)
969 GOTO 30
970 CASE(2)
971 GOTO 10
972 CASE(3)
973 GOTO 110
974 CASE(4)
975 GOTO 150
976 CASE(5)
977 GOTO 35
978 CASE DEFAULT
979 END SELECT
980 10 CONTINUE
981 30 CONTINUE
982 35 CONTINUE
983 imax = dmumps_ixamax(n, x, 1, keep(361))
984 dxmax = abs(x(imax))
985 DO i = 1, n
986 IF (iw(i, 1) .EQ. 1) THEN
987 r_w(i, 1) = r_w(i, 1) + abs(rhs(i))
988 r_w(i, 2) = zero
989 lcond1 = .true.
990 ELSE
991 r_w(i, 2) = r_w(i, 2) * dxmax + r_w(i, 1)
992 r_w(i, 1) = zero
993 lcond2 = .true.
994 ENDIF
995 ENDDO
996 DO i = 1, n
997 c_w(i) = x(i) * d(i)
998 ENDDO
999 imax = dmumps_ixamax(n, c_w(1), 1, keep(361))
1000 dximax = abs(c_w(imax))
1001 IF (.NOT.lcond1) GOTO 130
1002 100 CONTINUE
1003 CALL dmumps_sol_b(n, kase, y, cond(1), c_w, iw(1, 2), keep(361))
1004 IF (kase .EQ. 0) GOTO 120
1005 IF (kase .EQ. 1) CALL dmumps_sol_mulr(n, y, d)
1006 IF (kase .EQ. 2) CALL dmumps_sol_mulr(n, y, r_w)
1007 jump = 3
1008 RETURN
1009 110 CONTINUE
1010 IF (kase .EQ. 1) CALL dmumps_sol_mulr(n, y, r_w)
1011 IF (kase .EQ. 2) CALL dmumps_sol_mulr(n, y, d)
1012 GOTO 100
1013 120 CONTINUE
1014 IF (dximax .GT. zero) cond(1) = cond(1) / dximax
1015 erx = omega(1) * cond(1)
1016 130 CONTINUE
1017 IF (.NOT.lcond2) GOTO 170
1018 kase = 0
1019 140 CONTINUE
1020 CALL dmumps_sol_b(n, kase, y, cond(2), c_w, iw(1, 2), keep(361))
1021 IF (kase .EQ. 0) GOTO 160
1022 IF (kase .EQ. 1) CALL dmumps_sol_mulr(n, y, d)
1023 IF (kase .EQ. 2) CALL dmumps_sol_mulr(n, y, r_w(1, 2))
1024 jump = 4
1025 RETURN
1026 150 CONTINUE
1027 IF (kase .EQ. 1) CALL dmumps_sol_mulr(n, y, r_w(1, 2))
1028 IF (kase .EQ. 2) CALL dmumps_sol_mulr(n, y, d)
1029 GOTO 140
1030 160 IF (dximax .GT. zero) THEN
1031 cond(2) = cond(2) / dximax
1032 ENDIF
1033 erx = erx + omega(2) * cond(2)
1034 170 CONTINUE
1035 RETURN
1036 END SUBROUTINE dmumps_sol_lcond
1037 SUBROUTINE dmumps_sol_cpy_fs2rhscomp( JBDEB, JBFIN, NBROWS,
1038 & KEEP, RHSCOMP, NRHS, LRHSCOMP, FIRST_ROW_RHSCOMP, W, LD_W,
1039 & FIRST_ROW_W )
1040 INTEGER :: JBDEB, JBFIN, NBROWS
1041 INTEGER :: NRHS, LRHSCOMP
1042 INTEGER :: FIRST_ROW_RHSCOMP
1043 INTEGER, INTENT(IN) :: KEEP(500)
1044 DOUBLE PRECISION, INTENT(INOUT) :: RHSCOMP(LRHSCOMP,NRHS)
1045 INTEGER :: LD_W, FIRST_ROW_W
1046 DOUBLE PRECISION :: W(LD_W*(JBFIN-JBDEB+1))
1047 INTEGER :: JJ, K, ISHIFT
1048!$OMP PARALLEL DO PRIVATE(ISHIFT, JJ), IF
1049!$OMP& (JBFIN-JBDEB+1 > 2*KEEP(362) .AND.
1050!$OMP& NBROWS * (JBFIN-JBDEB+1) > 2*KEEP(363))
1051 DO k = jbdeb, jbfin
1052 ishift = first_row_w + ld_w * (k-jbdeb)
1053 DO jj = 0, nbrows-1
1054 rhscomp(first_row_rhscomp+jj,k) = w(ishift+jj)
1055 END DO
1056 END DO
1057!$OMP END PARALLEL DO
1058 RETURN
1059 END SUBROUTINE dmumps_sol_cpy_fs2rhscomp
1060 SUBROUTINE dmumps_sol_bwd_gthr( JBDEB, JBFIN, J1, J2,
1061 & RHSCOMP, NRHS, LRHSCOMP, W, LD_W, FIRST_ROW_W,
1062 & IW, LIW, KEEP, N, POSINRHSCOMP_BWD )
1063 INTEGER, INTENT(IN) :: JBDEB, JBFIN, J1, J2
1064 INTEGER, INTENT(IN) :: NRHS, LRHSCOMP
1065 INTEGER, INTENT(IN) :: FIRST_ROW_W, LD_W, LIW
1066 INTEGER, INTENT(IN) :: IW(LIW)
1067 INTEGER, INTENT(IN) :: KEEP(500)
1068 DOUBLE PRECISION, INTENT(INOUT) :: RHSCOMP(LRHSCOMP,NRHS)
1069 DOUBLE PRECISION :: W(LD_W*(JBFIN-JBDEB+1))
1070 INTEGER, INTENT(IN) :: N
1071 INTEGER, INTENT(IN) :: POSINRHSCOMP_BWD(N)
1072 INTEGER :: ISHIFT, JJ, K, IPOSINRHSCOMP
1073!$OMP PARALLEL DO PRIVATE(JJ,ISHIFT,IPOSINRHSCOMP), IF
1074!$OMP& ((JBFIN-JBDEB+1 > 2*KEEP(362) .AND.
1075!$OMP& (JBFIN-JBDEB+1)*(J2-KEEP(253)-J1+1)>2*KEEP(363)))
1076 DO k=jbdeb, jbfin
1077 ishift = first_row_w+(k-jbdeb)*ld_w
1078 DO jj = j1, j2-keep(253)
1079 iposinrhscomp = abs(posinrhscomp_bwd(iw(jj)))
1080 w(ishift+jj-j1)= rhscomp(iposinrhscomp,k)
1081 ENDDO
1082 ENDDO
1083!$OMP END PARALLEL DO
1084 RETURN
1085 END SUBROUTINE dmumps_sol_bwd_gthr
1086 SUBROUTINE dmumps_sol_q(MTYPE, IFLAG, N,
1087 & LHS, WRHS, W, RES, GIVNORM, ANORM, XNORM, SCLNRM,
1088 & MPRINT, ICNTL, KEEP,KEEP8)
1089 INTEGER MTYPE,N,IFLAG,ICNTL(60), KEEP(500)
1090 INTEGER(8) KEEP8(150)
1091 DOUBLE PRECISION RES(N),LHS(N)
1092 DOUBLE PRECISION WRHS(N)
1093 DOUBLE PRECISION W(N)
1094 DOUBLE PRECISION RESMAX,RESL2,XNORM, SCLNRM
1095 DOUBLE PRECISION ANORM,DZERO
1096 LOGICAL GIVNORM,PROK
1097 INTEGER MPRINT, MP
1098 INTEGER K
1099 INTRINSIC abs, max, sqrt
1100 mp = icntl(2)
1101 prok = (mprint .GT. 0)
1102 dzero = 0.0d0
1103 IF (.NOT.givnorm) anorm = dzero
1104 resmax = dzero
1105 resl2 = dzero
1106 DO 40 k = 1, n
1107 resmax = max(resmax, abs(res(k)))
1108 resl2 = resl2 + abs(res(k)) * abs(res(k))
1109 IF (.NOT.givnorm) anorm = max(anorm, w(k))
1110 40 CONTINUE
1111 xnorm = dzero
1112 DO 50 k = 1, n
1113 xnorm = max(xnorm, abs(lhs(k)))
1114 50 CONTINUE
1115 IF ( xnorm .EQ. dzero .OR. (exponent(xnorm) .LT.
1116 & minexponent(xnorm) + keep(122) )
1117 & .OR.
1118 & ( exponent(anorm)+exponent(xnorm) .LT.
1119 & minexponent(xnorm) + keep(122) )
1120 & .OR.
1121 & ( exponent(anorm) + exponent(xnorm) -exponent(resmax)
1122 & .LT. minexponent(xnorm) + keep(122) )
1123 & ) THEN
1124 IF (mod(iflag/2,2) .EQ. 0) THEN
1125 iflag = iflag + 2
1126 ENDIF
1127 IF ((mp .GT. 0) .AND. (icntl(4) .GE. 2)) WRITE( mp, * )
1128 & ' max-NORM of computed solut. is zero or close to zero. '
1129 ENDIF
1130 IF (resmax .EQ. dzero) THEN
1131 sclnrm = dzero
1132 ELSE
1133 sclnrm = resmax / (anorm * xnorm)
1134 ENDIF
1135 resl2 = sqrt(resl2)
1136 IF (prok) WRITE( mprint, 90 ) resmax, resl2, anorm, xnorm,
1137 & sclnrm
1138 90 FORMAT (/' RESIDUAL IS ............ (MAX-NORM) =',1pd9.2/
1139 & ' .. (2-NORM) =',1pd9.2/
1140 & ' RINFOG(4):NORM OF input Matrix (MAX-NORM)=',1pd9.2/
1141 & ' RINFOG(5):NORM OF Computed SOLUT (MAX-NORM)=',1pd9.2/
1142 & ' RINFOG(6):SCALED RESIDUAL ...... (MAX-NORM)=',1pd9.2)
1143 RETURN
1144 END SUBROUTINE dmumps_sol_q
1145 SUBROUTINE dmumps_solve_fwd_trsolve (A, LA, APOS, NPIV, LDADIAG,
1146 & NRHS_B, WCB, LWCB, LDA_WCB, PPIV_COURANT, MTYPE, KEEP)
1147 INTEGER, INTENT(IN) :: MTYPE, LDADIAG, NPIV, KEEP(500)
1148 INTEGER, INTENT(IN) :: NRHS_B, LDA_WCB
1149 INTEGER(8), INTENT(IN) :: LA, APOS, LWCB, PPIV_COURANT
1150 DOUBLE PRECISION, INTENT(IN) :: A(LA)
1151 DOUBLE PRECISION, INTENT(INOUT) :: WCB(LWCB)
1152 DOUBLE PRECISION ONE
1153 PARAMETER (ONE = 1.0d0)
1154 IF (keep(50).NE.0 .OR. mtype .eq. 1 ) THEN
1155#if defined(MUMPS_USE_BLAS2)
1156 IF ( nrhs_b == 1 ) THEN
1157 CALL dtrsv( 'U', 'T', 'U', npiv, a(apos), ldadiag,
1158 & wcb(ppiv_courant), 1 )
1159 ELSE
1160#endif
1161 CALL dtrsm( 'L','U','T','U', npiv, nrhs_b, one,
1162 & a(apos), ldadiag, wcb(ppiv_courant),
1163 & lda_wcb )
1164#if defined(MUMPS_USE_BLAS2)
1165 ENDIF
1166#endif
1167 ELSE
1168#if defined(MUMPS_USE_BLAS2)
1169 IF ( nrhs_b == 1 ) THEN
1170 CALL dtrsv( 'L', 'N', 'N', npiv, a(apos), ldadiag,
1171 & wcb(ppiv_courant), 1 )
1172 ELSE
1173#endif
1174 CALL dtrsm( 'L','L','N','N', npiv, nrhs_b, one,
1175 & a(apos), ldadiag, wcb(ppiv_courant),
1176 & lda_wcb )
1177#if defined(MUMPS_USE_BLAS2)
1178 ENDIF
1179#endif
1180 ENDIF
1181 RETURN
1182 END SUBROUTINE dmumps_solve_fwd_trsolve
1183 SUBROUTINE dmumps_solve_bwd_trsolve (A, LA, APOS, NPIV, LDADIAG,
1184 & NRHS_B, WCB, LWCB, LDA_WCB, PPIV_COURANT, MTYPE, KEEP)
1185 INTEGER, INTENT(IN) :: MTYPE, LDADIAG, NPIV, KEEP(500)
1186 INTEGER, INTENT(IN) :: NRHS_B, LDA_WCB
1187 INTEGER(8), INTENT(IN) :: LA, APOS, LWCB, PPIV_COURANT
1188 DOUBLE PRECISION, INTENT(IN) :: A(LA)
1189 DOUBLE PRECISION, INTENT(INOUT) :: WCB(LWCB)
1190 DOUBLE PRECISION ONE
1191 PARAMETER (ONE = 1.0d0)
1192 IF (mtype .eq. 1 ) THEN
1193#if defined(MUMPS_USE_BLAS2)
1194 IF ( nrhs_b == 1 ) THEN
1195 CALL dtrsv( 'L', 'T', 'N', npiv, a(apos), ldadiag,
1196 & wcb(ppiv_courant), 1 )
1197 ELSE
1198#endif
1199 CALL dtrsm( 'L','L','T','N', npiv, nrhs_b, one,
1200 & a(apos), ldadiag, wcb(ppiv_courant),
1201 & lda_wcb )
1202#if defined(MUMPS_USE_BLAS2)
1203 ENDIF
1204#endif
1205 ELSE
1206#if defined(MUMPS_USE_BLAS2)
1207 IF ( nrhs_b == 1 ) THEN
1208 CALL dtrsv( 'U', 'N', 'U', npiv, a(apos), ldadiag,
1209 & wcb(ppiv_courant), 1 )
1210 ELSE
1211#endif
1212 CALL dtrsm( 'L','U','N','U', npiv, nrhs_b, one,
1213 & a(apos), ldadiag, wcb(ppiv_courant),
1214 & lda_wcb )
1215#if defined(MUMPS_USE_BLAS2)
1216 ENDIF
1217#endif
1218 ENDIF
1219 RETURN
1220 END SUBROUTINE dmumps_solve_bwd_trsolve
1221 SUBROUTINE dmumps_solve_fwd_panels(
1222 & A, LA, APOS, NPIV, IW,
1223 & NRHS_B, WCB, LWCB, LDA_WCB,
1224 & PPIV_COURANT, MTYPE, KEEP)
1225 INTEGER, INTENT(IN) :: MTYPE, NPIV, KEEP(500)
1226 INTEGER, INTENT(IN) :: IW(NPIV)
1227 INTEGER, INTENT(IN) :: NRHS_B, LDA_WCB
1228 INTEGER(8), INTENT(IN) :: LA, APOS, LWCB, PPIV_COURANT
1229 DOUBLE PRECISION, INTENT(IN) :: A(LA)
1230 DOUBLE PRECISION, INTENT(INOUT) :: WCB(LWCB)
1231 INTEGER :: NB_TARGET
1232 INTEGER :: NBPANELS
1233 INTEGER :: NBROWS_PANEL, NBCOLS_PANEL, ICOL_BEG, ICOL_END
1234 INTEGER(8) :: PANEL_APOS, PPIV_PANEL
1235 DOUBLE PRECISION, PARAMETER :: ONE=1.0d0
1236 IF (keep(459) .LE. 1) THEN
1237 WRITE(*,*) " Internal error in DMUMPS_SOLVE_FWD_PANELS"
1238 CALL mumps_abort()
1239 ENDIF
1240 CALL mumps_ldltpanel_nbtarget( npiv, nb_target, keep )
1241 panel_apos = apos
1242 nbpanels = 0
1243 icol_beg = 1
1244 nbrows_panel = npiv
1245 ppiv_panel = ppiv_courant
1246 DO WHILE ( icol_beg .LE. npiv )
1247 nbpanels = nbpanels + 1
1248 icol_end = min(nb_target * nbpanels, npiv)
1249 IF ( iw(icol_end) .LT. 0 ) icol_end=icol_end+1
1250 nbcols_panel = icol_end - icol_beg + 1
1251 CALL dmumps_solve_fwd_trsolve (a, la, panel_apos,
1252 & nbcols_panel, nbcols_panel,
1253 & nrhs_b, wcb, lwcb, lda_wcb, ppiv_panel, mtype, keep)
1254 IF ( nbrows_panel .GT. nbcols_panel ) THEN
1255 CALL dmumps_solve_gemm_update( a, la,
1256 & panel_apos + int(nbcols_panel,8) * int(nbcols_panel,8),
1257 & nbcols_panel, nbcols_panel, nbrows_panel-nbcols_panel,
1258 & nrhs_b, wcb, lwcb, ppiv_panel, lda_wcb,
1259 & ppiv_panel+nbcols_panel, lda_wcb,
1260 & mtype, keep, one )
1261 ENDIF
1262 icol_beg = icol_end + 1
1263 panel_apos = panel_apos + int(nbcols_panel,8) *
1264 & int(nbrows_panel,8)
1265 nbrows_panel = nbrows_panel - nbcols_panel
1266 ppiv_panel = ppiv_panel + nbcols_panel
1267 ENDDO
1268 RETURN
1269 END SUBROUTINE dmumps_solve_fwd_panels
1270 SUBROUTINE dmumps_solve_bwd_panels(
1271 & A, LA, APOS, NPIV, IW,
1272 & NRHS_B, WCB, LWCB, LDA_WCB,
1273 & PPIV_COURANT, MTYPE, KEEP)
1274 INTEGER, INTENT(IN) :: MTYPE, NPIV, KEEP(500)
1275 INTEGER, INTENT(IN) :: IW(NPIV)
1276 INTEGER, INTENT(IN) :: NRHS_B, LDA_WCB
1277 INTEGER(8), INTENT(IN) :: LA, APOS, LWCB, PPIV_COURANT
1278 DOUBLE PRECISION, INTENT(IN) :: A(LA)
1279 DOUBLE PRECISION, INTENT(INOUT) :: WCB(LWCB)
1280 INTEGER, PARAMETER :: PANEL_TABSIZE = 20
1281 INTEGER(8) :: PANEL_POS(PANEL_TABSIZE)
1282 INTEGER :: PANEL_COL(PANEL_TABSIZE)
1283 INTEGER :: IPANEL, NBPANELS, NB_TARGET
1284 INTEGER :: NBROWS_PANEL, NBCOLS_PANEL
1285 INTEGER(8) :: PPIV_PANEL
1286 INTEGER :: MTYPE_TEMP
1287 DOUBLE PRECISION, PARAMETER :: ONE=1.0d0
1288 IF (keep(459) .LE. 1) THEN
1289 WRITE(*,*) " Internal error 1 in DMUMPS_SOLVE_BWD_PANELS"
1290 CALL mumps_abort()
1291 ENDIF
1292 IF ( keep(459)+1 .GT. panel_tabsize ) THEN
1293 WRITE(*,*) " Internal error 2 in DMUMPS_SOLVE_BWD_PANELS"
1294 CALL mumps_abort()
1295 ENDIF
1296 CALL mumps_ldltpanel_panelinfos( npiv, keep, iw,
1297 &nb_target, nbpanels, panel_col, panel_pos, panel_tabsize,
1298 &.false. )
1299 DO ipanel = nbpanels, 1, -1
1300 nbcols_panel = panel_col( ipanel+1 ) - panel_col( ipanel )
1301 nbrows_panel = npiv - panel_col( ipanel ) + 1
1302 ppiv_panel = ppiv_courant + panel_col( ipanel ) - 1
1303 IF ( nbrows_panel .GT. nbcols_panel ) THEN
1304 mtype_temp = 0
1305 CALL dmumps_solve_gemm_update( a, la,
1306 & apos-1_8+panel_pos(ipanel)+
1307 & int(nbcols_panel,8)*int(nbcols_panel,8),
1308 & nbrows_panel-nbcols_panel, nbcols_panel,
1309 & nbcols_panel,
1310 & nrhs_b, wcb, lwcb, ppiv_panel+nbcols_panel, lda_wcb,
1311 & ppiv_panel, lda_wcb,
1312 & mtype_temp, keep, one )
1313 ENDIF
1314 CALL dmumps_solve_bwd_trsolve (a, la,
1315 & apos+panel_pos(ipanel)-1_8,
1316 & nbcols_panel, nbcols_panel,
1317 & nrhs_b, wcb, lwcb, lda_wcb, ppiv_panel, mtype, keep)
1318 ENDDO
1319 RETURN
1320 END SUBROUTINE dmumps_solve_bwd_panels
1321 SUBROUTINE dmumps_solve_gemm_update
1322 & (a, la, apos1, nx, lda, ny,
1323 & nrhs_b, wcb, lwcb, ptrx, ldx,
1324 & ptry, ldy,
1325 & mtype, keep, coef_y )
1326 INTEGER, INTENT(IN) :: MTYPE, NY, NX, KEEP(500)
1327 INTEGER, INTENT(IN) :: NRHS_B, LDY, LDA, LDX
1328 INTEGER(8), INTENT(IN) :: LA, APOS1, LWCB, PTRX,
1329 & ptry
1330 DOUBLE PRECISION, INTENT(IN) :: A(LA)
1331 DOUBLE PRECISION, INTENT(INOUT) :: WCB(LWCB)
1332 DOUBLE PRECISION, INTENT(IN) :: COEF_Y
1333 DOUBLE PRECISION ALPHA, ZERO, ONE
1334 PARAMETER (ZERO = 0.0d0, one = 1.0d0, alpha=-1.0d0)
1335 IF ( nx .NE. 0 .AND. ny.NE.0 ) THEN
1336 IF ( mtype .eq. 1 ) THEN
1337#if defined(MUMPS_USE_BLAS2)
1338 IF ( nrhs_b == 1 ) THEN
1339 CALL dgemv('T', nx, ny, alpha, a(apos1),
1340 & lda, wcb(ptrx), 1, coef_y,
1341 & wcb(ptry), 1)
1342 ELSE
1343#endif
1344 CALL dgemm('T', 'N', ny, nrhs_b, nx, alpha,
1345 & a(apos1), lda, wcb(ptrx), ldx, coef_y,
1346 & wcb(ptry), ldy)
1347#if defined(MUMPS_USE_BLAS2)
1348 END IF
1349#endif
1350 ELSE
1351#if defined(MUMPS_USE_BLAS2)
1352 IF ( nrhs_b == 1 ) THEN
1353 CALL dgemv('N',ny, nx, alpha, a(apos1),
1354 & lda, wcb(ptrx), 1,
1355 & coef_y, wcb(ptry), 1 )
1356 ELSE
1357#endif
1358 CALL dgemm('N', 'n', NY, NRHS_B, NX, ALPHA,
1359 & A(APOS1), LDA, WCB(PTRX), LDX,
1360 & COEF_Y, WCB(PTRY), LDY)
1361#if defined(MUMPS_USE_BLAS2)
1362 END IF
1363#endif
1364 END IF
1365 END IF
1366 RETURN
1367 END SUBROUTINE DMUMPS_SOLVE_GEMM_UPDATE
1368 SUBROUTINE DMUMPS_SOL_LD_AND_RELOAD_PANEL (
1369 & INODE, N, NPIV, LIELL, NELIM, NSLAVES,
1370 & PPIV_COURANT,
1371 & IW, IPOS, LIW,
1372 & A, LA, APOS,
1373 & WCB, LWCB, LD_WCBPIV,
1374 & RHSCOMP, LRHSCOMP, NRHS,
1375 & POSINRHSCOMP_FWD, JBDEB, JBFIN,
1376 & MTYPE, KEEP, OOCWRITE_COMPATIBLE_WITH_BLR,
1377 & IGNORE_K459
1378 & )
1379 USE DMUMPS_OOC
1380 IMPLICIT NONE
1381 INTEGER, INTENT(IN) :: MTYPE, INODE, N, NPIV, LIELL,
1382 & NELIM, NSLAVES
1383 INTEGER, INTENT(IN) :: LRHSCOMP, NRHS, LIW, JBDEB, JBFIN
1384 INTEGER, INTENT(IN) :: IW(LIW), IPOS, POSINRHSCOMP_FWD(N)
1385 INTEGER(8), INTENT(IN) :: LWCB, APOS, LA, PPIV_COURANT
1386 INTEGER, INTENT(IN) :: LD_WCBPIV
1387 INTEGER, INTENT(IN) :: KEEP(500)
1388 DOUBLE PRECISION, INTENT(IN) :: WCB( LWCB )
1389 DOUBLE PRECISION, INTENT(IN) :: A( LA )
1390 DOUBLE PRECISION, INTENT(INOUT) :: RHSCOMP(LRHSCOMP, NRHS)
1391 LOGICAL, INTENT(IN) :: OOCWRITE_COMPATIBLE_WITH_BLR
1392 LOGICAL, INTENT(IN) :: IGNORE_K459
1393 INTEGER :: J1, J3
1394 INTEGER :: IPOSINRHSCOMP, JJ, K, NBK,
1395 & LDAJ, NRHS_B
1396 INTEGER(8) :: IFR8 , APOS1, APOS2, APOSOFF, IFR_ini8,
1397 & POSWCB1, POSWCB2
1398 DOUBLE PRECISION :: VALPIV, A11, A22, A12, DETPIV
1399 INTEGER, PARAMETER :: PANEL_TABSIZE = 20
1400 INTEGER(8) :: PANEL_POS(PANEL_TABSIZE)
1401 INTEGER :: PANEL_COL(PANEL_TABSIZE)
1402 INTEGER :: IPANEL, ICOL, NBPANELS, NB_TARGET
1403 LOGICAL :: SKIP_IT
1404 LOGICAL :: OMP_FLAG
1405 DOUBLE PRECISION ONE
1406 PARAMETER (ONE = 1.0D0)
1407.EQ. IF ( NPIV 0 ) RETURN
1408 NRHS_B = JBFIN-JBDEB+1
1409.EQ..OR..NE. IF ( MTYPE 1 KEEP(50) 0 ) THEN
1410 J1 = IPOS + 1
1411 J3 = IPOS + NPIV
1412 ELSE
1413 J1 = IPOS + LIELL + 1
1414 J3 = IPOS + LIELL + NPIV
1415 END IF
1416 IPOSINRHSCOMP = POSINRHSCOMP_FWD(IW(J1))
1417.eq. IF ( KEEP(50) 0 ) THEN
1418 OMP_FLAG = .FALSE.
1419.GE.!$ OMP_FLAG=(int(NRHS_B,8)*int(NPIV,8)int(KEEP(363),8))
1420 IF (OMP_FLAG) THEN
1421!$OMP PARALLEL DO PRIVATE(IFR8) COLLAPSE(2)
1422 DO K = JBDEB, JBFIN
1423 DO IFR8 = 0_8, int(NPIV-1,8)
1424 RHSCOMP(IPOSINRHSCOMP+IFR8, K) =
1425 & WCB(PPIV_COURANT+(K-JBDEB)*LD_WCBPIV+IFR8)
1426 ENDDO
1427 ENDDO
1428!$OMP END PARALLEL DO
1429 ELSE
1430 DO K = JBDEB, JBFIN
1431 DO IFR8 = 0_8, int(NPIV-1,8)
1432 RHSCOMP(IPOSINRHSCOMP+IFR8, K) =
1433 & WCB(PPIV_COURANT+(K-JBDEB)*LD_WCBPIV+IFR8)
1434 ENDDO
1435 ENDDO
1436 ENDIF
1437 ELSE
1438 CALL MUMPS_LDLTPANEL_PANELINFOS( NPIV, KEEP, IW(IPOS+LIELL+1),
1439 & NB_TARGET, NBPANELS, PANEL_COL, PANEL_POS, PANEL_TABSIZE,
1440 & IGNORE_K459 )
1441 IFR_ini8 = PPIV_COURANT
1442.GE..AND.!$ OMP_FLAG = ( JBFIN-JBDEB+1KEEP(362)
1443.GE.!$ & ((J3-J1+1)*(JBFIN-JBDEB+1) KEEP(363)))
1444!$OMP PARALLEL DO PRIVATE(JJ,IFR8,NBK,APOS1,APOS2,APOSOFF,VALPIV,
1445!$OMP& IPANEL,ICOL,
1446!$OMP& POSWCB1,POSWCB2,A11,A22,A12,DETPIV,LDAJ,SKIP_IT)
1447!$OMP& IF(OMP_FLAG)
1448 DO K = JBDEB, JBFIN
1449 DO JJ = J1, J3
1450 IPANEL = (JJ-J1)/NB_TARGET + 1
1451.LT. IF ( JJ-J1+1 PANEL_COL(IPANEL) ) IPANEL = IPANEL -1
1452 ICOL = JJ-J1+1 - PANEL_COL(IPANEL) + 1
1453 LDAJ = PANEL_COL(IPANEL+1) - PANEL_COL(IPANEL)
1454 APOS1 = APOS-1_8+PANEL_POS( IPANEL ) + int(ICOL-1,8) *
1455 & int(LDAJ+1,8)
1456 IFR8 = IFR_ini8 + int(K-JBDEB,8)*int(LD_WCBPIV,8) +
1457 & int(JJ-J1,8)
1458.NE. IF ( JJ J1 ) THEN
1459.LT. IF ( IW(LIELL+JJ-1) 0 ) THEN
1460 SKIP_IT = .TRUE.
1461 ELSE
1462 SKIP_IT = .FALSE.
1463 ENDIF
1464 ELSE
1465 SKIP_IT = .FALSE.
1466 ENDIF
1467 IF (SKIP_IT) THEN
1468.GT. ELSE IF ( IW(JJ+LIELL) 0 ) THEN
1469 VALPIV = ONE/A( APOS1 )
1470 RHSCOMP(IPOSINRHSCOMP+JJ-J1 , K ) =
1471 & WCB( IFR8 ) * VALPIV
1472 APOS1 = APOS1 + int(LDAJ + 1,8)
1473 ELSE
1474 APOS2 = APOS1+int(LDAJ+1,8)
1475 APOSOFF=APOS1+1_8
1476 A11 = A(APOS1)
1477 A22 = A(APOS2)
1478 A12 = A(APOSOFF)
1479 DETPIV = A11*A22 - A12**2
1480 A22 = A11/DETPIV
1481 A11 = A(APOS2)/DETPIV
1482 A12 = -A12/DETPIV
1483 POSWCB1 = IFR8
1484 POSWCB2 = POSWCB1+1_8
1485 RHSCOMP(IPOSINRHSCOMP+JJ-J1,K) =
1486 & WCB(POSWCB1)*A11
1487 & + WCB(POSWCB2)*A12
1488 RHSCOMP(IPOSINRHSCOMP+JJ-J1+1,K) =
1489 & WCB(POSWCB1)*A12
1490 & + WCB(POSWCB2)*A22
1491 ENDIF
1492 ENDDO
1493 ENDDO
1494!$OMP END PARALLEL DO
1495 END IF
1496 RETURN
1497 END SUBROUTINE DMUMPS_SOL_LD_AND_RELOAD_PANEL
1498 SUBROUTINE DMUMPS_SOL_LD_AND_RELOAD (
1499 & INODE, N, NPIV, LIELL, NELIM, NSLAVES,
1500 & PPIV_COURANT,
1501 & IW, IPOS, LIW,
1502 & A, LA, APOS,
1503 & WCB, LWCB, LD_WCBPIV,
1504 & RHSCOMP, LRHSCOMP, NRHS,
1505 & POSINRHSCOMP_FWD, JBDEB, JBFIN,
1506 & MTYPE, KEEP, OOCWRITE_COMPATIBLE_WITH_BLR,
1507 & IGNORE_K459
1508 & )
1509 USE DMUMPS_OOC
1510 INTEGER, INTENT(IN) :: MTYPE, INODE, N, NPIV, LIELL,
1511 & NELIM, NSLAVES
1512 INTEGER, INTENT(IN) :: LRHSCOMP, NRHS, LIW, JBDEB, JBFIN
1513 INTEGER, INTENT(IN) :: IW(LIW), IPOS, POSINRHSCOMP_FWD(N)
1514 INTEGER(8), INTENT(IN) :: LWCB, APOS, LA, PPIV_COURANT
1515 INTEGER, INTENT(IN) :: LD_WCBPIV
1516 INTEGER, INTENT(IN) :: KEEP(500)
1517 DOUBLE PRECISION, INTENT(IN) :: WCB( LWCB )
1518 DOUBLE PRECISION, INTENT(IN) :: A( LA )
1519 DOUBLE PRECISION, INTENT(INOUT) :: RHSCOMP(LRHSCOMP, NRHS)
1520 LOGICAL, INTENT(IN) :: OOCWRITE_COMPATIBLE_WITH_BLR
1521 LOGICAL, INTENT(IN) :: IGNORE_K459
1522 INTEGER :: TempNROW, J1, J3, PANEL_SIZE, TYPEF
1523 INTEGER :: IPOSINRHSCOMP, JJ, K, NBK, LDAJ,
1524 & LDAJ_ini, NBK_ini, LDAJ_FIRST_PANEL, NRHS_B
1525 INTEGER(8) :: IFR8 , APOS1, APOS2, APOSOFF, IFR_ini8,
1526 & POSWCB1, POSWCB2
1527 DOUBLE PRECISION :: VALPIV, A11, A22, A12, DETPIV
1528!$ LOGICAL :: OMP_FLAG
1529 DOUBLE PRECISION ONE
1530 PARAMETER (ONE = 1.0D0)
1531 NRHS_B = JBFIN-JBDEB+1
1532.EQ..OR..NE. IF ( MTYPE 1 KEEP(50) 0 ) THEN
1533 J1 = IPOS + 1
1534 J3 = IPOS + NPIV
1535 ELSE
1536 J1 = IPOS + LIELL + 1
1537 J3 = IPOS + LIELL + NPIV
1538 END IF
1539 IPOSINRHSCOMP = POSINRHSCOMP_FWD(IW(J1))
1540.eq. IF ( KEEP(50) 0 ) THEN
1541.GE..AND..GE.!$ OMP_FLAG=(NRHS_BKEEP(362)NRHS_B*NPIVKEEP(363))
1542!$OMP PARALLEL DO PRIVATE(IFR8) IF (OMP_FLAG)
1543 DO K=JBDEB,JBFIN
1544 IFR8 = PPIV_COURANT + (K-JBDEB)*LD_WCBPIV
1545 RHSCOMP(IPOSINRHSCOMP:IPOSINRHSCOMP+NPIV-1, K) =
1546 & WCB(IFR8:IFR8+int(NPIV-1,8))
1547 ENDDO
1548!$OMP END PARALLEL DO
1549 ELSE
1550 IFR8 = PPIV_COURANT - 1_8
1551.EQ..AND. IF (KEEP(201)1OOCWRITE_COMPATIBLE_WITH_BLR) THEN
1552.EQ. IF (MTYPE1) THEN
1553.EQ..AND..NE. IF ((MTYPE1)NSLAVES0) THEN
1554 TempNROW= NPIV+NELIM
1555 LDAJ_FIRST_PANEL=TempNROW
1556 ELSE
1557 TempNROW= LIELL
1558 LDAJ_FIRST_PANEL=TempNROW
1559 ENDIF
1560 TYPEF=TYPEF_L
1561 ELSE
1562 TempNROW= NPIV
1563 LDAJ_FIRST_PANEL=LIELL
1564 TYPEF= TYPEF_U
1565 ENDIF
1566 PANEL_SIZE = DMUMPS_OOC_PANEL_SIZE( LDAJ_FIRST_PANEL )
1567 LDAJ = TempNROW
1568 ELSE
1569.GT..AND..NE. IF ( KEEP(459) 1 KEEP(50) 0
1570.AND..NOT. & IGNORE_K459 ) THEN
1571 CALL MUMPS_LDLTPANEL_NBTARGET( NPIV, PANEL_SIZE, KEEP )
1572 LDAJ = PANEL_SIZE
1573 ELSE
1574 PANEL_SIZE = -1
1575 LDAJ = NPIV
1576 ENDIF
1577 ENDIF
1578.EQ..AND. IF (KEEP(201)1OOCWRITE_COMPATIBLE_WITH_BLR) THEN
1579 NBK = 0
1580 ENDIF
1581 IFR_ini8 = PPIV_COURANT - 1_8
1582 LDAJ_ini = LDAJ
1583.EQ..AND. IF (KEEP(201)1OOCWRITE_COMPATIBLE_WITH_BLR)
1584 & NBK_ini = NBK
1585.GE..AND.!$ OMP_FLAG = ( JBFIN-JBDEB+1KEEP(362)
1586.GE.!$ & ((J3-J1+1)*(JBFIN-JBDEB+1) KEEP(363)))
1587!$OMP PARALLEL DO PRIVATE(JJ,IFR8,NBK,APOS1,APOS2,APOSOFF,VALPIV,
1588!$OMP& POSWCB1, POSWCB2,A11,A22,A12,DETPIV,LDAJ) IF(OMP_FLAG)
1589 DO K = JBDEB, JBFIN
1590 IFR8 = IFR_ini8 + int(K-JBDEB,8)*int(LD_WCBPIV,8)
1591 NBK = NBK_ini
1592 APOS1 = APOS
1593 LDAJ = LDAJ_ini
1594 JJ = J1
1595 DO
1596.GT. IF (JJ J3) EXIT
1597 IFR8 = IFR8 + 1_8
1598.GT. IF (IW(JJ+LIELL) 0) THEN
1599 VALPIV = ONE/A( APOS1 )
1600 RHSCOMP(IPOSINRHSCOMP+JJ-J1 , K ) =
1601 & WCB( IFR8 ) * VALPIV
1602.EQ..AND. IF (KEEP(201)1OOCWRITE_COMPATIBLE_WITH_BLR)
1603 & THEN
1604 NBK = NBK+1
1605.EQ. IF (NBKPANEL_SIZE) THEN
1606 NBK = 0
1607 LDAJ = LDAJ - PANEL_SIZE
1608 ENDIF
1609 ENDIF
1610 APOS1 = APOS1 + int(LDAJ + 1,8)
1611 JJ = JJ+1
1612 ELSE
1613.EQ..AND. IF (KEEP(201)1OOCWRITE_COMPATIBLE_WITH_BLR)
1614 & THEN
1615 NBK = NBK+1
1616 ENDIF
1617 APOS2 = APOS1+int(LDAJ+1,8)
1618.EQ..AND. IF (KEEP(201)1OOCWRITE_COMPATIBLE_WITH_BLR)
1619 & THEN
1620 APOSOFF = APOS1+int(LDAJ,8)
1621 ELSE
1622 APOSOFF=APOS1+1_8
1623 ENDIF
1624 A11 = A(APOS1)
1625 A22 = A(APOS2)
1626 A12 = A(APOSOFF)
1627 DETPIV = A11*A22 - A12**2
1628 A22 = A11/DETPIV
1629 A11 = A(APOS2)/DETPIV
1630 A12 = -A12/DETPIV
1631 POSWCB1 = IFR8
1632 POSWCB2 = POSWCB1+1_8
1633 RHSCOMP(IPOSINRHSCOMP+JJ-J1,K) =
1634 & WCB(POSWCB1)*A11
1635 & + WCB(POSWCB2)*A12
1636 RHSCOMP(IPOSINRHSCOMP+JJ-J1+1,K) =
1637 & WCB(POSWCB1)*A12
1638 & + WCB(POSWCB2)*A22
1639.EQ..AND. IF (KEEP(201)1OOCWRITE_COMPATIBLE_WITH_BLR)
1640 & THEN
1641 NBK = NBK+1
1642.GE. IF (NBKPANEL_SIZE) THEN
1643 LDAJ = LDAJ - NBK
1644 NBK = 0
1645 ENDIF
1646 ENDIF
1647 APOS1 = APOS2 + int(LDAJ + 1,8)
1648 JJ = JJ+2
1649 IFR8 = IFR8+1_8
1650 ENDIF
1651 ENDDO
1652 ENDDO
1653!$OMP END PARALLEL DO
1654 END IF
1655 RETURN
1656 END SUBROUTINE DMUMPS_SOL_LD_AND_RELOAD
1657 SUBROUTINE DMUMPS_SET_SCALING_LOC( scaling_data, N, ILOC, LILOC,
1658 & COMM, MYID, I_AM_SLAVE, MASTER, NB_BYTES, NB_BYTES_MAX,
1659 & K16_8, LP, LPOK, ICNTL, INFO )
1660 IMPLICIT NONE
1661 type scaling_data_t
1662 SEQUENCE
1663 DOUBLE PRECISION, dimension(:), pointer :: SCALING
1664 DOUBLE PRECISION, dimension(:), pointer :: SCALING_LOC
1665 end type scaling_data_t
1666 type (scaling_data_t), INTENT(INOUT) :: scaling_data
1667 INTEGER, INTENT(IN) :: N, LILOC, COMM, MYID, MASTER, LP
1668 INTEGER, INTENT(IN) :: ILOC(LILOC)
1669 INTEGER(8), INTENT(INOUT) :: NB_BYTES, NB_BYTES_MAX
1670 INTEGER(8), INTENT(IN) :: K16_8
1671 LOGICAL, INTENT(IN) :: I_AM_SLAVE, LPOK
1672 INTEGER, INTENT(INOUT) :: INFO(80)
1673 INTEGER, INTENT(IN) :: ICNTL(60)
1674 DOUBLE PRECISION, POINTER, DIMENSION(:) :: SCALING
1675 INTEGER :: I, IERR_MPI, allocok
1676 INCLUDE 'mpif.h'
1677 NULLIFY(scaling_data%SCALING_LOC)
1678 IF (I_AM_SLAVE) THEN
1679 ALLOCATE(scaling_data%SCALING_LOC(max(1,LILOC)),
1680 & stat=allocok)
1681 IF (allocok > 0) THEN
1682 INFO(1)=-13
1683 INFO(2)=max(1,LILOC)
1684 GOTO 35
1685 ENDIF
1686 NB_BYTES = NB_BYTES + int(max(1,LILOC),8)*K16_8
1687 NB_BYTES_MAX = max(NB_BYTES_MAX,NB_BYTES)
1688 ENDIF
1689.NE. IF (MYID MASTER) THEN
1690 ALLOCATE(SCALING(N), stat=allocok)
1691 IF (allocok > 0) THEN
1692 IF (LPOK) THEN
1693 WRITE(LP,*) 'error allocating temporary scaling array'
1694 ENDIF
1695 INFO(1)=-13
1696 INFO(2)=N
1697 GOTO 35
1698 ENDIF
1699 NB_BYTES = NB_BYTES + int(N,8)*K16_8
1700 NB_BYTES_MAX = max(NB_BYTES_MAX,NB_BYTES)
1701 ELSE
1702 SCALING => scaling_data%SCALING
1703 ENDIF
1704 35 CONTINUE
1705 CALL MUMPS_PROPINFO( ICNTL(1), INFO(1),
1706 & COMM, MYID )
1707.LT. IF (INFO(1) 0) GOTO 90
1708 CALL MPI_BCAST( SCALING(1), N, MPI_DOUBLE_PRECISION,
1709 & MASTER, COMM, IERR_MPI)
1710 IF ( I_AM_SLAVE ) THEN
1711 DO I = 1, LILOC
1712.GE..AND..LE. IF (ILOC(I) 1 ILOC(I) N) THEN
1713 scaling_data%SCALING_LOC(I) = SCALING(ILOC(I))
1714 ENDIF
1715 ENDDO
1716 ENDIF
1717 90 CONTINUE
1718.NE. IF (MYID MASTER) THEN
1719 IF (associated(SCALING)) THEN
1720 DEALLOCATE(SCALING)
1721 NB_BYTES = NB_BYTES - int(N,8)*K16_8
1722 ENDIF
1723 ENDIF
1724 NULLIFY(SCALING)
1725.LT. IF (INFO(1) 0) THEN
1726 IF (associated(scaling_data%SCALING_LOC)) THEN
1727 DEALLOCATE(scaling_data%SCALING_LOC)
1728 NULLIFY(scaling_data%SCALING_LOC)
1729 ENDIF
1730 ENDIF
1731 RETURN
1732 END SUBROUTINE DMUMPS_SET_SCALING_LOC
mumps_abort
#define mumps_abort
Definition VE_Metis.h:25
dmumps_solve_fwd_trsolve
subroutine dmumps_solve_fwd_trsolve(a, la, apos, npiv, ldadiag, nrhs_b, wcb, lwcb, lda_wcb, ppiv_courant, mtype, keep)
Definition dsol_aux.F:1147
dmumps_solve_bwd_panels
subroutine dmumps_solve_bwd_panels(a, la, apos, npiv, iw, nrhs_b, wcb, lwcb, lda_wcb, ppiv_courant, mtype, keep)
Definition dsol_aux.F:1274
dmumps_eltyd
subroutine dmumps_eltyd(mtype, n, nelt, eltptr, leltvar, eltvar, na_elt8, a_elt, saverhs, x, y, w, k50)
Definition dsol_aux.F:653
dmumps_eltqd2
subroutine dmumps_eltqd2(mtype, n, nelt, eltptr, leltvar, eltvar, na_elt8, a_elt, lhs, wrhs, w, rhs, keep, keep8)
Definition dsol_aux.F:509
dmumps_sol_x_elt
subroutine dmumps_sol_x_elt(mtype, n, nelt, eltptr, leltvar, eltvar, na_elt8, a_elt, w, keep, keep8)
Definition dsol_aux.F:529
dmumps_freetopso
subroutine dmumps_freetopso(n, keep28, iwcb, liww, w, lwc, poswcb, iwposcb, ptricb, ptracb)
Definition dsol_aux.F:17
dmumps_solve_bwd_trsolve
subroutine dmumps_solve_bwd_trsolve(a, la, apos, npiv, ldadiag, nrhs_b, wcb, lwcb, lda_wcb, ppiv_courant, mtype, keep)
Definition dsol_aux.F:1185
dmumps_scal_x
subroutine dmumps_scal_x(a, nz8, n, irn, icn, z, keep, keep8, colsca, eff_size_schur, sym_perm)
Definition dsol_aux.F:174
dmumps_sol_omega
subroutine dmumps_sol_omega(n, rhs, x, y, r_w, c_w, iw, iflag, omega, noiter, testconv, lp, arret, grain)
Definition dsol_aux.F:861
dmumps_sol_scalx_elt
subroutine dmumps_sol_scalx_elt(mtype, n, nelt, eltptr, leltvar, eltvar, na_elt8, a_elt, w, keep, keep8, colsca)
Definition dsol_aux.F:588
dmumps_sol_x
subroutine dmumps_sol_x(a, nz8, n, irn, icn, z, keep, keep8, eff_size_schur, sym_perm)
Definition dsol_aux.F:88
dmumps_solve_gemm_update
subroutine dmumps_solve_gemm_update(a, la, apos1, nx, lda, ny, nrhs_b, wcb, lwcb, ptrx, ldx, ptry, ldy, mtype, keep, coef_y)
Definition dsol_aux.F:1326
dmumps_sol_cpy_fs2rhscomp
subroutine dmumps_sol_cpy_fs2rhscomp(jbdeb, jbfin, nbrows, keep, rhscomp, nrhs, lrhscomp, first_row_rhscomp, w, ld_w, first_row_w)
Definition dsol_aux.F:1040
dmumps_sol_mulr
subroutine dmumps_sol_mulr(n, r, w)
Definition dsol_aux.F:294
dmumps_solve_get_ooc_node
subroutine dmumps_solve_get_ooc_node(inode, ptrfac, keep, a, la, step, keep8, n, must_be_permuted, ierr)
Definition dsol_aux.F:732
dmumps_sol_b
subroutine dmumps_sol_b(n, kase, x, est, w, iw, grain)
Definition dsol_aux.F:304
dmumps_sol_bwd_gthr
subroutine dmumps_sol_bwd_gthr(jbdeb, jbfin, j1, j2, rhscomp, nrhs, lrhscomp, w, ld_w, first_row_w, iw, liw, keep, n, posinrhscomp_bwd)
Definition dsol_aux.F:1063
dmumps_qd2
subroutine dmumps_qd2(mtype, n, nz8, aspk, irn, icn, lhs, wrhs, w, rhs, keep, keep8)
Definition dsol_aux.F:422
dmumps_compso
subroutine dmumps_compso(n, keep28, iwcb, liww, w, lwc, poswcb, iwposcb, ptricb, ptracb)
Definition dsol_aux.F:37
dmumps_sol_y
subroutine dmumps_sol_y(a, nz8, n, irn, icn, rhs, x, r, w, keep, keep8)
Definition dsol_aux.F:224
dmumps_sol_q
subroutine dmumps_sol_q(mtype, iflag, n, lhs, wrhs, w, res, givnorm, anorm, xnorm, sclnrm, mprint, icntl, keep, keep8)
Definition dsol_aux.F:1089
dmumps_build_mapping_info
subroutine dmumps_build_mapping_info(id)
Definition dsol_aux.F:775
dmumps_sol_lcond
subroutine dmumps_sol_lcond(n, rhs, x, y, d, r_w, c_w, iw, kase, omega, erx, cond, lp, keep, keep8)
Definition dsol_aux.F:939
dmumps_solve_fwd_panels
subroutine dmumps_solve_fwd_panels(a, la, apos, npiv, iw, nrhs_b, wcb, lwcb, lda_wcb, ppiv_courant, mtype, keep)
Definition dsol_aux.F:1225
dmumps_mv_elt
subroutine dmumps_mv_elt(n, nelt, eltptr, eltvar, a_elt, x, y, k50, mtype)
Definition dsol_matvec.F:16
dtrsv
subroutine dtrsv(uplo, trans, diag, n, a, lda, x, incx)
DTRSV
Definition dtrsv.f:143
dgemv
subroutine dgemv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy)
DGEMV
Definition dgemv.f:156
dgemm
subroutine dgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
DGEMM
Definition dgemm.f:187
dtrsm
subroutine dtrsm(side, uplo, transa, diag, m, n, alpha, a, lda, b, ldb)
DTRSM
Definition dtrsm.f:181
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
dmumps_ooc
Definition dmumps_ooc.F:14
dmumps_ooc::ooc_node_permuted
integer ooc_node_permuted
Definition dmumps_ooc.F:24
dmumps_ooc::ooc_node_not_in_mem
integer ooc_node_not_in_mem
Definition dmumps_ooc.F:24
dmumps_ooc::dmumps_read_ooc
subroutine, public dmumps_read_ooc(dest, inode, ierr)
Definition dmumps_ooc.F:394
dmumps_ooc::dmumps_solve_modify_state_node
subroutine dmumps_solve_modify_state_node(inode)
Definition dmumps_ooc.F:1388
dmumps_ooc::dmumps_solve_is_inode_in_mem
integer function dmumps_solve_is_inode_in_mem(inode, ptrfac, nsteps, a, la, ierr)
Definition dmumps_ooc.F:1326
dmumps_ooc::dmumps_solve_alloc_factor_space
subroutine, public dmumps_solve_alloc_factor_space(inode, ptrfac, keep, keep8, a, ierr)
Definition dmumps_ooc.F:1475
dmumps_struc_def
Definition dmumps_struc_def.F:14
arret
subroutine arret(nn)
Definition arret.F:87
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