1C 2C This file is part of MUMPS 5.1.2, released 3C on Mon Oct 2 07:37:01 UTC 2017 4C 5C 6C Copyright 1991-2017 CERFACS, CNRS, ENS Lyon, INP Toulouse, Inria, 7C 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 http://www.cecill.info/licences/Licence_CeCILL-C_V1-en.html 12C 13 SUBROUTINE SMUMPS_FREETOPSO( N, KEEP28, IWCB, LIWW, 14 & W, LWC, 15 & POSWCB,IWPOSCB,PTRICB,PTRACB) 16 IMPLICIT NONE 17 INTEGER(8), INTENT(IN) :: LWC 18 INTEGER(8), INTENT(INOUT) :: POSWCB 19 INTEGER N,LIWW,IWPOSCB, KEEP28 20 INTEGER IWCB(LIWW),PTRICB(KEEP28) 21 INTEGER(8) :: PTRACB(KEEP28) 22 REAL W(LWC) 23 INTEGER SIZFI, SIZFR 24 IF ( IWPOSCB .eq. LIWW ) RETURN 25 DO WHILE ( IWCB( IWPOSCB + 2 ) .eq. 0 ) 26 SIZFR = IWCB( IWPOSCB + 1 ) 27 SIZFI = 2 28 IWPOSCB = IWPOSCB + SIZFI 29 POSWCB = POSWCB + SIZFR 30 IF ( IWPOSCB .eq. LIWW ) RETURN 31 END DO 32 RETURN 33 END SUBROUTINE SMUMPS_FREETOPSO 34 SUBROUTINE SMUMPS_COMPSO(N,KEEP28,IWCB,LIWW,W,LWC, 35 & POSWCB,IWPOSCB,PTRICB,PTRACB) 36 IMPLICIT NONE 37 INTEGER(8), INTENT(IN) :: LWC 38 INTEGER(8), INTENT(INOUT) :: POSWCB 39 INTEGER N,LIWW,IWPOSCB,KEEP28 40 INTEGER IWCB(LIWW),PTRICB(KEEP28) 41 INTEGER(8) :: PTRACB(KEEP28) 42 REAL W(LWC) 43 INTEGER IPTIW,SIZFI,LONGI 44 INTEGER(8) :: IPTA, LONGR, SIZFR, I8 45 INTEGER :: I 46 IPTIW = IWPOSCB 47 IPTA = POSWCB 48 LONGI = 0 49 LONGR = 0_8 50 IF ( IPTIW .EQ. LIWW ) RETURN 5110 CONTINUE 52 IF (IWCB(IPTIW+2).EQ.0) THEN 53 SIZFR = int(IWCB(IPTIW+1),8) 54 SIZFI = 2 55 IF (LONGI.NE.0) THEN 56 DO 20 I=0,LONGI-1 57 IWCB(IPTIW + SIZFI - I) = IWCB (IPTIW - I) 58 20 CONTINUE 59 DO 30 I8=0,LONGR-1 60 W(IPTA + SIZFR - I8) = W(IPTA - I8) 61 30 CONTINUE 62 ENDIF 63 DO 40 I=1,KEEP28 64 IF ((PTRICB(I).LE.(IPTIW+1)).AND. 65 & (PTRICB(I).GT.IWPOSCB) ) THEN 66 PTRICB(I) = PTRICB(I) + SIZFI 67 PTRACB(I) = PTRACB(I) + SIZFR 68 ENDIF 6940 CONTINUE 70 IWPOSCB = IWPOSCB + SIZFI 71 IPTIW = IPTIW + SIZFI 72 POSWCB = POSWCB + SIZFR 73 IPTA = IPTA + SIZFR 74 ELSE 75 SIZFR = int(IWCB(IPTIW+1),8) 76 SIZFI = 2 77 IPTIW = IPTIW + SIZFI 78 LONGI = LONGI + SIZFI 79 IPTA = IPTA + SIZFR 80 LONGR = LONGR + SIZFR 81 ENDIF 82 IF (IPTIW.NE.LIWW) GOTO 10 83 RETURN 84 END SUBROUTINE SMUMPS_COMPSO 85 SUBROUTINE SMUMPS_SOL_X(A, NZ8, N, IRN, ICN, Z, KEEP,KEEP8) 86 INTEGER N, I, J, KEEP(500) 87 INTEGER(8), INTENT(IN) :: NZ8 88 INTEGER(8) KEEP8(150) 89 INTEGER IRN(NZ8), ICN(NZ8) 90 REAL A(NZ8) 91 REAL Z(N) 92 REAL, PARAMETER :: ZERO = 0.0E0 93 INTEGER(8) :: K 94 INTRINSIC abs 95 DO 10 I = 1, N 96 Z(I) = ZERO 97 10 CONTINUE 98 IF (KEEP(264).EQ.0) THEN 99 IF (KEEP(50) .EQ.0) THEN 100 DO K = 1_8, NZ8 101 I = IRN(K) 102 J = ICN(K) 103 IF ((I .LT. 1) .OR. (I .GT. N)) CYCLE 104 IF ((J .LT. 1) .OR. (J .GT. N)) CYCLE 105 Z(I) = Z(I) + abs(A(K)) 106 ENDDO 107 ELSE 108 DO K = 1_8, NZ8 109 I = IRN(K) 110 J = ICN(K) 111 IF ((I .LT. 1) .OR. (I .GT. N)) CYCLE 112 IF ((J .LT. 1) .OR. (J .GT. N)) CYCLE 113 Z(I) = Z(I) + abs(A(K)) 114 IF (J.NE.I) THEN 115 Z(J) = Z(J) + abs(A(K)) 116 ENDIF 117 ENDDO 118 ENDIF 119 ELSE 120 IF (KEEP(50) .EQ.0) THEN 121 DO K = 1_8, NZ8 122 I = IRN(K) 123 J = ICN(K) 124 Z(I) = Z(I) + abs(A(K)) 125 ENDDO 126 ELSE 127 DO K = 1_8, NZ8 128 I = IRN(K) 129 J = ICN(K) 130 Z(I) = Z(I) + abs(A(K)) 131 IF (J.NE.I) THEN 132 Z(J) = Z(J) + abs(A(K)) 133 ENDIF 134 ENDDO 135 ENDIF 136 ENDIF 137 RETURN 138 END SUBROUTINE SMUMPS_SOL_X 139 SUBROUTINE SMUMPS_SCAL_X(A, NZ8, N, IRN, ICN, Z, 140 & KEEP, KEEP8, COLSCA) 141 INTEGER, INTENT(IN) :: N, KEEP(500) 142 INTEGER(8), INTENT(IN) :: NZ8 143 INTEGER(8), INTENT(IN) :: KEEP8(150) 144 INTEGER, INTENT(IN) :: IRN(NZ8), ICN(NZ8) 145 REAL, INTENT(IN) :: A(NZ8) 146 REAL, INTENT(IN) :: COLSCA(N) 147 REAL, INTENT(OUT) :: Z(N) 148 REAL, PARAMETER :: ZERO = 0.0E0 149 INTEGER :: I, J 150 INTEGER(8) :: K 151 DO 10 I = 1, N 152 Z(I) = ZERO 153 10 CONTINUE 154 IF (KEEP(50) .EQ.0) THEN 155 DO K = 1_8, NZ8 156 I = IRN(K) 157 J = ICN(K) 158 IF ((I .LT. 1) .OR. (I .GT. N)) CYCLE 159 IF ((J .LT. 1) .OR. (J .GT. N)) CYCLE 160 Z(I) = Z(I) + abs(A(K)*COLSCA(J)) 161 ENDDO 162 ELSE 163 DO K = 1, NZ8 164 I = IRN(K) 165 J = ICN(K) 166 IF ((I .LT. 1) .OR. (I .GT. N)) CYCLE 167 IF ((J .LT. 1) .OR. (J .GT. N)) CYCLE 168 Z(I) = Z(I) + abs(A(K)*COLSCA(J)) 169 IF (J.NE.I) THEN 170 Z(J) = Z(J) + abs(A(K)*COLSCA(I)) 171 ENDIF 172 ENDDO 173 ENDIF 174 RETURN 175 END SUBROUTINE SMUMPS_SCAL_X 176 SUBROUTINE SMUMPS_SOL_Y(A, NZ8, N, IRN, ICN, RHS, X, R, W, 177 & KEEP,KEEP8) 178 IMPLICIT NONE 179 INTEGER, INTENT(IN) :: N, KEEP(500) 180 INTEGER(8), INTENT(IN) :: NZ8 181 INTEGER(8), INTENT(IN) :: KEEP8(150) 182 INTEGER, INTENT(IN) :: IRN(NZ8), ICN(NZ8) 183 REAL, INTENT(IN) :: A(NZ8), RHS(N), X(N) 184 REAL, INTENT(OUT) :: W(N) 185 REAL, INTENT(OUT) :: R(N) 186 INTEGER I, J 187 INTEGER(8) :: K8 188 REAL, PARAMETER :: ZERO = 0.0E0 189 REAL D 190 DO I = 1, N 191 R(I) = RHS(I) 192 W(I) = ZERO 193 ENDDO 194 IF (KEEP(264).EQ.0) THEN 195 IF (KEEP(50) .EQ.0) THEN 196 DO K8 = 1_8, NZ8 197 I = IRN(K8) 198 J = ICN(K8) 199 IF ((I .GT. N) .OR. (J .GT. N) .OR. (I .LT. 1) .OR. 200 & (J .LT. 1)) CYCLE 201 D = A(K8) * X(J) 202 R(I) = R(I) - D 203 W(I) = W(I) + abs(D) 204 ENDDO 205 ELSE 206 DO K8 = 1_8, NZ8 207 I = IRN(K8) 208 J = ICN(K8) 209 IF ((I .GT. N) .OR. (J .GT. N) .OR. (I .LT. 1) .OR. 210 & (J .LT. 1)) CYCLE 211 D = A(K8) * X(J) 212 R(I) = R(I) - D 213 W(I) = W(I) + abs(D) 214 IF (I.NE.J) THEN 215 D = A(K8) * X(I) 216 R(J) = R(J) - D 217 W(J) = W(J) + abs(D) 218 ENDIF 219 ENDDO 220 ENDIF 221 ELSE 222 IF (KEEP(50) .EQ.0) THEN 223 DO K8 = 1_8, NZ8 224 I = IRN(K8) 225 J = ICN(K8) 226 D = A(K8) * X(J) 227 R(I) = R(I) - D 228 W(I) = W(I) + abs(D) 229 ENDDO 230 ELSE 231 DO K8 = 1_8, NZ8 232 I = IRN(K8) 233 J = ICN(K8) 234 D = A(K8) * X(J) 235 R(I) = R(I) - D 236 W(I) = W(I) + abs(D) 237 IF (I.NE.J) THEN 238 D = A(K8) * X(I) 239 R(J) = R(J) - D 240 W(J) = W(J) + abs(D) 241 ENDIF 242 ENDDO 243 ENDIF 244 ENDIF 245 RETURN 246 END SUBROUTINE SMUMPS_SOL_Y 247 SUBROUTINE SMUMPS_SOL_MULR(N, R, W) 248 INTEGER, intent(in) :: N 249 REAL, intent(in) :: W(N) 250 REAL, intent(inout) :: R(N) 251 INTEGER I 252 DO 10 I = 1, N 253 R(I) = R(I) * W(I) 254 10 CONTINUE 255 RETURN 256 END SUBROUTINE SMUMPS_SOL_MULR 257 SUBROUTINE SMUMPS_SOL_B(N, KASE, X, EST, W, IW) 258 INTEGER, intent(in) :: N 259 INTEGER, intent(inout) :: KASE 260 INTEGER IW(N) 261 REAL W(N), X(N) 262 REAL, intent(inout) :: EST 263 INTRINSIC abs, nint, real, sign 264 INTEGER SMUMPS_IXAMAX 265 EXTERNAL SMUMPS_IXAMAX 266 INTEGER ITMAX 267 PARAMETER (ITMAX = 5) 268 INTEGER I, ITER, J, JLAST, JUMP 269 REAL ALTSGN 270 REAL TEMP 271 SAVE ITER, J, JLAST, JUMP 272 REAL ZERO, ONE 273 PARAMETER( ZERO = 0.0E0 ) 274 PARAMETER( ONE = 1.0E0 ) 275 REAL, PARAMETER :: RZERO = 0.0E0 276 REAL, PARAMETER :: RONE = 1.0E0 277 IF (KASE .EQ. 0) THEN 278 DO 10 I = 1, N 279 X(I) = ONE / real(N) 280 10 CONTINUE 281 KASE = 1 282 JUMP = 1 283 RETURN 284 ENDIF 285 SELECT CASE (JUMP) 286 CASE (1) 287 GOTO 20 288 CASE(2) 289 GOTO 40 290 CASE(3) 291 GOTO 70 292 CASE(4) 293 GOTO 120 294 CASE(5) 295 GOTO 160 296 CASE DEFAULT 297 END SELECT 298 20 CONTINUE 299 IF (N .EQ. 1) THEN 300 W(1) = X(1) 301 EST = abs(W(1)) 302 GOTO 190 303 ENDIF 304 DO 30 I = 1, N 305 X(I) = sign( RONE,real(X(I)) ) 306 IW(I) = nint(real(X(I))) 307 30 CONTINUE 308 KASE = 2 309 JUMP = 2 310 RETURN 311 40 CONTINUE 312 J = SMUMPS_IXAMAX(N, X, 1) 313 ITER = 2 314 50 CONTINUE 315 DO 60 I = 1, N 316 X(I) = ZERO 317 60 CONTINUE 318 X(J) = ONE 319 KASE = 1 320 JUMP = 3 321 RETURN 322 70 CONTINUE 323 DO 80 I = 1, N 324 W(I) = X(I) 325 80 CONTINUE 326 DO 90 I = 1, N 327 IF (nint(sign(RONE, real(X(I)))) .NE. IW(I)) GOTO 100 328 90 CONTINUE 329 GOTO 130 330 100 CONTINUE 331 DO 110 I = 1, N 332 X(I) = sign(RONE, real(X(I))) 333 IW(I) = nint(real(X(I))) 334 110 CONTINUE 335 KASE = 2 336 JUMP = 4 337 RETURN 338 120 CONTINUE 339 JLAST = J 340 J = SMUMPS_IXAMAX(N, X, 1) 341 IF ((abs(X(JLAST)) .NE. abs(X(J))) .AND. (ITER .LT. ITMAX)) THEN 342 ITER = ITER + 1 343 GOTO 50 344 ENDIF 345 130 CONTINUE 346 EST = RZERO 347 DO 140 I = 1, N 348 EST = EST + abs(W(I)) 349 140 CONTINUE 350 ALTSGN = RONE 351 DO 150 I = 1, N 352 X(I) = ALTSGN * (RONE + real(I - 1) / real(N - 1)) 353 ALTSGN = -ALTSGN 354 150 CONTINUE 355 KASE = 1 356 JUMP = 5 357 RETURN 358 160 CONTINUE 359 TEMP = RZERO 360 DO 170 I = 1, N 361 TEMP = TEMP + abs(X(I)) 362 170 CONTINUE 363 TEMP = 2.0E0 * TEMP / real(3 * N) 364 IF (TEMP .GT. EST) THEN 365 DO 180 I = 1, N 366 W(I) = X(I) 367 180 CONTINUE 368 EST = TEMP 369 ENDIF 370 190 KASE = 0 371 RETURN 372 END SUBROUTINE SMUMPS_SOL_B 373 SUBROUTINE SMUMPS_QD2( MTYPE, N, NZ8, ASPK, IRN, ICN, 374 & LHS, WRHS, W, RHS, KEEP,KEEP8) 375 IMPLICIT NONE 376 INTEGER MTYPE, N 377 INTEGER(8), INTENT(IN) :: NZ8 378 INTEGER, INTENT(IN) :: IRN( NZ8 ), ICN( NZ8 ) 379 INTEGER KEEP(500) 380 INTEGER(8) KEEP8(150) 381 REAL, INTENT(IN) :: ASPK( NZ8 ) 382 REAL, INTENT(IN) :: LHS( N ), WRHS( N ) 383 REAL, INTENT(OUT):: RHS( N ) 384 REAL, INTENT(OUT):: W( N ) 385 INTEGER I, J 386 INTEGER(8) :: K8 387 REAL, PARAMETER :: DZERO = 0.0E0 388 DO I = 1, N 389 W(I) = DZERO 390 RHS(I) = WRHS(I) 391 ENDDO 392 IF ( KEEP(50) .EQ. 0 ) THEN 393 IF (MTYPE .EQ. 1) THEN 394 IF (KEEP(264).EQ.0) THEN 395 DO K8 = 1_8, NZ8 396 I = IRN(K8) 397 J = ICN(K8) 398 IF ((I .LE. 0) .OR. (I .GT. N) .OR. (J .LE. 0) .OR. 399 & (J .GT. N)) CYCLE 400 RHS(I) = RHS(I) - ASPK(K8) * LHS(J) 401 W(I) = W(I) + abs(ASPK(K8)) 402 ENDDO 403 ELSE 404 DO K8 = 1_8, NZ8 405 I = IRN(K8) 406 J = ICN(K8) 407 RHS(I) = RHS(I) - ASPK(K8) * LHS(J) 408 W(I) = W(I) + abs(ASPK(K8)) 409 ENDDO 410 ENDIF 411 ELSE 412 IF (KEEP(264).EQ.0) THEN 413 DO K8 = 1_8, NZ8 414 I = IRN(K8) 415 J = ICN(K8) 416 IF ((I .LE. 0) .OR. (I .GT. N) .OR. (J .LE. 0) .OR. 417 & (J .GT. N)) CYCLE 418 RHS(J) = RHS(J) - ASPK(K8) * LHS(I) 419 W(J) = W(J) + abs(ASPK(K8)) 420 ENDDO 421 ELSE 422 DO K8 = 1_8, NZ8 423 I = IRN(K8) 424 J = ICN(K8) 425 RHS(J) = RHS(J) - ASPK(K8) * LHS(I) 426 W(J) = W(J) + abs(ASPK(K8)) 427 ENDDO 428 ENDIF 429 ENDIF 430 ELSE 431 IF (KEEP(264).EQ.0) THEN 432 DO K8 = 1_8, NZ8 433 I = IRN(K8) 434 J = ICN(K8) 435 IF ((I .LE. 0) .OR. (I .GT. N) .OR. (J .LE. 0) .OR. 436 & (J .GT. N)) CYCLE 437 RHS(I) = RHS(I) - ASPK(K8) * LHS(J) 438 W(I) = W(I) + abs(ASPK(K8)) 439 IF (J.NE.I) THEN 440 RHS(J) = RHS(J) - ASPK(K8) * LHS(I) 441 W(J) = W(J) + abs(ASPK(K8)) 442 ENDIF 443 ENDDO 444 ELSE 445 DO K8 = 1_8, NZ8 446 I = IRN(K8) 447 J = ICN(K8) 448 RHS(I) = RHS(I) - ASPK(K8) * LHS(J) 449 W(I) = W(I) + abs(ASPK(K8)) 450 IF (J.NE.I) THEN 451 RHS(J) = RHS(J) - ASPK(K8) * LHS(I) 452 W(J) = W(J) + abs(ASPK(K8)) 453 ENDIF 454 ENDDO 455 ENDIF 456 ENDIF 457 RETURN 458 END SUBROUTINE SMUMPS_QD2 459 SUBROUTINE SMUMPS_ELTQD2( MTYPE, N, 460 & NELT, ELTPTR, LELTVAR, ELTVAR, NA_ELT8, A_ELT, 461 & LHS, WRHS, W, RHS, KEEP,KEEP8 ) 462 IMPLICIT NONE 463 INTEGER MTYPE, N, NELT, LELTVAR 464 INTEGER(8), INTENT(IN) :: NA_ELT8 465 INTEGER ELTPTR(NELT+1), ELTVAR(LELTVAR) 466 INTEGER KEEP(500) 467 INTEGER(8) KEEP8(150) 468 REAL A_ELT(NA_ELT8) 469 REAL LHS( N ), WRHS( N ), RHS( N ) 470 REAL W(N) 471 CALL SMUMPS_MV_ELT(N, NELT, ELTPTR, ELTVAR, A_ELT, 472 & LHS, RHS, KEEP(50), MTYPE ) 473 RHS = WRHS - RHS 474 CALL SMUMPS_SOL_X_ELT( MTYPE, N, 475 & NELT, ELTPTR, LELTVAR, ELTVAR, NA_ELT8, A_ELT, 476 & W, KEEP,KEEP8 ) 477 RETURN 478 END SUBROUTINE SMUMPS_ELTQD2 479 SUBROUTINE SMUMPS_SOL_X_ELT( MTYPE, N, 480 & NELT, ELTPTR, LELTVAR, ELTVAR, NA_ELT8, A_ELT, 481 & W, KEEP,KEEP8 ) 482 IMPLICIT NONE 483 INTEGER MTYPE, N, NELT, LELTVAR 484 INTEGER(8), INTENT(IN) :: NA_ELT8 485 INTEGER ELTPTR(NELT+1), ELTVAR(LELTVAR) 486 INTEGER KEEP(500) 487 INTEGER(8) KEEP8(150) 488 REAL A_ELT(NA_ELT8) 489 REAL TEMP 490 REAL W(N) 491 INTEGER I, J, IEL, SIZEI, IELPTR 492 INTEGER(8) :: K8 493 REAL DZERO 494 PARAMETER(DZERO = 0.0E0) 495 W = DZERO 496 K8 = 1_8 497 DO IEL = 1, NELT 498 SIZEI = ELTPTR( IEL + 1 ) - ELTPTR( IEL ) 499 IELPTR = ELTPTR( IEL ) - 1 500 IF ( KEEP(50).EQ.0 ) THEN 501 IF (MTYPE.EQ.1) THEN 502 DO J = 1, SIZEI 503 DO I = 1, SIZEI 504 W( ELTVAR( IELPTR + I) ) = 505 & W( ELTVAR( IELPTR + I) ) 506 & + abs(A_ELT( K8 )) 507 K8 = K8 + 1_8 508 END DO 509 END DO 510 ELSE 511 DO J = 1, SIZEI 512 TEMP = W( ELTVAR( IELPTR + J ) ) 513 DO I = 1, SIZEI 514 TEMP = TEMP + abs( A_ELT(K8)) 515 K8 = K8 + 1_8 516 END DO 517 W(ELTVAR( IELPTR + J )) = 518 & W(ELTVAR( IELPTR + J )) + TEMP 519 END DO 520 ENDIF 521 ELSE 522 DO J = 1, SIZEI 523 W(ELTVAR( IELPTR + J )) = 524 & W(ELTVAR( IELPTR + J )) + abs(A_ELT( K8 )) 525 K8 = K8 + 1_8 526 DO I = J+1, SIZEI 527 W(ELTVAR( IELPTR + J )) = 528 & W(ELTVAR( IELPTR + J )) + abs(A_ELT( K8 )) 529 W(ELTVAR( IELPTR + I ) ) = 530 & W(ELTVAR( IELPTR + I )) + abs(A_ELT( K8 )) 531 K8 = K8 + 1_8 532 END DO 533 ENDDO 534 ENDIF 535 ENDDO 536 RETURN 537 END SUBROUTINE SMUMPS_SOL_X_ELT 538 SUBROUTINE SMUMPS_SOL_SCALX_ELT(MTYPE, N, 539 & NELT, ELTPTR, LELTVAR, ELTVAR, NA_ELT8, A_ELT, 540 & W, KEEP,KEEP8, COLSCA ) 541 IMPLICIT NONE 542 INTEGER MTYPE, N, NELT, LELTVAR 543 INTEGER(8), INTENT(IN) :: NA_ELT8 544 INTEGER ELTPTR(NELT+1), ELTVAR(LELTVAR) 545 INTEGER KEEP(500) 546 INTEGER(8) KEEP8(150) 547 REAL COLSCA(N) 548 REAL A_ELT(NA_ELT8) 549 REAL W(N) 550 REAL TEMP, TEMP2 551 INTEGER I, J, IEL, SIZEI, IELPTR 552 INTEGER(8) :: K8 553 REAL DZERO 554 PARAMETER(DZERO = 0.0E0) 555 W = DZERO 556 K8 = 1_8 557 DO IEL = 1, NELT 558 SIZEI = ELTPTR( IEL + 1 ) - ELTPTR( IEL ) 559 IELPTR = ELTPTR( IEL ) - 1 560 IF ( KEEP(50).EQ.0 ) THEN 561 IF (MTYPE.EQ.1) THEN 562 DO J = 1, SIZEI 563 TEMP2 = abs(COLSCA(ELTVAR( IELPTR + J) )) 564 DO I = 1, SIZEI 565 W( ELTVAR( IELPTR + I) ) = 566 & W( ELTVAR( IELPTR + I) ) 567 & + abs(A_ELT( K8 )) * TEMP2 568 K8 = K8 + 1_8 569 END DO 570 END DO 571 ELSE 572 DO J = 1, SIZEI 573 TEMP = W( ELTVAR( IELPTR + J ) ) 574 TEMP2= abs(COLSCA(ELTVAR( IELPTR + J) )) 575 DO I = 1, SIZEI 576 TEMP = TEMP + abs(A_ELT( K8 )) * TEMP2 577 K8 = K8 + 1_8 578 END DO 579 W(ELTVAR( IELPTR + J )) = 580 & W(ELTVAR( IELPTR + J )) + TEMP 581 END DO 582 ENDIF 583 ELSE 584 DO J = 1, SIZEI 585 W(ELTVAR( IELPTR + J )) = 586 & W(ELTVAR( IELPTR + J )) + 587 & abs( A_ELT( K8 )*COLSCA(ELTVAR( IELPTR + J)) ) 588 K8 = K8 + 1_8 589 DO I = J+1, SIZEI 590 W(ELTVAR( IELPTR + J )) = 591 & W(ELTVAR( IELPTR + J )) + 592 & abs(A_ELT( K8 )*COLSCA(ELTVAR( IELPTR + J))) 593 W(ELTVAR( IELPTR + I ) ) = 594 & W(ELTVAR( IELPTR + I )) + 595 & abs(A_ELT( K8 )*COLSCA(ELTVAR( IELPTR + I))) 596 K8 = K8 + 1_8 597 END DO 598 ENDDO 599 ENDIF 600 ENDDO 601 RETURN 602 END SUBROUTINE SMUMPS_SOL_SCALX_ELT 603 SUBROUTINE SMUMPS_ELTYD( MTYPE, N, NELT, ELTPTR, 604 & LELTVAR, ELTVAR, NA_ELT8, A_ELT, 605 & SAVERHS, X, Y, W, K50 ) 606 IMPLICIT NONE 607 INTEGER N, NELT, K50, MTYPE, LELTVAR 608 INTEGER(8) :: NA_ELT8 609 INTEGER ELTPTR( NELT + 1 ), ELTVAR( LELTVAR ) 610 REAL A_ELT( NA_ELT8 ), X( N ), Y( N ), 611 & SAVERHS(N) 612 REAL W(N) 613 INTEGER IEL, I , J, K, SIZEI, IELPTR 614 REAL ZERO 615 REAL TEMP 616 REAL TEMP2 617 PARAMETER( ZERO = 0.0E0 ) 618 Y = SAVERHS 619 W = ZERO 620 K = 1 621 DO IEL = 1, NELT 622 SIZEI = ELTPTR( IEL + 1 ) - ELTPTR( IEL ) 623 IELPTR = ELTPTR( IEL ) - 1 624 IF ( K50 .eq. 0 ) THEN 625 IF ( MTYPE .eq. 1 ) THEN 626 DO J = 1, SIZEI 627 TEMP = X( ELTVAR( IELPTR + J ) ) 628 DO I = 1, SIZEI 629 Y( ELTVAR( IELPTR + I ) ) = 630 & Y( ELTVAR( IELPTR + I ) ) - 631 & A_ELT( K ) * TEMP 632 W( ELTVAR( IELPTR + I ) ) = 633 & W( ELTVAR( IELPTR + I ) ) + 634 & abs( A_ELT( K ) * TEMP ) 635 K = K + 1 636 END DO 637 END DO 638 ELSE 639 DO J = 1, SIZEI 640 TEMP = Y( ELTVAR( IELPTR + J ) ) 641 TEMP2 = W( ELTVAR( IELPTR + J ) ) 642 DO I = 1, SIZEI 643 TEMP = TEMP - 644 & A_ELT( K ) * X( ELTVAR( IELPTR + I ) ) 645 TEMP2 = TEMP2 + abs( 646 & A_ELT( K ) * X( ELTVAR( IELPTR + I ) ) ) 647 K = K + 1 648 END DO 649 Y( ELTVAR( IELPTR + J ) ) = TEMP 650 W( ELTVAR( IELPTR + J ) ) = TEMP2 651 END DO 652 END IF 653 ELSE 654 DO J = 1, SIZEI 655 Y( ELTVAR( IELPTR + J ) ) = 656 & Y( ELTVAR( IELPTR + J ) ) - 657 & A_ELT( K ) * X( ELTVAR( IELPTR + J ) ) 658 W( ELTVAR( IELPTR + J ) ) = 659 & W( ELTVAR( IELPTR + J ) ) + abs( 660 & A_ELT( K ) * X( ELTVAR( IELPTR + J ) ) ) 661 K = K + 1 662 DO I = J+1, SIZEI 663 Y( ELTVAR( IELPTR + I ) ) = 664 & Y( ELTVAR( IELPTR + I ) ) - 665 & A_ELT( K ) * X( ELTVAR( IELPTR + J ) ) 666 Y( ELTVAR( IELPTR + J ) ) = 667 & Y( ELTVAR( IELPTR + J ) ) - 668 & A_ELT( K ) * X( ELTVAR( IELPTR + I ) ) 669 W( ELTVAR( IELPTR + I ) ) = 670 & W( ELTVAR( IELPTR + I ) ) + abs( 671 & A_ELT( K ) * X( ELTVAR( IELPTR + J ) ) ) 672 W( ELTVAR( IELPTR + J ) ) = 673 & W( ELTVAR( IELPTR + J ) ) + abs( 674 & A_ELT( K ) * X( ELTVAR( IELPTR + I ) ) ) 675 K = K + 1 676 END DO 677 END DO 678 END IF 679 END DO 680 RETURN 681 END SUBROUTINE SMUMPS_ELTYD 682 SUBROUTINE SMUMPS_SOLVE_GET_OOC_NODE( 683 & INODE,PTRFAC,KEEP,A,LA,STEP, 684 & KEEP8,N,MUST_BE_PERMUTED,IERR) 685 USE SMUMPS_OOC 686 IMPLICIT NONE 687 INTEGER INODE,KEEP(500),N 688 INTEGER(8) KEEP8(150) 689 INTEGER(8) :: LA 690 INTEGER(8) :: PTRFAC(KEEP(28)) 691 INTEGER STEP(N) 692 INTEGER IERR 693 REAL A(LA) 694 INTEGER RETURN_VALUE 695 LOGICAL MUST_BE_PERMUTED 696 RETURN_VALUE=SMUMPS_SOLVE_IS_INODE_IN_MEM(INODE,PTRFAC, 697 & KEEP(28),A,LA,IERR) 698 IF(RETURN_VALUE.EQ.OOC_NODE_NOT_IN_MEM)THEN 699 IF(IERR.LT.0)THEN 700 RETURN 701 ENDIF 702 CALL SMUMPS_SOLVE_ALLOC_FACTOR_SPACE(INODE,PTRFAC, 703 & KEEP,KEEP8,A,IERR) 704 IF(IERR.LT.0)THEN 705 RETURN 706 ENDIF 707 CALL SMUMPS_READ_OOC( 708 & A(PTRFAC(STEP(INODE))), 709 & INODE,IERR 710 & ) 711 IF(IERR.LT.0)THEN 712 RETURN 713 ENDIF 714 ELSE 715 IF(IERR.LT.0)THEN 716 RETURN 717 ENDIF 718 ENDIF 719 IF(RETURN_VALUE.NE.OOC_NODE_PERMUTED)THEN 720 MUST_BE_PERMUTED=.TRUE. 721 CALL SMUMPS_SOLVE_MODIFY_STATE_NODE(INODE) 722 ELSE 723 MUST_BE_PERMUTED=.FALSE. 724 ENDIF 725 RETURN 726 END SUBROUTINE SMUMPS_SOLVE_GET_OOC_NODE 727 SUBROUTINE SMUMPS_BUILD_MAPPING_INFO(id) 728 USE SMUMPS_STRUC_DEF 729 IMPLICIT NONE 730 INCLUDE 'mpif.h' 731 TYPE(SMUMPS_STRUC), TARGET :: id 732 INTEGER, ALLOCATABLE, DIMENSION(:) :: LOCAL_LIST 733 INTEGER :: I,IERR,TMP,NSTEPS,N_LOCAL_LIST 734 INTEGER :: MASTER,TAG_SIZE,TAG_LIST 735 INTEGER :: STATUS(MPI_STATUS_SIZE) 736 LOGICAL :: I_AM_SLAVE 737 PARAMETER(MASTER=0, TAG_SIZE=85,TAG_LIST=86) 738 I_AM_SLAVE = (id%MYID .NE. MASTER 739 & .OR. ((id%MYID.EQ.MASTER).AND.(id%KEEP(46).EQ.1))) 740 NSTEPS = id%KEEP(28) 741 ALLOCATE(LOCAL_LIST(NSTEPS),STAT=IERR) 742 IF(IERR.GT.0) THEN 743 WRITE(*,*)'Problem in solve: error allocating LOCAL_LIST' 744 CALL MUMPS_ABORT() 745 END IF 746 N_LOCAL_LIST = 0 747 IF(I_AM_SLAVE) THEN 748 DO I=1,NSTEPS 749 IF(id%PTLUST_S(I).NE.0) THEN 750 N_LOCAL_LIST = N_LOCAL_LIST + 1 751 LOCAL_LIST(N_LOCAL_LIST) = I 752 END IF 753 END DO 754 IF(id%MYID.NE.MASTER) THEN 755 CALL MPI_SEND(N_LOCAL_LIST, 1, 756 & MPI_INTEGER, MASTER, TAG_SIZE, id%COMM,IERR) 757 CALL MPI_SEND(LOCAL_LIST, N_LOCAL_LIST, 758 & MPI_INTEGER, MASTER, TAG_LIST, id%COMM,IERR) 759 DEALLOCATE(LOCAL_LIST) 760 ALLOCATE(id%IPTR_WORKING(1), 761 & id%WORKING(1), 762 & STAT=IERR) 763 IF(IERR.GT.0) THEN 764 WRITE(*,*)'Problem in solve: error allocating ', 765 & 'IPTR_WORKING and WORKING' 766 CALL MUMPS_ABORT() 767 END IF 768 END IF 769 END IF 770 IF(id%MYID.EQ.MASTER) THEN 771 ALLOCATE(id%IPTR_WORKING(id%NPROCS+1), STAT=IERR) 772 IF(IERR.GT.0) THEN 773 WRITE(*,*)'Problem in solve: error allocating IPTR_WORKING' 774 CALL MUMPS_ABORT() 775 END IF 776 id%IPTR_WORKING = 0 777 id%IPTR_WORKING(1) = 1 778 id%IPTR_WORKING(MASTER+2) = N_LOCAL_LIST 779 DO I=1, id%NPROCS-1 780 CALL MPI_RECV(TMP, 1, MPI_INTEGER, MPI_ANY_SOURCE, 781 & TAG_SIZE, id%COMM, STATUS, IERR) 782 id%IPTR_WORKING(STATUS(MPI_SOURCE)+2) = TMP 783 END DO 784 DO I=2, id%NPROCS+1 785 id%IPTR_WORKING(I) = id%IPTR_WORKING(I) 786 & + id%IPTR_WORKING(I-1) 787 END DO 788 ALLOCATE(id%WORKING(id%IPTR_WORKING(id%NPROCS+1)-1),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 TMP = MASTER + 1 794 IF (I_AM_SLAVE) THEN 795 id%WORKING(id%IPTR_WORKING(TMP):id%IPTR_WORKING(TMP+1)-1) 796 & = LOCAL_LIST(1:id%IPTR_WORKING(TMP+1) 797 & -id%IPTR_WORKING(TMP)) 798 ENDIF 799 DO I=1,id%NPROCS-1 800 CALL MPI_RECV(LOCAL_LIST, NSTEPS, MPI_INTEGER, 801 & MPI_ANY_SOURCE, TAG_LIST, id%COMM, STATUS, IERR) 802 TMP = STATUS(MPI_SOURCE)+1 803 id%WORKING(id%IPTR_WORKING(TMP):id%IPTR_WORKING(TMP+1)-1) 804 & = LOCAL_LIST(1:id%IPTR_WORKING(TMP+1)- 805 & id%IPTR_WORKING(TMP)) 806 END DO 807 DEALLOCATE(LOCAL_LIST) 808 END IF 809 END SUBROUTINE SMUMPS_BUILD_MAPPING_INFO 810 SUBROUTINE SMUMPS_SOL_OMEGA(N, RHS, 811 & X, Y, R_W, C_W, IW, IFLAG, 812 & OMEGA, NOITER, TESTConv, 813 & LP, ARRET ) 814 IMPLICIT NONE 815 INTEGER N, IFLAG 816 INTEGER IW(N,2) 817 REAL RHS(N) 818 REAL X(N), Y(N) 819 REAL R_W(N,2) 820 REAL C_W(N) 821 INTEGER LP, NOITER 822 LOGICAL TESTConv 823 REAL OMEGA(2) 824 REAL ARRET 825 REAL, PARAMETER :: CGCE=0.2E0 826 REAL, PARAMETER :: CTAU=1.0E3 827 INTEGER I, IMAX 828 REAL OM1, OM2, DXMAX 829 REAL TAU, DD 830 REAL OLDOMG(2) 831 REAL, PARAMETER :: ZERO=0.0E0 832 REAL, PARAMETER :: ONE=1.0E0 833 INTEGER SMUMPS_IXAMAX 834 INTRINSIC abs, max 835 SAVE OM1, OLDOMG 836 IMAX = SMUMPS_IXAMAX(N, X, 1) 837 DXMAX = abs(X(IMAX)) 838 OMEGA(1) = ZERO 839 OMEGA(2) = ZERO 840 DO I = 1, N 841 TAU = (R_W(I, 2) * DXMAX + abs(RHS(I))) * real(N) * CTAU 842 DD = R_W(I, 1) + abs(RHS(I)) 843 IF (DD .GT. TAU * epsilon(CTAU)) THEN 844 OMEGA(1) = max(OMEGA(1), abs(Y(I)) / DD) 845 IW(I, 1) = 1 846 ELSE 847 IF (TAU .GT. ZERO) THEN 848 OMEGA(2) = max(OMEGA(2), 849 & abs(Y(I)) / (DD + R_W(I, 2) * DXMAX)) 850 ENDIF 851 IW(I, 1) = 2 852 ENDIF 853 ENDDO 854 IF (TESTConv) THEN 855 OM2 = OMEGA(1) + OMEGA(2) 856 IF (OM2 .LT. ARRET ) THEN 857 IFLAG = 1 858 GOTO 70 859 ENDIF 860 IF (NOITER .GE. 1) THEN 861 IF (OM2 .GT. OM1 * CGCE) THEN 862 IF (OM2 .GT. OM1) THEN 863 OMEGA(1) = OLDOMG(1) 864 OMEGA(2) = OLDOMG(2) 865 DO I = 1, N 866 X(I) = C_W(I) 867 ENDDO 868 IFLAG = 2 869 GOTO 70 870 ENDIF 871 IFLAG = 3 872 GOTO 70 873 ENDIF 874 ENDIF 875 DO I = 1, N 876 C_W(I) = X(I) 877 ENDDO 878 OLDOMG(1) = OMEGA(1) 879 OLDOMG(2) = OMEGA(2) 880 OM1 = OM2 881 ENDIF 882 IFLAG = 0 883 RETURN 884 70 CONTINUE 885 RETURN 886 END SUBROUTINE SMUMPS_SOL_OMEGA 887 SUBROUTINE SMUMPS_SOL_LCOND(N, RHS, 888 & X, Y, D, R_W, C_W, IW, KASE, 889 & OMEGA, ERX, COND, 890 & LP, KEEP,KEEP8 ) 891 IMPLICIT NONE 892 INTEGER N, KASE, KEEP(500) 893 INTEGER(8) KEEP8(150) 894 INTEGER IW(N,2) 895 REAL RHS(N) 896 REAL X(N), Y(N) 897 REAL D(N) 898 REAL R_W(N,2) 899 REAL C_W(N) 900 INTEGER LP 901 REAL COND(2),OMEGA(2) 902 LOGICAL LCOND1, LCOND2 903 INTEGER JUMP, I, IMAX 904 REAL ERX, DXMAX 905 REAL DXIMAX 906 REAL, PARAMETER :: ZERO = 0.0E0 907 REAL, PARAMETER :: ONE = 1.0E0 908 INTEGER SMUMPS_IXAMAX 909 INTRINSIC abs, max 910 SAVE LCOND1, LCOND2, JUMP, DXIMAX, DXMAX 911 IF (KASE .EQ. 0) THEN 912 LCOND1 = .FALSE. 913 LCOND2 = .FALSE. 914 COND(1) = ONE 915 COND(2) = ONE 916 ERX = ZERO 917 JUMP = 1 918 ENDIF 919 SELECT CASE (JUMP) 920 CASE (1) 921 GOTO 30 922 CASE(2) 923 GOTO 10 924 CASE(3) 925 GOTO 110 926 CASE(4) 927 GOTO 150 928 CASE(5) 929 GOTO 35 930 CASE DEFAULT 931 END SELECT 932 10 CONTINUE 933 30 CONTINUE 934 35 CONTINUE 935 IMAX = SMUMPS_IXAMAX(N, X, 1) 936 DXMAX = abs(X(IMAX)) 937 DO I = 1, N 938 IF (IW(I, 1) .EQ. 1) THEN 939 R_W(I, 1) = R_W(I, 1) + abs(RHS(I)) 940 R_W(I, 2) = ZERO 941 LCOND1 = .TRUE. 942 ELSE 943 R_W(I, 2) = R_W(I, 2) * DXMAX + R_W(I, 1) 944 R_W(I, 1) = ZERO 945 LCOND2 = .TRUE. 946 ENDIF 947 ENDDO 948 DO I = 1, N 949 C_W(I) = X(I) * D(I) 950 ENDDO 951 IMAX = SMUMPS_IXAMAX(N, C_W(1), 1) 952 DXIMAX = abs(C_W(IMAX)) 953 IF (.NOT.LCOND1) GOTO 130 954 100 CONTINUE 955 CALL SMUMPS_SOL_B(N, KASE, Y, COND(1), C_W, IW(1, 2)) 956 IF (KASE .EQ. 0) GOTO 120 957 IF (KASE .EQ. 1) CALL SMUMPS_SOL_MULR(N, Y, D) 958 IF (KASE .EQ. 2) CALL SMUMPS_SOL_MULR(N, Y, R_W) 959 JUMP = 3 960 RETURN 961 110 CONTINUE 962 IF (KASE .EQ. 1) CALL SMUMPS_SOL_MULR(N, Y, R_W) 963 IF (KASE .EQ. 2) CALL SMUMPS_SOL_MULR(N, Y, D) 964 GOTO 100 965 120 CONTINUE 966 IF (DXIMAX .GT. ZERO) COND(1) = COND(1) / DXIMAX 967 ERX = OMEGA(1) * COND(1) 968 130 CONTINUE 969 IF (.NOT.LCOND2) GOTO 170 970 KASE = 0 971 140 CONTINUE 972 CALL SMUMPS_SOL_B(N, KASE, Y, COND(2), C_W, IW(1, 2)) 973 IF (KASE .EQ. 0) GOTO 160 974 IF (KASE .EQ. 1) CALL SMUMPS_SOL_MULR(N, Y, D) 975 IF (KASE .EQ. 2) CALL SMUMPS_SOL_MULR(N, Y, R_W(1, 2)) 976 JUMP = 4 977 RETURN 978 150 CONTINUE 979 IF (KASE .EQ. 1) CALL SMUMPS_SOL_MULR(N, Y, R_W(1, 2)) 980 IF (KASE .EQ. 2) CALL SMUMPS_SOL_MULR(N, Y, D) 981 GOTO 140 982 160 IF (DXIMAX .GT. ZERO) THEN 983 COND(2) = COND(2) / DXIMAX 984 ENDIF 985 ERX = ERX + OMEGA(2) * COND(2) 986 170 CONTINUE 987 RETURN 988 END SUBROUTINE SMUMPS_SOL_LCOND 989 SUBROUTINE SMUMPS_SOL_CPY_FS2RHSCOMP( JBDEB, JBFIN, NBROWS, 990 & KEEP, RHSCOMP, NRHS, LRHSCOMP, FIRST_ROW_RHSCOMP, W, LD_W, 991 & FIRST_ROW_W ) 992 INTEGER :: JBDEB, JBFIN, NBROWS 993 INTEGER :: NRHS, LRHSCOMP 994 INTEGER :: FIRST_ROW_RHSCOMP 995 INTEGER, INTENT(IN) :: KEEP(500) 996# if defined(RHSCOMP_BYROWS) 997 REAL, INTENT(INOUT) :: RHSCOMP(NRHS,LRHSCOMP) 998# else 999 REAL, INTENT(INOUT) :: RHSCOMP(LRHSCOMP,NRHS) 1000# endif 1001 INTEGER :: LD_W, FIRST_ROW_W 1002 REAL :: W(LD_W*(JBFIN-JBDEB+1)) 1003 INTEGER :: JJ, K, ISHIFT 1004#if defined(RHSCOMP_BYROWS) 1005!$OMP PARALLEL DO PRIVATE (ISHIFT, K), IF 1006!$OMP& ((NBROWS) * (JBFIN-JBDEB+1) > KEEP(363)) 1007 DO JJ = 0, NBROWS-1 1008 ISHIFT = FIRST_ROW_W+JJ 1009 DO K = JBDEB, JBFIN 1010 RHSCOMP(K,FIRST_ROW_RHSCOMP+JJ) = 1011 & W(ISHIFT+LD_W*(K-JBDEB)) 1012 END DO 1013 END DO 1014!$OMP END PARALLEL DO 1015#else 1016!$OMP PARALLEL DO PRIVATE(ISHIFT, JJ), IF 1017!$OMP& (JBFIN-JBDEB+1 > 2*KEEP(362) .AND. 1018!$OMP& NBROWS * (JBFIN-JBDEB+1) > 2*KEEP(363)) 1019 DO K = JBDEB, JBFIN 1020 ISHIFT = FIRST_ROW_W + LD_W * (K-JBDEB) 1021 DO JJ = 0, NBROWS-1 1022 RHSCOMP(FIRST_ROW_RHSCOMP+JJ,K) = W(ISHIFT+JJ) 1023 END DO 1024 END DO 1025!$OMP END PARALLEL DO 1026#endif 1027 RETURN 1028 END SUBROUTINE SMUMPS_SOL_CPY_FS2RHSCOMP 1029 SUBROUTINE SMUMPS_SOL_BWD_GTHR( JBDEB, JBFIN, J1, J2, 1030 & RHSCOMP, NRHS, LRHSCOMP, W, LD_W, FIRST_ROW_W, 1031 & IW, LIW, KEEP, N, POSINRHSCOMP_BWD ) 1032 INTEGER, INTENT(IN) :: JBDEB, JBFIN, J1, J2 1033 INTEGER, INTENT(IN) :: NRHS, LRHSCOMP 1034 INTEGER, INTENT(IN) :: FIRST_ROW_W, LD_W, LIW 1035 INTEGER, INTENT(IN) :: IW(LIW) 1036 INTEGER, INTENT(IN) :: KEEP(500) 1037# if defined(RHSCOMP_BYROWS) 1038 REAL, INTENT(INOUT) :: RHSCOMP(NRHS,LRHSCOMP) 1039# else 1040 REAL, INTENT(INOUT) :: RHSCOMP(LRHSCOMP,NRHS) 1041# endif 1042 REAL :: W(LD_W*(JBFIN-JBDEB+1)) 1043 INTEGER, INTENT(IN) :: N 1044 INTEGER, INTENT(IN) :: POSINRHSCOMP_BWD(N) 1045 INTEGER :: ISHIFT, JJ, K, IPOSINRHSCOMP 1046#if defined(RHSCOMP_BYROWS) 1047!$OMP PARALLEL DO PRIVATE(K,ISHIFT,IPOSINRHSCOMP), IF 1048!$OMP& ((JBFIN-JBDEB+1)*(J2-KEEP(253)-J1+1)>KEEP(363)) 1049 DO JJ = J1, J2-KEEP(253) 1050 ISHIFT = FIRST_ROW_W+JJ-J1 1051 IPOSINRHSCOMP = abs(POSINRHSCOMP_BWD(IW(JJ))) 1052 DO K=JBDEB, JBFIN 1053 W(ISHIFT+(K-JBDEB)*LD_W) = RHSCOMP(K,IPOSINRHSCOMP) 1054 ENDDO 1055 ENDDO 1056!$OMP END PARALLEL DO 1057#else 1058!$OMP PARALLEL DO PRIVATE(JJ,ISHIFT,IPOSINRHSCOMP), IF 1059!$OMP& ((JBFIN-JBDEB+1 > 2*KEEP(362) .AND. 1060!$OMP& (JBFIN-JBDEB+1)*(J2-KEEP(253)-J1+1)>2*KEEP(363))) 1061 DO K=JBDEB, JBFIN 1062 ISHIFT = FIRST_ROW_W+(K-JBDEB)*LD_W 1063 DO JJ = J1, J2-KEEP(253) 1064 IPOSINRHSCOMP = abs(POSINRHSCOMP_BWD(IW(JJ))) 1065 W(ISHIFT+JJ-J1)= RHSCOMP(IPOSINRHSCOMP,K) 1066 ENDDO 1067 ENDDO 1068!$OMP END PARALLEL DO 1069#endif 1070 RETURN 1071 END SUBROUTINE SMUMPS_SOL_BWD_GTHR 1072 SUBROUTINE SMUMPS_SOL_Q(MTYPE, IFLAG, N, 1073 & LHS, WRHS, W, RES, GIVNORM, ANORM, XNORM, SCLNRM, 1074 & MPRINT, ICNTL, KEEP,KEEP8) 1075 INTEGER MTYPE,N,IFLAG,ICNTL(40), KEEP(500) 1076 INTEGER(8) KEEP8(150) 1077 REAL RES(N),LHS(N) 1078 REAL WRHS(N) 1079 REAL W(N) 1080 REAL RESMAX,RESL2,XNORM, SCLNRM 1081 REAL ANORM,DZERO 1082 LOGICAL GIVNORM,PROK 1083 INTEGER MPRINT, MP 1084 INTEGER K 1085 INTRINSIC abs, max, sqrt 1086 MP = ICNTL(2) 1087 PROK = (MPRINT .GT. 0) 1088 DZERO = 0.0E0 1089 IF (.NOT.GIVNORM) ANORM = DZERO 1090 RESMAX = DZERO 1091 RESL2 = DZERO 1092 DO 40 K = 1, N 1093 RESMAX = max(RESMAX, abs(RES(K))) 1094 RESL2 = RESL2 + abs(RES(K)) * abs(RES(K)) 1095 IF (.NOT.GIVNORM) ANORM = max(ANORM, W(K)) 1096 40 CONTINUE 1097 XNORM = DZERO 1098 DO 50 K = 1, N 1099 XNORM = max(XNORM, abs(LHS(K))) 1100 50 CONTINUE 1101 IF ( XNORM .EQ. DZERO .OR. (exponent(XNORM) .LT. 1102 & minexponent(XNORM) + KEEP(122) ) 1103 & .OR. 1104 & ( exponent(ANORM)+exponent(XNORM) .LT. 1105 & minexponent(XNORM) + KEEP(122) ) 1106 & .OR. 1107 & ( exponent(ANORM) + exponent(XNORM) -exponent(RESMAX) 1108 & .LT. minexponent(XNORM) + KEEP(122) ) 1109 & ) THEN 1110 IF (mod(IFLAG/2,2) .EQ. 0) THEN 1111 IFLAG = IFLAG + 2 1112 ENDIF 1113 IF ((MP .GT. 0) .AND. (ICNTL(4) .GE. 2)) WRITE( MP, * ) 1114 & ' max-NORM of computed solut. is zero or close to zero. ' 1115 ENDIF 1116 IF (RESMAX .EQ. DZERO) THEN 1117 SCLNRM = DZERO 1118 ELSE 1119 SCLNRM = RESMAX / (ANORM * XNORM) 1120 ENDIF 1121 RESL2 = sqrt(RESL2) 1122 IF (PROK) WRITE( MPRINT, 90 ) RESMAX, RESL2, ANORM, XNORM, 1123 & SCLNRM 1124 90 FORMAT (/' RESIDUAL IS ............ (MAX-NORM) =',1PD9.2/ 1125 & ' .. (2-NORM) =',1PD9.2/ 1126 & ' RINFOG(4):NORM OF input Matrix (MAX-NORM)=',1PD9.2/ 1127 & ' RINFOG(5):NORM OF Computed SOLUT (MAX-NORM)=',1PD9.2/ 1128 & ' RINFOG(6):SCALED RESIDUAL ...... (MAX-NORM)=',1PD9.2) 1129 RETURN 1130 END SUBROUTINE SMUMPS_SOL_Q 1131