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 13C 14C History: 15C ------- 16C This maximum transversal set of routines are 17C based on the work done by Jacko Koster at CERFACS for 18C his PhD thesis from Institut National Polytechnique de Toulouse 19C at CERFACS (1995-1997) and includes modifications provided 20C by the author as well as work done by Stephane Pralet 21C first at CERFACS during his PhD thesis (2003-2004) then 22C at INPT-IRIT (2004-2005) during his post-doctoral position. 23C 24C The main research publication references for this work are: 25C [1] I. S. Duff, (1981), 26C "Algorithm 575. Permutations for a zero-free diagonal", 27C ACM Trans. Math. Software 7(3), 387-390. 28C [2] I. S. Duff and J. Koster, (1998), 29C "The design and use of algorithms for permuting large 30C entries to the diagonal of sparse matrices", 31C SIAM J. Matrix Anal. Appl., vol. 20, no. 4, pp. 889-901. 32C [3] I. S. Duff and J. Koster, (2001), 33C "On algorithms for permuting large entries to the diagonal 34C of sparse matrices", 35C SIAM J. Matrix Anal. Appl., vol. 22, no. 4, pp. 973-996. 36C 37 SUBROUTINE ZMUMPS_MTRANSI(ICNTL,CNTL) 38 IMPLICIT NONE 39 INTEGER NICNTL, NCNTL 40 PARAMETER (NICNTL=10, NCNTL=10) 41 INTEGER ICNTL(NICNTL) 42 DOUBLE PRECISION CNTL(NCNTL) 43 INTEGER I 44 ICNTL(1) = 6 45 ICNTL(2) = 6 46 ICNTL(3) = -1 47 ICNTL(4) = -1 48 ICNTL(5) = 0 49 DO 10 I = 6,NICNTL 50 ICNTL(I) = 0 51 10 CONTINUE 52 CNTL(1) = 0.0D0 53 CNTL(2) = 0.0D0 54 DO 20 I = 3,NCNTL 55 CNTL(I) = 0.0D0 56 20 CONTINUE 57 RETURN 58 END SUBROUTINE ZMUMPS_MTRANSI 59 SUBROUTINE ZMUMPS_MTRANSB 60 & (M,N,NE,IP,IRN,A,IPERM,NUM,JPERM,PR,Q,L,D,RINF) 61 IMPLICIT NONE 62 INTEGER :: M,N,NUM 63 INTEGER(8), INTENT(IN) :: NE 64 INTEGER :: IRN(NE),IPERM(M),JPERM(N),Q(M),L(M) 65 INTEGER(8), INTENT(IN) :: IP(N+1) 66 INTEGER(8), INTENT(OUT) :: PR(N) 67 DOUBLE PRECISION :: A(NE) 68 DOUBLE PRECISION :: D(M), RINF 69 INTEGER :: I,II,J,JJ,JORD,Q0,QLEN,IDUM,JDUM,ISP,JSP, 70 & I0,UP,LOW, IK 71 INTEGER(8) :: K,KK,KK1,KK2 72 DOUBLE PRECISION CSP,DI,DNEW,DQ0,AI,A0,BV,TBV,RLX 73 DOUBLE PRECISION ZERO,MINONE,ONE 74 PARAMETER (ZERO=0.0D0,MINONE=-1.0D0,ONE=1.0D0) 75 INTRINSIC abs,min 76 EXTERNAL ZMUMPS_MTRANSD, ZMUMPS_MTRANSE, 77 & ZMUMPS_MTRANSF, ZMUMPS_MTRANSX 78 RLX = D(1) 79 NUM = 0 80 BV = RINF 81 DO 10 I = 1,N 82 JPERM(I) = 0 83 PR(I) = IP(I) 84 10 CONTINUE 85 DO 12 I = 1,M 86 IPERM(I) = 0 87 D(I) = ZERO 88 12 CONTINUE 89 DO 30 J = 1,N 90 A0 = MINONE 91 DO 20 K = IP(J),IP(J+1)-1_8 92 I = IRN(K) 93 AI = abs(A(K)) 94 IF (AI.GT.D(I)) D(I) = AI 95 IF (JPERM(J).NE.0) GO TO 20 96 IF (AI.GE.BV) THEN 97 A0 = BV 98 IF (IPERM(I).NE.0) GO TO 20 99 JPERM(J) = I 100 IPERM(I) = J 101 NUM = NUM + 1 102 ELSE 103 IF (AI.LE.A0) GO TO 20 104 A0 = AI 105 I0 = I 106 ENDIF 107 20 CONTINUE 108 IF (A0.NE.MINONE .AND. A0.LT.BV) THEN 109 BV = A0 110 IF (IPERM(I0).NE.0) GO TO 30 111 IPERM(I0) = J 112 JPERM(J) = I0 113 NUM = NUM + 1 114 ENDIF 115 30 CONTINUE 116 IF (M.EQ.N) THEN 117 DO 35 I = 1,M 118 BV = min(BV,D(I)) 119 35 CONTINUE 120 ENDIF 121 IF (NUM.EQ.N) GO TO 1000 122 DO 95 J = 1,N 123 IF (JPERM(J).NE.0) GO TO 95 124 DO 50 K = IP(J),IP(J+1)-1_8 125 I = IRN(K) 126 AI = abs(A(K)) 127 IF (AI.LT.BV) GO TO 50 128 IF (IPERM(I).EQ.0) GO TO 90 129 JJ = IPERM(I) 130 KK1 = PR(JJ) 131 KK2 = IP(JJ+1) - 1_8 132 IF (KK1.GT.KK2) GO TO 50 133 DO 70 KK = KK1,KK2 134 II = IRN(KK) 135 IF (IPERM(II).NE.0) GO TO 70 136 IF (abs(A(KK)).GE.BV) GO TO 80 137 70 CONTINUE 138 PR(JJ) = KK2 + 1_8 139 50 CONTINUE 140 GO TO 95 141 80 JPERM(JJ) = II 142 IPERM(II) = JJ 143 PR(JJ) = KK + 1_8 144 90 NUM = NUM + 1 145 JPERM(J) = I 146 IPERM(I) = J 147 PR(J) = K + 1_8 148 95 CONTINUE 149 IF (NUM.EQ.N) GO TO 1000 150 DO 99 I = 1,M 151 D(I) = MINONE 152 L(I) = 0 153 99 CONTINUE 154 TBV = BV * (ONE-RLX) 155 DO 100 JORD = 1,N 156 IF (JPERM(JORD).NE.0) GO TO 100 157 QLEN = 0 158 LOW = M + 1 159 UP = M + 1 160 CSP = MINONE 161 J = JORD 162 PR(J) = -1_8 163 DO 115 K = IP(J),IP(J+1)-1_8 164 I = IRN(K) 165 DNEW = abs(A(K)) 166 IF (CSP.GE.DNEW) GO TO 115 167 IF (IPERM(I).EQ.0) THEN 168 CSP = DNEW 169 ISP = I 170 JSP = J 171 IF (CSP.GE.TBV) GO TO 160 172 ELSE 173 D(I) = DNEW 174 IF (DNEW.GE.TBV) THEN 175 LOW = LOW - 1 176 Q(LOW) = I 177 ELSE 178 QLEN = QLEN + 1 179 L(I) = QLEN 180 CALL ZMUMPS_MTRANSD(I,M,Q,D,L,1) 181 ENDIF 182 JJ = IPERM(I) 183 PR(JJ) = int(J,8) 184 ENDIF 185 115 CONTINUE 186 DO 150 JDUM = 1,NUM 187 IF (LOW.EQ.UP) THEN 188 IF (QLEN.EQ.0) GO TO 160 189 I = Q(1) 190 IF (CSP.GE.D(I)) GO TO 160 191 BV = D(I) 192 TBV = BV * (ONE-RLX) 193 DO 152 IDUM = 1,M 194 CALL ZMUMPS_MTRANSE(QLEN,M,Q,D,L,1) 195 L(I) = 0 196 LOW = LOW - 1 197 Q(LOW) = I 198 IF (QLEN.EQ.0) GO TO 153 199 I = Q(1) 200 IF (D(I).LT.TBV) GO TO 153 201 152 CONTINUE 202 ENDIF 203 153 UP = UP - 1 204 Q0 = Q(UP) 205 DQ0 = D(Q0) 206 L(Q0) = UP 207 J = IPERM(Q0) 208 DO 155 K = IP(J),IP(J+1)-1_8 209 I = IRN(K) 210 IF (L(I).GE.UP) GO TO 155 211 DNEW = min(DQ0,abs(A(K))) 212 IF (CSP.GE.DNEW) GO TO 155 213 IF (IPERM(I).EQ.0) THEN 214 CSP = DNEW 215 ISP = I 216 JSP = J 217 IF (CSP.GE.TBV) GO TO 160 218 ELSE 219 DI = D(I) 220 IF (DI.GE.TBV .OR. DI.GE.DNEW) GO TO 155 221 D(I) = DNEW 222 IF (DNEW.GE.TBV) THEN 223 IF (DI.NE.MINONE) THEN 224 CALL ZMUMPS_MTRANSF(L(I),QLEN,M,Q,D,L,1) 225 ENDIF 226 L(I) = 0 227 LOW = LOW - 1 228 Q(LOW) = I 229 ELSE 230 IF (DI.EQ.MINONE) THEN 231 QLEN = QLEN + 1 232 L(I) = QLEN 233 ENDIF 234 CALL ZMUMPS_MTRANSD(I,M,Q,D,L,1) 235 ENDIF 236 JJ = IPERM(I) 237 PR(JJ) = int(J,8) 238 ENDIF 239 155 CONTINUE 240 150 CONTINUE 241 160 IF (CSP.EQ.MINONE) GO TO 190 242 BV = min(BV,CSP) 243 TBV = BV * (ONE-RLX) 244 NUM = NUM + 1 245 I = ISP 246 J = JSP 247 DO 170 JDUM = 1,NUM+1 248 I0 = JPERM(J) 249 JPERM(J) = I 250 IPERM(I) = J 251 J = int(PR(J)) 252 IF (J.EQ.-1) GO TO 190 253 I = I0 254 170 CONTINUE 255 190 DO 191 IK = UP,M 256 I = Q(IK) 257 D(I) = MINONE 258 L(I) = 0 259 191 CONTINUE 260 DO 192 IK = LOW,UP-1 261 I = Q(IK) 262 D(I) = MINONE 263 192 CONTINUE 264 DO 193 IK = 1,QLEN 265 I = Q(IK) 266 D(I) = MINONE 267 L(I) = 0 268 193 CONTINUE 269 100 CONTINUE 270 1000 IF (M.EQ.N .and. NUM.EQ.N) GO TO 2000 271 CALL ZMUMPS_MTRANSX(M,N,IPERM,L,JPERM) 272 2000 RETURN 273 END SUBROUTINE ZMUMPS_MTRANSB 274 SUBROUTINE ZMUMPS_MTRANSD(I,N,Q,D,L,IWAY) 275 IMPLICIT NONE 276 INTEGER I,N,IWAY 277 INTEGER Q(N),L(N) 278 DOUBLE PRECISION D(N) 279 INTEGER IDUM,K,POS,POSK,QK 280 PARAMETER (K=2) 281 DOUBLE PRECISION DI 282 POS = L(I) 283 IF (POS.LE.1) GO TO 20 284 DI = D(I) 285 IF (IWAY.EQ.1) THEN 286 DO 10 IDUM = 1,N 287 POSK = POS/K 288 QK = Q(POSK) 289 IF (DI.LE.D(QK)) GO TO 20 290 Q(POS) = QK 291 L(QK) = POS 292 POS = POSK 293 IF (POS.LE.1) GO TO 20 294 10 CONTINUE 295 ELSE 296 DO 15 IDUM = 1,N 297 POSK = POS/K 298 QK = Q(POSK) 299 IF (DI.GE.D(QK)) GO TO 20 300 Q(POS) = QK 301 L(QK) = POS 302 POS = POSK 303 IF (POS.LE.1) GO TO 20 304 15 CONTINUE 305 ENDIF 306 20 Q(POS) = I 307 L(I) = POS 308 RETURN 309 END SUBROUTINE ZMUMPS_MTRANSD 310 SUBROUTINE ZMUMPS_MTRANSE(QLEN,N,Q,D,L,IWAY) 311 IMPLICIT NONE 312 INTEGER QLEN,N,IWAY 313 INTEGER Q(N),L(N) 314 DOUBLE PRECISION D(N) 315 INTEGER I,IDUM,K,POS,POSK 316 PARAMETER (K=2) 317 DOUBLE PRECISION DK,DR,DI 318 I = Q(QLEN) 319 DI = D(I) 320 QLEN = QLEN - 1 321 POS = 1 322 IF (IWAY.EQ.1) THEN 323 DO 10 IDUM = 1,N 324 POSK = K*POS 325 IF (POSK.GT.QLEN) GO TO 20 326 DK = D(Q(POSK)) 327 IF (POSK.LT.QLEN) THEN 328 DR = D(Q(POSK+1)) 329 IF (DK.LT.DR) THEN 330 POSK = POSK + 1 331 DK = DR 332 ENDIF 333 ENDIF 334 IF (DI.GE.DK) GO TO 20 335 Q(POS) = Q(POSK) 336 L(Q(POS)) = POS 337 POS = POSK 338 10 CONTINUE 339 ELSE 340 DO 15 IDUM = 1,N 341 POSK = K*POS 342 IF (POSK.GT.QLEN) GO TO 20 343 DK = D(Q(POSK)) 344 IF (POSK.LT.QLEN) THEN 345 DR = D(Q(POSK+1)) 346 IF (DK.GT.DR) THEN 347 POSK = POSK + 1 348 DK = DR 349 ENDIF 350 ENDIF 351 IF (DI.LE.DK) GO TO 20 352 Q(POS) = Q(POSK) 353 L(Q(POS)) = POS 354 POS = POSK 355 15 CONTINUE 356 ENDIF 357 20 Q(POS) = I 358 L(I) = POS 359 RETURN 360 END SUBROUTINE ZMUMPS_MTRANSE 361 SUBROUTINE ZMUMPS_MTRANSF(POS0,QLEN,N,Q,D,L,IWAY) 362 IMPLICIT NONE 363 INTEGER POS0,QLEN,N,IWAY 364 INTEGER Q(N),L(N) 365 DOUBLE PRECISION D(N) 366 INTEGER I,IDUM,K,POS,POSK,QK 367 PARAMETER (K=2) 368 DOUBLE PRECISION DK,DR,DI 369 IF (QLEN.EQ.POS0) THEN 370 QLEN = QLEN - 1 371 RETURN 372 ENDIF 373 I = Q(QLEN) 374 DI = D(I) 375 QLEN = QLEN - 1 376 POS = POS0 377 IF (IWAY.EQ.1) THEN 378 IF (POS.LE.1) GO TO 20 379 DO 10 IDUM = 1,N 380 POSK = POS/K 381 QK = Q(POSK) 382 IF (DI.LE.D(QK)) GO TO 20 383 Q(POS) = QK 384 L(QK) = POS 385 POS = POSK 386 IF (POS.LE.1) GO TO 20 387 10 CONTINUE 388 20 Q(POS) = I 389 L(I) = POS 390 IF (POS.NE.POS0) RETURN 391 DO 30 IDUM = 1,N 392 POSK = K*POS 393 IF (POSK.GT.QLEN) GO TO 40 394 DK = D(Q(POSK)) 395 IF (POSK.LT.QLEN) THEN 396 DR = D(Q(POSK+1)) 397 IF (DK.LT.DR) THEN 398 POSK = POSK + 1 399 DK = DR 400 ENDIF 401 ENDIF 402 IF (DI.GE.DK) GO TO 40 403 QK = Q(POSK) 404 Q(POS) = QK 405 L(QK) = POS 406 POS = POSK 407 30 CONTINUE 408 ELSE 409 IF (POS.LE.1) GO TO 34 410 DO 32 IDUM = 1,N 411 POSK = POS/K 412 QK = Q(POSK) 413 IF (DI.GE.D(QK)) GO TO 34 414 Q(POS) = QK 415 L(QK) = POS 416 POS = POSK 417 IF (POS.LE.1) GO TO 34 418 32 CONTINUE 419 34 Q(POS) = I 420 L(I) = POS 421 IF (POS.NE.POS0) RETURN 422 DO 36 IDUM = 1,N 423 POSK = K*POS 424 IF (POSK.GT.QLEN) GO TO 40 425 DK = D(Q(POSK)) 426 IF (POSK.LT.QLEN) THEN 427 DR = D(Q(POSK+1)) 428 IF (DK.GT.DR) THEN 429 POSK = POSK + 1 430 DK = DR 431 ENDIF 432 ENDIF 433 IF (DI.LE.DK) GO TO 40 434 QK = Q(POSK) 435 Q(POS) = QK 436 L(QK) = POS 437 POS = POSK 438 36 CONTINUE 439 ENDIF 440 40 Q(POS) = I 441 L(I) = POS 442 RETURN 443 END SUBROUTINE ZMUMPS_MTRANSF 444 SUBROUTINE ZMUMPS_MTRANSQ(IP,LENL,LENH,W,WLEN,A,NVAL,VAL) 445 IMPLICIT NONE 446 INTEGER ::WLEN,NVAL 447 INTEGER :: LENL(*),LENH(*),W(*) 448 INTEGER(8) :: IP(*) 449 DOUBLE PRECISION :: A(*),VAL 450 INTEGER XX,J,K,S,POS 451 INTEGER(8) :: II 452 PARAMETER (XX=10) 453 DOUBLE PRECISION SPLIT(XX),HA 454 NVAL = 0 455 DO 10 K = 1,WLEN 456 J = W(K) 457 DO 15 II = IP(J)+int(LENL(J),8),IP(J)+int(LENH(J)-1,8) 458 HA = A(II) 459 IF (NVAL.EQ.0) THEN 460 SPLIT(1) = HA 461 NVAL = 1 462 ELSE 463 DO 20 S = NVAL,1,-1 464 IF (SPLIT(S).EQ.HA) GO TO 15 465 IF (SPLIT(S).GT.HA) THEN 466 POS = S + 1 467 GO TO 21 468 ENDIF 469 20 CONTINUE 470 POS = 1 471 21 DO 22 S = NVAL,POS,-1 472 SPLIT(S+1) = SPLIT(S) 473 22 CONTINUE 474 SPLIT(POS) = HA 475 NVAL = NVAL + 1 476 ENDIF 477 IF (NVAL.EQ.XX) GO TO 11 478 15 CONTINUE 479 10 CONTINUE 480 11 IF (NVAL.GT.0) VAL = SPLIT((NVAL+1)/2) 481 RETURN 482 END SUBROUTINE ZMUMPS_MTRANSQ 483 SUBROUTINE ZMUMPS_MTRANSR(N,NE,IP,IRN,A) 484 IMPLICIT NONE 485 INTEGER, INTENT(IN) :: N 486 INTEGER(8), INTENT(IN) :: NE 487 INTEGER(8), INTENT(IN) :: IP(N+1) 488 INTEGER, INTENT(INOUT) :: IRN(NE) 489 DOUBLE PRECISION, INTENT(INOUT) :: A(NE) 490 INTEGER :: THRESH,TDLEN 491 PARAMETER (THRESH=15,TDLEN=50) 492 INTEGER :: J, LEN, HI 493 INTEGER(8) :: K, IPJ, TD, FIRST, LAST, MID, R, S 494 DOUBLE PRECISION :: HA, KEY 495 INTEGER(8) :: TODO(TDLEN) 496 DO 100 J = 1,N 497 LEN = int(IP(J+1) - IP(J)) 498 IF (LEN.LE.1) GO TO 100 499 IPJ = IP(J) 500 IF (LEN.LT.THRESH) GO TO 400 501 TODO(1) = IPJ 502 TODO(2) = IPJ +int(LEN,8) 503 TD = 2_8 504 500 CONTINUE 505 FIRST = TODO(TD-1) 506 LAST = TODO(TD) 507 KEY = A((FIRST+LAST)/2) 508 DO 475 K = FIRST,LAST-1 509 HA = A(K) 510 IF (HA.EQ.KEY) GO TO 475 511 IF (HA.GT.KEY) GO TO 470 512 KEY = HA 513 GO TO 470 514 475 CONTINUE 515 TD = TD - 2_8 516 GO TO 425 517 470 MID = FIRST 518 DO 450 K = FIRST,LAST-1 519 IF (A(K).LE.KEY) GO TO 450 520 HA = A(MID) 521 A(MID) = A(K) 522 A(K) = HA 523 HI = IRN(MID) 524 IRN(MID) = IRN(K) 525 IRN(K) = HI 526 MID = MID + 1 527 450 CONTINUE 528 IF (MID-FIRST.GE.LAST-MID) THEN 529 TODO(TD+2) = LAST 530 TODO(TD+1) = MID 531 TODO(TD) = MID 532 ELSE 533 TODO(TD+2) = MID 534 TODO(TD+1) = FIRST 535 TODO(TD) = LAST 536 TODO(TD-1) = MID 537 ENDIF 538 TD = TD + 2_8 539 425 CONTINUE 540 IF (TD.EQ.0_8) GO TO 400 541 IF (TODO(TD)-TODO(TD-1).GE.int(THRESH,8)) GO TO 500 542 TD = TD - 2_8 543 GO TO 425 544 400 DO 200 R = IPJ+1_8,IPJ+int(LEN-1,8) 545 IF (A(R-1) .LT. A(R)) THEN 546 HA = A(R) 547 HI = IRN(R) 548 A(R) = A(R-1_8) 549 IRN(R) = IRN(R-1_8) 550 DO 300 S = R-1,IPJ+1_8,-1_8 551 IF (A(S-1) .LT. HA) THEN 552 A(S) = A(S-1) 553 IRN(S) = IRN(S-1) 554 ELSE 555 A(S) = HA 556 IRN(S) = HI 557 GO TO 200 558 END IF 559 300 CONTINUE 560 A(IPJ) = HA 561 IRN(IPJ) = HI 562 END IF 563 200 CONTINUE 564 100 CONTINUE 565 RETURN 566 END SUBROUTINE ZMUMPS_MTRANSR 567 SUBROUTINE ZMUMPS_MTRANSS(M,N,NE,IP,IRN,A,IPERM,NUMX, 568 & W,LEN,LENL,LENH,FC,IW,IW4,RLX,RINF) 569 IMPLICIT NONE 570 INTEGER, INTENT(IN) :: M,N 571 INTEGER(8), INTENT(IN) :: NE 572 INTEGER, INTENT(OUT) :: NUMX 573 INTEGER(8), INTENT(IN) :: IP(N+1) 574 INTEGER :: IRN(NE),IPERM(N), 575 & W(N),LEN(N),LENL(N),LENH(N),FC(N),IW(M),IW4(3*N+M) 576 DOUBLE PRECISION A(NE),RLX,RINF 577 INTEGER :: NUM,NVAL,WLEN,I,J,L,CNT,MOD, IDUM 578 INTEGER(8) :: K, II, KDUM1, KDUM2 579 DOUBLE PRECISION :: BVAL,BMIN,BMAX 580 EXTERNAL ZMUMPS_MTRANSQ,ZMUMPS_MTRANSU,ZMUMPS_MTRANSX 581 DO 20 J = 1,N 582 FC(J) = J 583 LEN(J) = int(IP(J+1) - IP(J)) 584 20 CONTINUE 585 DO 21 I = 1,M 586 IW(I) = 0 587 21 CONTINUE 588 CNT = 1 589 MOD = 1 590 NUMX = 0 591 CALL ZMUMPS_MTRANSU(CNT,MOD,M,N,IRN,NE,IP,LEN,FC,IW, 592 & NUMX,N, 593 & IW4(1),IW4(N+1),IW4(2*N+1),IW4(2*N+M+1)) 594 NUM = NUMX 595 IF (NUM.NE.N) THEN 596 BMAX = RINF 597 ELSE 598 BMAX = RINF 599 DO 30 J = 1,N 600 BVAL = 0.0D0 601 DO 25 K = IP(J),IP(J+1)-1_8 602 IF (A(K).GT.BVAL) BVAL = A(K) 603 25 CONTINUE 604 IF (BVAL.LT.BMAX) BMAX = BVAL 605 30 CONTINUE 606 BMAX = 1.001D0 * BMAX 607 ENDIF 608 BVAL = 0.0D0 609 BMIN = 0.0D0 610 WLEN = 0 611 DO 48 J = 1,N 612 L = int(IP(J+1) - IP(J)) 613 LENH(J) = L 614 LEN(J) = L 615 DO 45 K = IP(J),IP(J+1)-1_8 616 IF (A(K).LT.BMAX) GO TO 46 617 45 CONTINUE 618 K = IP(J+1) 619 46 LENL(J) = int(K - IP(J)) 620 IF (LENL(J).EQ.L) GO TO 48 621 WLEN = WLEN + 1 622 W(WLEN) = J 623 48 CONTINUE 624 DO 90 KDUM1 = 1_8,NE 625 IF (NUM.EQ.NUMX) THEN 626 DO 50 I = 1,M 627 IPERM(I) = IW(I) 628 50 CONTINUE 629 DO 80 KDUM2 = 1_8,NE 630 BMIN = BVAL 631 IF (BMAX-BMIN .LE. RLX) GO TO 1000 632 CALL ZMUMPS_MTRANSQ(IP,LENL,LEN,W,WLEN,A,NVAL,BVAL) 633 IF (NVAL.LE.1) GO TO 1000 634 K = 1 635 DO 70 IDUM = 1,N 636 IF (K.GT.WLEN) GO TO 71 637 J = W(K) 638 DO 55 II = IP(J)+int(LEN(J)-1,8), 639 & IP(J)+int(LENL(J),8),-1_8 640 IF (A(II).GE.BVAL) GO TO 60 641 I = IRN(II) 642 IF (IW(I).NE.J) GO TO 55 643 IW(I) = 0 644 NUM = NUM - 1 645 FC(N-NUM) = J 646 55 CONTINUE 647 60 LENH(J) = LEN(J) 648 LEN(J) = int(II - IP(J) + 1) 649 IF (LENL(J).EQ.LENH(J)) THEN 650 W(K) = W(WLEN) 651 WLEN = WLEN - 1 652 ELSE 653 K = K + 1 654 ENDIF 655 70 CONTINUE 656 71 IF (NUM.LT.NUMX) GO TO 81 657 80 CONTINUE 658 81 MOD = 1 659 ELSE 660 BMAX = BVAL 661 IF (BMAX-BMIN .LE. RLX) GO TO 1000 662 CALL ZMUMPS_MTRANSQ(IP,LEN,LENH,W,WLEN,A,NVAL,BVAL) 663 IF (NVAL.EQ.0. OR. BVAL.EQ.BMIN) GO TO 1000 664 K = 1 665 DO 87 IDUM = 1,N 666 IF (K.GT.WLEN) GO TO 88 667 J = W(K) 668 DO 85 II = IP(J)+int(LEN(J),8),IP(J)+int(LENH(J)-1,8) 669 IF (A(II).LT.BVAL) GO TO 86 670 85 CONTINUE 671 86 LENL(J) = LEN(J) 672 LEN(J) = int(II - IP(J)) 673 IF (LENL(J).EQ.LENH(J)) THEN 674 W(K) = W(WLEN) 675 WLEN = WLEN - 1 676 ELSE 677 K = K + 1 678 ENDIF 679 87 CONTINUE 680 88 MOD = 0 681 ENDIF 682 CNT = CNT + 1 683 CALL ZMUMPS_MTRANSU(CNT,MOD,M,N,IRN,NE,IP,LEN,FC,IW, 684 & NUM,NUMX, 685 & IW4(1),IW4(N+1),IW4(2*N+1),IW4(2*N+M+1)) 686 90 CONTINUE 687 1000 IF (M.EQ.N .and. NUMX.EQ.N) GO TO 2000 688 CALL ZMUMPS_MTRANSX(M,N,IPERM,IW,W) 689 2000 RETURN 690 END SUBROUTINE ZMUMPS_MTRANSS 691C 692 SUBROUTINE ZMUMPS_MTRANSU 693 & (ID,MOD,M,N,IRN,LIRN,IP,LENC,FC,IPERM,NUM,NUMX, 694 & PR,ARP,CV,OUT) 695 IMPLICIT NONE 696 INTEGER :: ID,MOD,M,N,NUM,NUMX 697 INTEGER(8), INTENT(IN) :: LIRN 698 INTEGER :: ARP(N),CV(M),IRN(LIRN), 699 & FC(N),IPERM(M),LENC(N),OUT(N),PR(N) 700 INTEGER(8), INTENT(IN) :: IP(N) 701 INTEGER I,J,J1,JORD,NFC,K,KK, 702 & NUM0,NUM1,NUM2,ID0,ID1,LAST 703 INTEGER(8) :: IN1, IN2, II 704 IF (ID.EQ.1) THEN 705 DO 5 I = 1,M 706 CV(I) = 0 707 5 CONTINUE 708 DO 6 J = 1,N 709 ARP(J) = 0 710 6 CONTINUE 711 NUM1 = N 712 NUM2 = N 713 ELSE 714 IF (MOD.EQ.1) THEN 715 DO 8 J = 1,N 716 ARP(J) = 0 717 8 CONTINUE 718 ENDIF 719 NUM1 = NUMX 720 NUM2 = N - NUMX 721 ENDIF 722 NUM0 = NUM 723 NFC = 0 724 ID0 = (ID-1)*N 725 DO 100 JORD = NUM0+1,N 726 ID1 = ID0 + JORD 727 J = FC(JORD-NUM0) 728 PR(J) = -1 729 DO 70 K = 1,JORD 730 IF (ARP(J).GE.LENC(J)) GO TO 30 731 IN1 = IP(J) + int(ARP(J),8) 732 IN2 = IP(J) + int(LENC(J) - 1,8) 733 DO 20 II = IN1,IN2 734 I = IRN(II) 735 IF (IPERM(I).EQ.0) GO TO 80 736 20 CONTINUE 737 ARP(J) = LENC(J) 738 30 OUT(J) = LENC(J) - 1 739 DO 60 KK = 1,JORD 740 IN1 = int(OUT(J),8) 741 IF (IN1.LT.0) GO TO 50 742 IN2 = IP(J) + int(LENC(J) - 1,8) 743 IN1 = IN2 - IN1 744 DO 40 II = IN1,IN2 745 I = IRN(II) 746 IF (CV(I).EQ.ID1) GO TO 40 747 J1 = J 748 J = IPERM(I) 749 CV(I) = ID1 750 PR(J) = J1 751 OUT(J1) = int(IN2 - II) - 1 752 GO TO 70 753 40 CONTINUE 754 50 J1 = PR(J) 755 IF (J1.EQ.-1) THEN 756 NFC = NFC + 1 757 FC(NFC) = J 758 IF (NFC.GT.NUM2) THEN 759 LAST = JORD 760 GO TO 101 761 ENDIF 762 GO TO 100 763 ENDIF 764 J = J1 765 60 CONTINUE 766 70 CONTINUE 767 80 IPERM(I) = J 768 ARP(J) = int(II - IP(J)) + 1 769 NUM = NUM + 1 770 DO 90 K = 1,JORD 771 J = PR(J) 772 IF (J.EQ.-1) GO TO 95 773 II = IP(J) + int(LENC(J) - OUT(J) - 2,8) 774 I = IRN(II) 775 IPERM(I) = J 776 90 CONTINUE 777 95 IF (NUM.EQ.NUM1) THEN 778 LAST = JORD 779 GO TO 101 780 ENDIF 781 100 CONTINUE 782 LAST = N 783 101 DO 110 JORD = LAST+1,N 784 NFC = NFC + 1 785 FC(NFC) = FC(JORD-NUM0) 786 110 CONTINUE 787 RETURN 788 END SUBROUTINE ZMUMPS_MTRANSU 789C 790 SUBROUTINE ZMUMPS_MTRANSW(M,N,NE,IP,IRN,A,IPERM,NUM, 791 & JPERM,L32,OUT,PR,Q,L,U,D,RINF) 792 IMPLICIT NONE 793 INTEGER :: M,N,NUM 794 INTEGER(8), INTENT(IN) :: NE 795 INTEGER :: IRN(NE),IPERM(M),Q(M),L32(max(M,N)) 796 INTEGER(8) :: IP(N+1), PR(N), L(M), JPERM(N), OUT(N) 797 DOUBLE PRECISION A(NE),U(M),D(M),RINF,RINF3 798 INTEGER :: I,I0,II,J,JJ,JORD,Q0,QLEN,JDUM,JSP, 799 & UP,LOW,IK 800 INTEGER(8) :: K, KK, KK1, KK2, K0, K1, K2, ISP 801 DOUBLE PRECISION :: CSP,DI,DMIN,DNEW,DQ0,VJ,RLX 802 LOGICAL :: LORD 803 DOUBLE PRECISION :: ZERO, ONE 804 PARAMETER (ZERO=0.0D0,ONE=1.0D0) 805 EXTERNAL ZMUMPS_MTRANSD, ZMUMPS_MTRANSE, 806 & ZMUMPS_MTRANSF, ZMUMPS_MTRANSX 807 RLX = U(1) 808 RINF3 = U(2) 809 LORD = (JPERM(1).EQ.6) 810 NUM = 0 811 DO 10 I = 1,N 812 JPERM(I) = 0_8 813 PR(I) = IP(I) 814 D(I) = RINF 815 10 CONTINUE 816 DO 15 I = 1,M 817 U(I) = RINF3 818 IPERM(I) = 0 819 L(I) = 0_8 820 15 CONTINUE 821 DO 30 J = 1,N 822 IF (int(IP(J+1)-IP(J)) .GT. N/10 .AND. N.GT.50) GO TO 30 823 DO 20 K = IP(J),IP(J+1)-1 824 I = IRN(K) 825 IF (A(K).GT.U(I)) GO TO 20 826 U(I) = A(K) 827 IPERM(I) = J 828 L(I) = K 829 20 CONTINUE 830 30 CONTINUE 831 DO 40 I = 1,M 832 J = IPERM(I) 833 IF (J.EQ.0) GO TO 40 834 IF (JPERM(J).EQ.0_8) THEN 835 JPERM(J) = L(I) 836 D(J) = U(I) 837 NUM = NUM + 1 838 ELSEIF (D(J).GT.U(I)) THEN 839 K = JPERM(J) 840 II = IRN(K) 841 IPERM(II) = 0 842 JPERM(J) = L(I) 843 D(J) = U(I) 844 ELSE 845 IPERM(I) = 0 846 ENDIF 847 40 CONTINUE 848 IF (NUM.EQ.N) GO TO 1000 849 DO 45 I = 1,M 850 D(I) = ZERO 851 45 CONTINUE 852 DO 95 J = 1,N 853 IF (JPERM(J).NE.0) GO TO 95 854 K1 = IP(J) 855 K2 = IP(J+1) - 1_8 856 IF (K1.GT.K2) GO TO 95 857 VJ = RINF 858 DO 50 K = K1,K2 859 I = IRN(K) 860 DI = A(K) - U(I) 861 IF (DI.GT.VJ) GO TO 50 862 IF (DI.LT.VJ .OR. DI.EQ.RINF) GO TO 55 863 IF (IPERM(I).NE.0 .OR. IPERM(I0).EQ.0) GO TO 50 864 55 VJ = DI 865 I0 = I 866 K0 = K 867 50 CONTINUE 868 D(J) = VJ 869 K = K0 870 I = I0 871 IF (IPERM(I).EQ.0) GO TO 90 872 DO 60 K = K0,K2 873 I = IRN(K) 874 IF (A(K)-U(I).GT.VJ) GO TO 60 875 JJ = IPERM(I) 876 KK1 = PR(JJ) 877 KK2 = IP(JJ+1) - 1_8 878 IF (KK1.GT.KK2) GO TO 60 879 DO 70 KK = KK1,KK2 880 II = IRN(KK) 881 IF (IPERM(II).GT.0) GO TO 70 882 IF (A(KK)-U(II).LE.D(JJ)) GO TO 80 883 70 CONTINUE 884 PR(JJ) = KK2 + 1_8 885 60 CONTINUE 886 GO TO 95 887 80 JPERM(JJ) = KK 888 IPERM(II) = JJ 889 PR(JJ) = KK + 1_8 890 90 NUM = NUM + 1 891 JPERM(J) = K 892 IPERM(I) = J 893 PR(J) = K + 1_8 894 95 CONTINUE 895 IF (NUM.EQ.N) GO TO 1000 896 DO 99 I = 1,M 897 D(I) = RINF 898 Q(I) = 0 899 99 CONTINUE 900 DO 100 JORD = 1,N 901 IF (JPERM(JORD).NE.0) GO TO 100 902 DMIN = RINF 903 QLEN = 0 904 LOW = M + 1 905 UP = M + 1 906 CSP = RINF 907 J = JORD 908 PR(J) = -1_8 909 DO 115 K = IP(J),IP(J+1)-1_8 910 I = IRN(K) 911 DNEW = A(K) - U(I) 912 IF (DNEW.GE.CSP) GO TO 115 913 IF (IPERM(I).EQ.0) THEN 914 CSP = DNEW 915 ISP = K 916 JSP = J 917 ELSE 918 IF (DNEW.LT.DMIN) DMIN = DNEW 919 D(I) = DNEW 920 QLEN = QLEN + 1 921 L(QLEN) = K 922 ENDIF 923 115 CONTINUE 924 Q0 = QLEN 925 QLEN = 0 926 DO 120 IK = 1,Q0 927 K = L(IK) 928 I = IRN(K) 929 IF (CSP.LE.D(I)) THEN 930 D(I) = RINF 931 GO TO 120 932 ENDIF 933 IF (D(I).LE.DMIN) THEN 934 LOW = LOW - 1 935 L32(LOW) = I 936 Q(I) = LOW 937 ELSE 938 QLEN = QLEN + 1 939 Q(I) = QLEN 940 CALL ZMUMPS_MTRANSD(I,M,L32,D,Q,2) 941 ENDIF 942 JJ = IPERM(I) 943 OUT(JJ) = K 944 PR(JJ) = int(J,8) 945 120 CONTINUE 946 DO 150 JDUM = 1,NUM 947 IF (LOW.EQ.UP) THEN 948 IF (QLEN.EQ.0) GO TO 160 949 I = L32(1) 950 IF (D(I).LT.RINF) DMIN = D(I)*(ONE+RLX) 951 IF (DMIN.GE.CSP) GO TO 160 952 152 CALL ZMUMPS_MTRANSE(QLEN,M,L32,D,Q,2) 953 LOW = LOW - 1 954 L32(LOW) = I 955 Q(I) = LOW 956 IF (QLEN.EQ.0) GO TO 153 957 I = L32(1) 958 IF (D(I).GT.DMIN) GO TO 153 959 GO TO 152 960 ENDIF 961 153 Q0 = L32(UP-1) 962 DQ0 = D(Q0) 963 IF (DQ0.GE.CSP) GO TO 160 964 IF (DMIN.GE.CSP) GO TO 160 965 UP = UP - 1 966 J = IPERM(Q0) 967 VJ = DQ0 - A(JPERM(J)) + U(Q0) 968 K1 = IP(J+1)-1_8 969 IF (LORD) THEN 970 IF (CSP.NE.RINF) THEN 971 DI = CSP - VJ 972 IF (A(K1).GE.DI) THEN 973 K0 = JPERM(J) 974 IF (K0.GE.K1-6) GO TO 178 975 177 CONTINUE 976 K = (K0+K1)/2 977 IF (A(K).GE.DI) THEN 978 K1 = K 979 ELSE 980 K0 = K 981 ENDIF 982 IF (K0.GE.K1-6) GO TO 178 983 GO TO 177 984 178 DO 179 K = K0+1,K1 985 IF (A(K).LT.DI) GO TO 179 986 K1 = K - 1 987 GO TO 181 988 179 CONTINUE 989 ENDIF 990 ENDIF 991 181 IF (K1.EQ.JPERM(J)) K1 = K1 - 1 992 ENDIF 993 K0 = IP(J) 994 DI = CSP - VJ 995 DO 155 K = K0,K1 996 I = IRN(K) 997 IF (Q(I).GE.LOW) GO TO 155 998 DNEW = A(K) - U(I) 999 IF (DNEW.GE.DI) GO TO 155 1000 DNEW = DNEW + VJ 1001 IF (DNEW.GT.D(I)) GO TO 155 1002 IF (IPERM(I).EQ.0) THEN 1003 CSP = DNEW 1004 ISP = K 1005 JSP = J 1006 DI = CSP - VJ 1007 ELSE 1008 IF (DNEW.GE.D(I)) GO TO 155 1009 D(I) = DNEW 1010 IF (DNEW.LE.DMIN) THEN 1011 IF (Q(I).NE.0) THEN 1012 CALL ZMUMPS_MTRANSF(Q(I),QLEN,M,L32,D,Q,2) 1013 ENDIF 1014 LOW = LOW - 1 1015 L32(LOW) = I 1016 Q(I) = LOW 1017 ELSE 1018 IF (Q(I).EQ.0) THEN 1019 QLEN = QLEN + 1 1020 Q(I) = QLEN 1021 ENDIF 1022 CALL ZMUMPS_MTRANSD(I,M,L32,D,Q,2) 1023 ENDIF 1024 JJ = IPERM(I) 1025 OUT(JJ) = K 1026 PR(JJ) = int(J,8) 1027 ENDIF 1028 155 CONTINUE 1029 150 CONTINUE 1030 160 IF (CSP.EQ.RINF) GO TO 190 1031 NUM = NUM + 1 1032 I = IRN(ISP) 1033 J = JSP 1034 IPERM(I) = J 1035 JPERM(J) = ISP 1036 DO 170 JDUM = 1,NUM 1037 JJ = int(PR(J)) 1038 IF (JJ.EQ.-1) GO TO 180 1039 K = OUT(J) 1040 I = IRN(K) 1041 IPERM(I) = JJ 1042 JPERM(JJ) = K 1043 J = JJ 1044 170 CONTINUE 1045 180 DO 182 JJ = UP,M 1046 I = L32(JJ) 1047 U(I) = U(I) + D(I) - CSP 1048 182 CONTINUE 1049 190 DO 191 JJ = UP,M 1050 I = L32(JJ) 1051 D(I) = RINF 1052 Q(I) = 0 1053 191 CONTINUE 1054 DO 192 JJ = LOW,UP-1 1055 I = L32(JJ) 1056 D(I) = RINF 1057 Q(I) = 0 1058 192 CONTINUE 1059 DO 193 JJ = 1,QLEN 1060 I = L32(JJ) 1061 D(I) = RINF 1062 Q(I) = 0 1063 193 CONTINUE 1064 100 CONTINUE 1065 1000 CONTINUE 1066 DO 1200 J = 1,N 1067 K = JPERM(J) 1068 IF (K.NE.0) THEN 1069 D(J) = A(K) - U(IRN(K)) 1070 ELSE 1071 D(J) = ZERO 1072 ENDIF 1073 1200 CONTINUE 1074 DO 1201 I = 1,M 1075 IF (IPERM(I).EQ.0) U(I) = ZERO 1076 1201 CONTINUE 1077 IF (M.EQ.N .and. NUM.EQ.N) GO TO 2000 1078 CALL ZMUMPS_MTRANSX(M,N,IPERM,Q,L32) 1079 2000 RETURN 1080 END SUBROUTINE ZMUMPS_MTRANSW 1081 SUBROUTINE ZMUMPS_MTRANSZ 1082 & (M,N,IRN,LIRN,IP,LENC,IPERM,NUM,PR,ARP,CV,OUT) 1083 IMPLICIT NONE 1084 INTEGER :: M,N,NUM 1085 INTEGER(8), INTENT(IN) :: LIRN 1086 INTEGER :: ARP(N),CV(M),IRN(LIRN),IPERM(M),LENC(N),OUT(N),PR(N) 1087 INTEGER(8), INTENT(IN) :: IP(N) 1088C Local variables 1089 INTEGER :: I,J,J1,JORD,K,KK 1090 INTEGER(8) :: II, IN1, IN2 1091 EXTERNAL ZMUMPS_MTRANSX 1092 DO 10 I = 1,M 1093 CV(I) = 0 1094 IPERM(I) = 0 1095 10 CONTINUE 1096 DO 12 J = 1,N 1097 ARP(J) = LENC(J) - 1 1098 12 CONTINUE 1099 NUM = 0 1100 DO 1000 JORD = 1,N 1101 J = JORD 1102 PR(J) = -1 1103 DO 70 K = 1,JORD 1104 IN1 = int(ARP(J),8) 1105 IF (IN1.LT.0_8) GO TO 30 1106 IN2 = IP(J) + int(LENC(J) - 1,8) 1107 IN1 = IN2 - IN1 1108 DO 20 II = IN1,IN2 1109 I = IRN(II) 1110 IF (IPERM(I).EQ.0) GO TO 80 1111 20 CONTINUE 1112 ARP(J) = -1 1113 30 CONTINUE 1114 OUT(J) = LENC(J) - 1 1115 DO 60 KK = 1,JORD 1116 IN1 = int(OUT(J),8) 1117 IF (IN1.LT.0_8) GO TO 50 1118 IN2 = IP(J) + int(LENC(J) - 1,8) 1119 IN1 = IN2 - IN1 1120 DO 40 II = IN1,IN2 1121 I = IRN(II) 1122 IF (CV(I).EQ.JORD) GO TO 40 1123 J1 = J 1124 J = IPERM(I) 1125 CV(I) = JORD 1126 PR(J) = J1 1127 OUT(J1) = int(IN2 - II - 1_8) 1128 GO TO 70 1129 40 CONTINUE 1130 50 CONTINUE 1131 J = PR(J) 1132 IF (J.EQ.-1) GO TO 1000 1133 60 CONTINUE 1134 70 CONTINUE 1135 80 CONTINUE 1136 IPERM(I) = J 1137 ARP(J) = int(IN2 - II - 1_8) 1138 NUM = NUM + 1 1139 DO 90 K = 1,JORD 1140 J = PR(J) 1141 IF (J.EQ.-1) GO TO 1000 1142 II = IP(J) + int(LENC(J) - OUT(J) - 2,8) 1143 I = IRN(II) 1144 IPERM(I) = J 1145 90 CONTINUE 1146 1000 CONTINUE 1147 IF (M.EQ.N .and. NUM.EQ.N) GO TO 2000 1148 CALL ZMUMPS_MTRANSX(M,N,IPERM,CV,ARP) 1149 2000 RETURN 1150 END SUBROUTINE ZMUMPS_MTRANSZ 1151 SUBROUTINE ZMUMPS_MTRANSX(M,N,IPERM,RW,CW) 1152 IMPLICIT NONE 1153 INTEGER M,N 1154 INTEGER RW(M),CW(N),IPERM(M) 1155 INTEGER I,J,K 1156 DO 10 J = 1,N 1157 CW(J) = 0 1158 10 CONTINUE 1159 K = 0 1160 DO 20 I = 1,M 1161 IF (IPERM(I).EQ.0) THEN 1162 K = K + 1 1163 RW(K) = I 1164 ELSE 1165 J = IPERM(I) 1166 CW(J) = I 1167 ENDIF 1168 20 CONTINUE 1169 K = 0 1170 DO 30 J = 1,N 1171 IF (CW(J).NE.0) GO TO 30 1172 K = K + 1 1173 I = RW(K) 1174 IPERM(I) = -J 1175 30 CONTINUE 1176 DO 40 J = N+1,M 1177 K = K + 1 1178 I = RW(K) 1179 IPERM(I) = -J 1180 40 CONTINUE 1181 RETURN 1182 END SUBROUTINE ZMUMPS_MTRANSX 1183