1C 2C This file is part of MUMPS 4.10.0, built on Tue May 10 12:56:32 UTC 2011 3C 4C 5C This version of MUMPS is provided to you free of charge. It is public 6C domain, based on public domain software developed during the Esprit IV 7C European project PARASOL (1996-1999). Since this first public domain 8C version in 1999, research and developments have been supported by the 9C following institutions: CERFACS, CNRS, ENS Lyon, INPT(ENSEEIHT)-IRIT, 10C INRIA, and University of Bordeaux. 11C 12C The MUMPS team at the moment of releasing this version includes 13C Patrick Amestoy, Maurice Bremond, Alfredo Buttari, Abdou Guermouche, 14C Guillaume Joslin, Jean-Yves L'Excellent, Francois-Henry Rouet, Bora 15C Ucar and Clement Weisbecker. 16C 17C We are also grateful to Emmanuel Agullo, Caroline Bousquet, Indranil 18C Chowdhury, Philippe Combes, Christophe Daniel, Iain Duff, Vincent Espirat, 19C Aurelia Fevre, Jacko Koster, Stephane Pralet, Chiara Puglisi, Gregoire 20C Richard, Tzvetomila Slavova, Miroslav Tuma and Christophe Voemel who 21C have been contributing to this project. 22C 23C Up-to-date copies of the MUMPS package can be obtained 24C from the Web pages: 25C http://mumps.enseeiht.fr/ or http://graal.ens-lyon.fr/MUMPS 26C 27C 28C THIS MATERIAL IS PROVIDED AS IS, WITH ABSOLUTELY NO WARRANTY 29C EXPRESSED OR IMPLIED. ANY USE IS AT YOUR OWN RISK. 30C 31C 32C User documentation of any code that uses this software can 33C include this complete notice. You can acknowledge (using 34C references [1] and [2]) the contribution of this package 35C in any scientific publication dependent upon the use of the 36C package. You shall use reasonable endeavours to notify 37C the authors of the package of this publication. 38C 39C [1] P. R. Amestoy, I. S. Duff, J. Koster and J.-Y. L'Excellent, 40C A fully asynchronous multifrontal solver using distributed dynamic 41C scheduling, SIAM Journal of Matrix Analysis and Applications, 42C Vol 23, No 1, pp 15-41 (2001). 43C 44C [2] P. R. Amestoy and A. Guermouche and J.-Y. L'Excellent and 45C S. Pralet, Hybrid scheduling for the parallel solution of linear 46C systems. Parallel Computing Vol 32 (2), pp 136-156 (2006). 47C 48 SUBROUTINE SMUMPS_635(N,KEEP,ICNTL,MPG) 49 IMPLICIT NONE 50 INTEGER N, KEEP(500), ICNTL(40), MPG 51 KEEP(19)=0 52 RETURN 53 END SUBROUTINE SMUMPS_635 54 SUBROUTINE SMUMPS_634(ICNTL,KEEP,MPG,INFO) 55 IMPLICIT NONE 56 INTEGER KEEP(500), MPG, INFO(40), ICNTL(40) 57 IF (KEEP(19).EQ.0.AND.KEEP(110).EQ.0) THEN 58 IF (KEEP(111).NE.0) THEN 59 INFO(1) = -37 60 INFO(2) = 16 61 IF (KEEP(110).EQ.0) INFO(2) = 24 62 IF(MPG.GT.0) THEN 63 WRITE( MPG,'(A)') 64 &'** ERROR : Null space computation requirement' 65 WRITE( MPG,'(A)') 66 &'** not consistent with factorization options' 67 ENDIF 68 GOTO 333 69 ENDIF 70 ENDIF 71 IF (ICNTL(9).NE.1) THEN 72 IF (KEEP(111).NE.0) THEN 73 INFO(1) = -37 74 INFO(2) = 9 75 IF (MPG.GT.0) THEN 76 WRITE(MPG,'(A)') 77 &'** ERROR ICNTL(25) incompatible with ' 78 WRITE( MPG,'(A)') 79 &'** option transposed system (ICNLT(9)=1) ' 80 ENDIF 81 ENDIF 82 GOTO 333 83 ENDIF 84 333 CONTINUE 85 RETURN 86 END SUBROUTINE SMUMPS_634 87 SUBROUTINE SMUMPS_637(id) 88 USE SMUMPS_STRUC_DEF 89 IMPLICIT NONE 90 TYPE (SMUMPS_STRUC) id 91 NULLIFY(id%root%QR_TAU) 92 RETURN 93 END SUBROUTINE SMUMPS_637 94 SUBROUTINE SMUMPS_636(id) 95 USE SMUMPS_STRUC_DEF 96 IMPLICIT NONE 97 TYPE (SMUMPS_STRUC) id 98 IF (associated(id%root%QR_TAU)) THEN 99 DEALLOCATE(id%root%QR_TAU) 100 NULLIFY(id%root%QR_TAU) 101 ENDIF 102 RETURN 103 END SUBROUTINE SMUMPS_636 104 SUBROUTINE SMUMPS_146( MYID, root, N, IROOT, 105 & COMM, IW, LIW, IFREE, 106 & A, LA, PTRAST, PTLUST_S, PTRFAC, 107 & STEP, INFO, LDLT, QR, 108 & WK, LWK, KEEP,KEEP8,DKEEP) 109 IMPLICIT NONE 110 INCLUDE 'smumps_root.h' 111 INCLUDE 'mpif.h' 112 TYPE ( SMUMPS_ROOT_STRUC ) :: root 113 INTEGER N, IROOT, COMM, LIW, MYID, IFREE 114 INTEGER(8) :: LA 115 INTEGER(8) :: LWK 116 REAL WK( LWK ) 117 INTEGER KEEP(500) 118 REAL DKEEP(30) 119 INTEGER(8) KEEP8(150) 120 INTEGER(8) :: PTRAST(KEEP(28)) 121 INTEGER(8) :: PTRFAC(KEEP(28)) 122 INTEGER PTLUST_S(KEEP(28)), STEP(N), IW( LIW ) 123 INTEGER INFO( 2 ), LDLT, QR 124 REAL A( LA ) 125 INTEGER IOLDPS 126 INTEGER(8) :: IAPOS 127 INTEGER LOCAL_M, LOCAL_N, LPIV, IERR, allocok 128 INTEGER FWD_LOCAL_N_RHS, FWD_MTYPE 129 INCLUDE 'mumps_headers.h' 130 EXTERNAL numroc 131 INTEGER numroc 132 IF ( .NOT. root%yes ) RETURN 133 IF ( KEEP(60) .NE. 0 ) THEN 134 IF ((LDLT == 1 .OR. LDLT == 2) .AND. KEEP(60) == 3 ) THEN 135 CALL SMUMPS_320( WK, root%MBLOCK, 136 & root%MYROW, root%MYCOL, root%NPROW, root%NPCOL, 137 & root%SCHUR_POINTER(1), 138 & root%SCHUR_LLD, root%SCHUR_NLOC, 139 & root%TOT_ROOT_SIZE, MYID, COMM ) 140 ENDIF 141 RETURN 142 ENDIF 143 IOLDPS = PTLUST_S(STEP(IROOT))+KEEP(IXSZ) 144 IAPOS = PTRAST(STEP(IROOT)) 145 LOCAL_M = IW( IOLDPS + 2 ) 146 LOCAL_N = IW( IOLDPS + 1 ) 147 IAPOS = PTRFAC(IW ( IOLDPS + 4 )) 148 IF ( LDLT.EQ.0 .OR. LDLT.EQ.2 .OR. QR.ne.0 ) THEN 149 LPIV = LOCAL_M + root%MBLOCK 150 ELSE 151 LPIV = 1 152 END IF 153 IF (associated( root%IPIV )) DEALLOCATE(root%IPIV) 154 root%LPIV = LPIV 155 ALLOCATE( root%IPIV( LPIV ), stat = allocok ) 156 IF ( allocok .GT. 0 ) THEN 157 INFO(1) = -13 158 INFO(2) = LPIV 159 WRITE(*,*) MYID,': problem allocating IPIV(',LPIV,') in root' 160 CALL MUMPS_ABORT() 161 END IF 162 CALL DESCINIT( root%DESCRIPTOR(1), root%TOT_ROOT_SIZE, 163 & root%TOT_ROOT_SIZE, root%MBLOCK, root%NBLOCK, 164 & 0, 0, root%CNTXT_BLACS, LOCAL_M, IERR ) 165 IF ( LDLT.EQ.2 ) THEN 166 IF(root%MBLOCK.NE.root%NBLOCK) THEN 167 WRITE(*,*) ' Error: symmetrization only works for' 168 WRITE(*,*) ' square block sizes, MBLOCK/NBLOCK=', 169 & root%MBLOCK, root%NBLOCK 170 CALL MUMPS_ABORT() 171 END IF 172 IF ( LWK .LT. min( 173 & int(root%MBLOCK,8) * int(root%NBLOCK,8), 174 & int(root%TOT_ROOT_SIZE,8)* int(root%TOT_ROOT_SIZE,8 ) 175 & )) THEN 176 WRITE(*,*) 'Not enough workspace for symmetrization.' 177 CALL MUMPS_ABORT() 178 END IF 179 CALL SMUMPS_320( WK, root%MBLOCK, 180 & root%MYROW, root%MYCOL, root%NPROW, root%NPCOL, 181 & A( IAPOS ), LOCAL_M, LOCAL_N, 182 & root%TOT_ROOT_SIZE, MYID, COMM ) 183 END IF 184 IF (LDLT.EQ.0.OR.LDLT.EQ.2) THEN 185 CALL psgetrf( root%TOT_ROOT_SIZE, root%TOT_ROOT_SIZE, 186 & A( IAPOS ), 187 & 1, 1, root%DESCRIPTOR(1), root%IPIV(1), IERR ) 188 IF ( IERR .GT. 0 ) THEN 189 INFO(1)=-10 190 INFO(2)=IERR-1 191 END IF 192 ELSE 193 CALL pspotrf('L',root%TOT_ROOT_SIZE,A(IAPOS), 194 & 1,1,root%DESCRIPTOR(1),IERR) 195 IF ( IERR .GT. 0 ) THEN 196 INFO(1)=-40 197 INFO(2)=IERR-1 198 END IF 199 END IF 200 IF (KEEP(258).NE.0) THEN 201 IF (root%MBLOCK.NE.root%NBLOCK) THEN 202 write(*,*) "Internal error in SMUMPS_146:", 203 & "Block size different for rows and columns", 204 & root%MBLOCK, root%NBLOCK 205 CALL MUMPS_ABORT() 206 ENDIF 207 CALL SMUMPS_763(root%MBLOCK, root%IPIV(1),root%MYROW, 208 & root%MYCOL, root%NPROW, root%NPCOL, A(IAPOS), LOCAL_M, 209 & LOCAL_N, root%TOT_ROOT_SIZE, MYID, DKEEP(6), KEEP(259), 210 & LDLT) 211 ENDIF 212 IF (KEEP(252) .NE. 0) THEN 213 FWD_LOCAL_N_RHS = numroc(KEEP(253), root%NBLOCK, 214 & root%MYCOL, 0, root%NPCOL) 215 FWD_LOCAL_N_RHS = max(1,FWD_LOCAL_N_RHS) 216 FWD_MTYPE = 1 217 CALL SMUMPS_768( 218 & root%TOT_ROOT_SIZE, 219 & KEEP(253), 220 & FWD_MTYPE, 221 & A(IAPOS), 222 & root%DESCRIPTOR(1), 223 & LOCAL_M, LOCAL_N, FWD_LOCAL_N_RHS, 224 & root%IPIV(1), LPIV, 225 & root%RHS_ROOT(1,1), LDLT, 226 & root%MBLOCK, root%NBLOCK, 227 & root%CNTXT_BLACS, IERR) 228 ENDIF 229 RETURN 230 END SUBROUTINE SMUMPS_146 231 SUBROUTINE SMUMPS_556( 232 & N,PIV,FRERE,FILS,NFSIZ,IKEEP, 233 & NCST,KEEP,KEEP8,id) 234 USE SMUMPS_STRUC_DEF 235 IMPLICIT NONE 236 TYPE (SMUMPS_STRUC) :: id 237 INTEGER N,NCST 238 INTEGER PIV(N),FRERE(N),FILS(N),NFSIZ(N),IKEEP(N,3) 239 INTEGER KEEP(500) 240 INTEGER(8) KEEP8(150) 241 INTEGER I,P11,P1,P2,K1,K2,NLOCKED 242 LOGICAL V1,V2 243 NCST = 0 244 NLOCKED = 0 245 P11 = KEEP(93) 246 DO I=KEEP(93)-1,1,-2 247 P1 = PIV(I) 248 P2 = PIV(I+1) 249 K1 = IKEEP(P1,1) 250 IF(K1 .GT. 0) THEN 251 V1 = (abs(id%A(K1))*(id%ROWSCA(P1)**2).GE.1.0E-1) 252 ELSE 253 V1 = .FALSE. 254 ENDIF 255 K2 = IKEEP(P2,1) 256 IF(K2 .GT. 0) THEN 257 V2 = (abs(id%A(K2))*(id%ROWSCA(P2)**2).GE.1.0E-1) 258 ELSE 259 V2 = .FALSE. 260 ENDIF 261 IF(V1 .AND. V2) THEN 262 PIV(P11) = P1 263 P11 = P11 - 1 264 PIV(P11) = P2 265 P11 = P11 - 1 266 ELSE IF(V1) THEN 267 NCST = NCST+1 268 FRERE(NCST) = P1 269 NCST = NCST+1 270 FRERE(NCST) = P2 271 ELSE IF(V2) THEN 272 NCST = NCST+1 273 FRERE(NCST) = P2 274 NCST = NCST+1 275 FRERE(NCST) = P1 276 ELSE 277 NLOCKED = NLOCKED + 1 278 FILS(NLOCKED) = P1 279 NLOCKED = NLOCKED + 1 280 FILS(NLOCKED) = P2 281 ENDIF 282 ENDDO 283 DO I=1,NLOCKED 284 PIV(I) = FILS(I) 285 ENDDO 286 KEEP(94) = KEEP(94) + KEEP(93) - NLOCKED 287 KEEP(93) = NLOCKED 288 DO I=1,NCST 289 NLOCKED = NLOCKED + 1 290 PIV(NLOCKED) = FRERE(I) 291 ENDDO 292 DO I=1,KEEP(93)/2 293 NFSIZ(I) = 0 294 ENDDO 295 DO I=(KEEP(93)/2)+1,(KEEP(93)/2)+NCST,2 296 NFSIZ(I) = I+1 297 NFSIZ(I+1) = -1 298 ENDDO 299 DO I=(KEEP(93)/2)+NCST+1,(KEEP(93)/2)+KEEP(94) 300 NFSIZ(I) = 0 301 ENDDO 302 END SUBROUTINE SMUMPS_556 303 SUBROUTINE SMUMPS_550(N,NCMP,N11,N22,PIV, 304 & INVPERM,PERM) 305 IMPLICIT NONE 306 INTEGER N11,N22,N,NCMP 307 INTEGER, intent(in) :: PIV(N),PERM(N) 308 INTEGER, intent (out):: INVPERM(N) 309 INTEGER CMP_POS,EXP_POS,I,J,N2,K 310 N2 = N22/2 311 EXP_POS = 1 312 DO CMP_POS=1,NCMP 313 J = PERM(CMP_POS) 314 IF(J .LE. N2) THEN 315 K = 2*J-1 316 I = PIV(K) 317 INVPERM(I) = EXP_POS 318 EXP_POS = EXP_POS+1 319 K = K+1 320 I = PIV(K) 321 INVPERM(I) = EXP_POS 322 EXP_POS = EXP_POS+1 323 ELSE 324 K = N2 + J 325 I = PIV(K) 326 INVPERM(I) = EXP_POS 327 EXP_POS = EXP_POS+1 328 ENDIF 329 ENDDO 330 DO K=N22+N11+1,N 331 I = PIV(K) 332 INVPERM(I) = EXP_POS 333 EXP_POS = EXP_POS+1 334 ENDDO 335 RETURN 336 END SUBROUTINE SMUMPS_550 337 SUBROUTINE SMUMPS_547( 338 & N,NZ, IRN, ICN, PIV, 339 & NCMP, IW, LW, IPE, LEN, IQ, 340 & FLAG, ICMP, IWFR, 341 & IERROR, KEEP,KEEP8, ICNTL) 342 IMPLICIT NONE 343 INTEGER N,NZ,NCMP,LW,IWFR,IERROR 344 INTEGER ICNTL(40),KEEP(500) 345 INTEGER(8) KEEP8(150) 346 INTEGER IRN(NZ),ICN(NZ),IW(LW),PIV(N),IPE(N+1) 347 INTEGER LEN(N),IQ(N),FLAG(N),ICMP(N) 348 INTEGER MP,N11,N22,NDUP 349 INTEGER I,K,J,N1,LAST,K1,K2,L 350 MP = ICNTL(2) 351 IERROR = 0 352 N22 = KEEP(93) 353 N11 = KEEP(94) 354 NCMP = N22/2 + N11 355 DO I=1,NCMP 356 IPE(I) = 0 357 ENDDO 358 K = 1 359 DO I=1,N22/2 360 J = PIV(K) 361 ICMP(J) = I 362 K = K + 1 363 J = PIV(K) 364 ICMP(J) = I 365 K = K + 1 366 ENDDO 367 K = N22/2 + 1 368 DO I=N22+1,N22+N11 369 J = PIV(I) 370 ICMP(J) = K 371 K = K + 1 372 ENDDO 373 DO I=N11+N22+1,N 374 J = PIV(I) 375 ICMP(J) = 0 376 ENDDO 377 DO K=1,NZ 378 I = IRN(K) 379 J = ICN(K) 380 I = ICMP(I) 381 J = ICMP(J) 382 IF ((I.GT.N).OR.(J.GT.N).OR.(I.LT.1) 383 & .OR.(J.LT.1)) THEN 384 IERROR = IERROR + 1 385 ELSE 386 IF (I.NE.J) THEN 387 IPE(I) = IPE(I) + 1 388 IPE(J) = IPE(J) + 1 389 ENDIF 390 ENDIF 391 ENDDO 392 IQ(1) = 1 393 N1 = NCMP - 1 394 IF (N1.GT.0) THEN 395 DO I=1,N1 396 IQ(I+1) = IPE(I) + IQ(I) 397 ENDDO 398 ENDIF 399 LAST = max(IPE(NCMP)+IQ(NCMP)-1,IQ(NCMP)) 400 DO I = 1,NCMP 401 FLAG(I) = 0 402 IPE(I) = IQ(I) 403 ENDDO 404 DO K=1,LAST 405 IW(K) = 0 406 ENDDO 407 IWFR = LAST + 1 408 DO K=1,NZ 409 I = IRN(K) 410 J = ICN(K) 411 I = ICMP(I) 412 J = ICMP(J) 413 IF (I.NE.J) THEN 414 IF (I.LT.J) THEN 415 IF ((I.GE.1).AND.(J.LE.N)) THEN 416 IW(IQ(I)) = -J 417 IQ(I) = IQ(I) + 1 418 ENDIF 419 ELSE 420 IF ((J.GE.1).AND.(I.LE.N)) THEN 421 IW(IQ(J)) = -I 422 IQ(J) = IQ(J) + 1 423 ENDIF 424 ENDIF 425 ENDIF 426 ENDDO 427 NDUP = 0 428 DO I=1,NCMP 429 K1 = IPE(I) 430 K2 = IQ(I) -1 431 IF (K1.GT.K2) THEN 432 LEN(I) = 0 433 IQ(I) = 0 434 ELSE 435 DO K=K1,K2 436 J = -IW(K) 437 IF (J.LE.0) GO TO 250 438 L = IQ(J) 439 IQ(J) = L + 1 440 IF (FLAG(J).EQ.I) THEN 441 NDUP = NDUP + 1 442 IW(L) = 0 443 IW(K) = 0 444 ELSE 445 IW(L) = I 446 IW(K) = J 447 FLAG(J) = I 448 ENDIF 449 ENDDO 450 250 IQ(I) = IQ(I) - IPE(I) 451 IF (NDUP.EQ.0) LEN(I) = IQ(I) 452 ENDIF 453 ENDDO 454 IF (NDUP.NE.0) THEN 455 IWFR = 1 456 DO I=1,NCMP 457 K1 = IPE(I) 458 IF (IQ(I).EQ.0) THEN 459 LEN(I) = 0 460 IPE(I) = IWFR 461 CYCLE 462 ENDIF 463 K2 = K1 + IQ(I) - 1 464 L = IWFR 465 IPE(I) = IWFR 466 DO K=K1,K2 467 IF (IW(K).NE.0) THEN 468 IW(IWFR) = IW(K) 469 IWFR = IWFR + 1 470 ENDIF 471 ENDDO 472 LEN(I) = IWFR - L 473 ENDDO 474 ENDIF 475 IPE(NCMP+1) = IPE(NCMP) + LEN(NCMP) 476 IWFR = IPE(NCMP+1) 477 RETURN 478 END SUBROUTINE SMUMPS_547 479 SUBROUTINE SMUMPS_551( 480 & N, NE, IP, IRN, SCALING,LSC,CPERM, DIAG, 481 & ICNTL, WEIGHT,MARKED,FLAG, 482 & PIV_OUT, INFO) 483 IMPLICIT NONE 484 INTEGER N, NE, ICNTL(10), INFO(10),LSC 485 INTEGER CPERM(N),PIV_OUT(N),IP(N+1),IRN(NE), DIAG(N) 486 REAL SCALING(LSC),WEIGHT(N+2) 487 INTEGER MARKED(N),FLAG(N) 488 INTEGER NUM1,NUM2,NUMTOT,PATH_LENGTH,NLAST 489 INTEGER I,CUR_EL,CUR_EL_PATH,CUR_EL_PATH_NEXT,BEST_BEG 490 INTEGER L1,L2,PTR_SET1,PTR_SET2,TUP,T22 491 REAL BEST_SCORE,CUR_VAL,TMP,VAL 492 REAL INITSCORE, SMUMPS_739, 493 & SMUMPS_740, SMUMPS_741 494 LOGICAL VRAI,FAUX,MAX_CARD_DIAG,USE_SCALING 495 INTEGER SUM 496 REAL ZERO,ONE 497 PARAMETER (SUM = 1, VRAI = .TRUE., FAUX = .FALSE.) 498 PARAMETER(ZERO = 0.0E0, ONE = 1.0E0) 499 MAX_CARD_DIAG = .TRUE. 500 NUM1 = 0 501 NUM2 = 0 502 NUMTOT = 0 503 NLAST = N 504 INFO = 0 505 MARKED = 1 506 FLAG = 0 507 VAL = ONE 508 IF(LSC .GT. 1) THEN 509 USE_SCALING = .TRUE. 510 ELSE 511 USE_SCALING = .FALSE. 512 ENDIF 513 TUP = ICNTL(2) 514 IF(TUP .EQ. SUM) THEN 515 INITSCORE = ZERO 516 ELSE 517 INITSCORE = ONE 518 ENDIF 519 IF(ICNTL(2) .GT. 2 .OR. ICNTL(2) .LE. 0) THEN 520 WRITE(*,*) 521 & 'ERROR: WRONG VALUE FOR ICNTL(2) = ',ICNTL(2) 522 INFO(1) = -1 523 RETURN 524 ENDIF 525 T22 = ICNTL(1) 526 IF(ICNTL(1) .LT. 0 .OR. ICNTL(1) .GT. 2) THEN 527 WRITE(*,*) 528 & 'ERROR: WRONG VALUE FOR ICNTL(1) = ',ICNTL(1) 529 INFO(1) = -1 530 RETURN 531 ENDIF 532 DO CUR_EL=1,N 533 IF(MARKED(CUR_EL) .LE. 0) THEN 534 CYCLE 535 ENDIF 536 IF(CPERM(CUR_EL) .LT. 0) THEN 537 MARKED(CUR_EL) = -1 538 CYCLE 539 ENDIF 540 PATH_LENGTH = 2 541 CUR_EL_PATH = CPERM(CUR_EL) 542 IF(CUR_EL_PATH .EQ. CUR_EL) THEN 543 MARKED(CUR_EL) = -1 544 CYCLE 545 ENDIF 546 MARKED(CUR_EL) = 0 547 WEIGHT(1) = INITSCORE 548 WEIGHT(2) = INITSCORE 549 L1 = IP(CUR_EL+1)-IP(CUR_EL) 550 L2 = IP(CUR_EL_PATH+1)-IP(CUR_EL_PATH) 551 PTR_SET1 = IP(CUR_EL) 552 PTR_SET2 = IP(CUR_EL_PATH) 553 IF(USE_SCALING) THEN 554 VAL = -SCALING(CUR_EL_PATH) - SCALING(CUR_EL+N) 555 ENDIF 556 CUR_VAL = SMUMPS_741( 557 & CUR_EL,CUR_EL_PATH, 558 & IRN(PTR_SET1),IRN(PTR_SET2), 559 & L1,L2, 560 & VAL,DIAG,N,FLAG,FAUX,T22) 561 WEIGHT(PATH_LENGTH+1) = 562 & SMUMPS_739(WEIGHT(1),CUR_VAL,TUP) 563 DO 564 IF(CUR_EL_PATH .EQ. CUR_EL) EXIT 565 PATH_LENGTH = PATH_LENGTH+1 566 MARKED(CUR_EL_PATH) = 0 567 CUR_EL_PATH_NEXT = CPERM(CUR_EL_PATH) 568 L1 = IP(CUR_EL_PATH+1)-IP(CUR_EL_PATH) 569 L2 = IP(CUR_EL_PATH_NEXT+1)-IP(CUR_EL_PATH_NEXT) 570 PTR_SET1 = IP(CUR_EL_PATH) 571 PTR_SET2 = IP(CUR_EL_PATH_NEXT) 572 IF(USE_SCALING) THEN 573 VAL = -SCALING(CUR_EL_PATH_NEXT) 574 & - SCALING(CUR_EL_PATH+N) 575 ENDIF 576 CUR_VAL = SMUMPS_741( 577 & CUR_EL_PATH,CUR_EL_PATH_NEXT, 578 & IRN(PTR_SET1),IRN(PTR_SET2), 579 & L1,L2, 580 & VAL,DIAG,N,FLAG,VRAI,T22) 581 WEIGHT(PATH_LENGTH+1) = 582 & SMUMPS_739(WEIGHT(PATH_LENGTH-1),CUR_VAL,TUP) 583 CUR_EL_PATH = CUR_EL_PATH_NEXT 584 ENDDO 585 IF(mod(PATH_LENGTH,2) .EQ. 1) THEN 586 IF(WEIGHT(PATH_LENGTH+1) .GE. WEIGHT(PATH_LENGTH)) THEN 587 CUR_EL_PATH = CPERM(CUR_EL) 588 ELSE 589 CUR_EL_PATH = CUR_EL 590 ENDIF 591 DO I=1,(PATH_LENGTH-1)/2 592 NUM2 = NUM2+1 593 PIV_OUT(NUM2) = CUR_EL_PATH 594 CUR_EL_PATH = CPERM(CUR_EL_PATH) 595 NUM2 = NUM2+1 596 PIV_OUT(NUM2) = CUR_EL_PATH 597 CUR_EL_PATH = CPERM(CUR_EL_PATH) 598 ENDDO 599 NUMTOT = NUMTOT + PATH_LENGTH - 1 600 ELSE 601 IF(MAX_CARD_DIAG) THEN 602 CUR_EL_PATH = CPERM(CUR_EL) 603 IF(DIAG(CUR_EL) .NE. 0) THEN 604 BEST_BEG = CUR_EL_PATH 605 GOTO 1000 606 ENDIF 607 DO I=1,(PATH_LENGTH/2) 608 CUR_EL_PATH_NEXT = CPERM(CUR_EL_PATH) 609 IF(DIAG(CUR_EL_PATH) .NE. 0) THEN 610 BEST_BEG = CUR_EL_PATH_NEXT 611 GOTO 1000 612 ENDIF 613 ENDDO 614 ENDIF 615 BEST_BEG = CUR_EL 616 BEST_SCORE = WEIGHT(PATH_LENGTH-1) 617 CUR_EL_PATH = CPERM(CUR_EL) 618 DO I=1,(PATH_LENGTH/2)-1 619 TMP = SMUMPS_739(WEIGHT(PATH_LENGTH), 620 & WEIGHT(2*I-1),TUP) 621 TMP = SMUMPS_740(TMP,WEIGHT(2*I),TUP) 622 IF(TMP .GT. BEST_SCORE) THEN 623 BEST_SCORE = TMP 624 BEST_BEG = CUR_EL_PATH 625 ENDIF 626 CUR_EL_PATH = CPERM(CUR_EL_PATH) 627 TMP = SMUMPS_739(WEIGHT(PATH_LENGTH+1), 628 & WEIGHT(2*I),TUP) 629 TMP = SMUMPS_740(TMP,WEIGHT(2*I+1),TUP) 630 IF(TMP .GT. BEST_SCORE) THEN 631 BEST_SCORE = TMP 632 BEST_BEG = CUR_EL_PATH 633 ENDIF 634 CUR_EL_PATH = CPERM(CUR_EL_PATH) 635 ENDDO 636 1000 CUR_EL_PATH = BEST_BEG 637 DO I=1,(PATH_LENGTH/2)-1 638 NUM2 = NUM2+1 639 PIV_OUT(NUM2) = CUR_EL_PATH 640 CUR_EL_PATH = CPERM(CUR_EL_PATH) 641 NUM2 = NUM2+1 642 PIV_OUT(NUM2) = CUR_EL_PATH 643 CUR_EL_PATH = CPERM(CUR_EL_PATH) 644 ENDDO 645 NUMTOT = NUMTOT + PATH_LENGTH - 2 646 MARKED(CUR_EL_PATH) = -1 647 ENDIF 648 ENDDO 649 DO I=1,N 650 IF(MARKED(I) .LT. 0) THEN 651 IF(DIAG(I) .EQ. 0) THEN 652 PIV_OUT(NLAST) = I 653 NLAST = NLAST - 1 654 ELSE 655 NUM1 = NUM1 + 1 656 PIV_OUT(NUM2+NUM1) = I 657 NUMTOT = NUMTOT + 1 658 ENDIF 659 ENDIF 660 ENDDO 661 INFO(2) = NUMTOT 662 INFO(3) = NUM1 663 INFO(4) = NUM2 664 RETURN 665 END SUBROUTINE SMUMPS_551 666 FUNCTION SMUMPS_739(A,B,T) 667 IMPLICIT NONE 668 REAL SMUMPS_739 669 REAL A,B 670 INTEGER T 671 INTEGER SUM 672 PARAMETER(SUM = 1) 673 IF(T .EQ. SUM) THEN 674 SMUMPS_739 = A+B 675 ELSE 676 SMUMPS_739 = A*B 677 ENDIF 678 END FUNCTION SMUMPS_739 679 FUNCTION SMUMPS_740(A,B,T) 680 IMPLICIT NONE 681 REAL SMUMPS_740 682 REAL A,B 683 INTEGER T 684 INTEGER SUM 685 PARAMETER(SUM = 1) 686 IF(T .EQ. SUM) THEN 687 SMUMPS_740 = A-B 688 ELSE 689 SMUMPS_740 = A/B 690 ENDIF 691 END FUNCTION SMUMPS_740 692 FUNCTION SMUMPS_741(CUR_EL,CUR_EL_PATH, 693 & SET1,SET2,L1,L2,VAL,DIAG,N,FLAG,FLAGON,T) 694 IMPLICIT NONE 695 REAL SMUMPS_741 696 INTEGER CUR_EL,CUR_EL_PATH,L1,L2,N 697 INTEGER SET1(L1),SET2(L2),DIAG(N),FLAG(N) 698 REAL VAL 699 LOGICAL FLAGON 700 INTEGER T 701 INTEGER I,INTER,MERGE 702 INTEGER STRUCT,MA47 703 PARAMETER(STRUCT=0,MA47=1) 704 IF(T .EQ. STRUCT) THEN 705 IF(.NOT. FLAGON) THEN 706 DO I=1,L1 707 FLAG(SET1(I)) = CUR_EL 708 ENDDO 709 ENDIF 710 INTER = 0 711 DO I=1,L2 712 IF(FLAG(SET2(I)) .EQ. CUR_EL) THEN 713 INTER = INTER + 1 714 FLAG(SET2(I)) = CUR_EL_PATH 715 ENDIF 716 ENDDO 717 MERGE = L1 + L2 - INTER 718 SMUMPS_741 = real(INTER) / real(MERGE) 719 ELSE IF (T .EQ. MA47) THEN 720 MERGE = 3 721 IF(DIAG(CUR_EL) .NE. 0) MERGE = 2 722 IF(DIAG(CUR_EL_PATH) .NE. 0) MERGE = MERGE - 2 723 IF(MERGE .EQ. 0) THEN 724 SMUMPS_741 = real(L1+L2-2) 725 SMUMPS_741 = -(SMUMPS_741**2)/2.0E0 726 ELSE IF(MERGE .EQ. 1) THEN 727 SMUMPS_741 = - real(L1+L2-4) * real(L1-2) 728 ELSE IF(MERGE .EQ. 2) THEN 729 SMUMPS_741 = - real(L1+L2-4) * real(L2-2) 730 ELSE 731 SMUMPS_741 = - real(L1-2) * real(L2-2) 732 ENDIF 733 ELSE 734 SMUMPS_741 = VAL 735 ENDIF 736 RETURN 737 END FUNCTION 738 SUBROUTINE SMUMPS_622(NA, NCMP, 739 & INVPERM,PERM, 740 & LISTVAR_SCHUR, SIZE_SCHUR, AOTOA) 741 IMPLICIT NONE 742 INTEGER, INTENT(IN):: SIZE_SCHUR, LISTVAR_SCHUR(SIZE_SCHUR) 743 INTEGER, INTENT(IN):: NA, NCMP 744 INTEGER, INTENT(IN):: AOTOA(NCMP), PERM(NCMP) 745 INTEGER, INTENT(OUT):: INVPERM(NA) 746 INTEGER CMP_POS, IO, I, K, IPOS 747 DO CMP_POS=1, NCMP 748 IO = PERM(CMP_POS) 749 INVPERM(AOTOA(IO)) = CMP_POS 750 ENDDO 751 IPOS = NCMP 752 DO K =1, SIZE_SCHUR 753 I = LISTVAR_SCHUR(K) 754 IPOS = IPOS+1 755 INVPERM(I) = IPOS 756 ENDDO 757 RETURN 758 END SUBROUTINE SMUMPS_622 759 SUBROUTINE SMUMPS_623 760 & (NA,N,NZ, IRN, ICN, IW, LW, IPE, LEN, 761 & IQ, FLAG, IWFR, 762 & NRORM, NIORM, IFLAG,IERROR, ICNTL, 763 & symmetry, SYM, MedDens, NBQD, AvgDens, 764 & LISTVAR_SCHUR, SIZE_SCHUR, ATOAO, AOTOA) 765 IMPLICIT NONE 766 INTEGER, INTENT(IN) :: NA,N,NZ,LW 767 INTEGER, INTENT(IN) :: SIZE_SCHUR, LISTVAR_SCHUR(SIZE_SCHUR) 768 INTEGER, INTENT(IN) :: IRN(NZ), ICN(NZ) 769 INTEGER, INTENT(IN) :: ICNTL(40), SYM 770 INTEGER, INTENT(INOUT) :: IFLAG 771 INTEGER, INTENT(OUT) :: IERROR,NRORM,NIORM,IWFR 772 INTEGER, INTENT(OUT) :: AOTOA(N) 773 INTEGER, INTENT(OUT) :: ATOAO(NA) 774 INTEGER, INTENT(OUT) :: LEN(N), IPE(N+1) 775 INTEGER, INTENT(OUT) :: symmetry, 776 & MedDens, NBQD, AvgDens 777 INTEGER, INTENT(OUT) :: FLAG(N), IW(LW), IQ(N) 778 INTEGER MP, MPG 779 INTEGER I,K,J,N1,LAST,NDUP,K1,K2,L 780 INTEGER NBERR, THRESH, IAO 781 INTEGER NZOFFA, NDIAGA 782 REAL RSYM 783 INTRINSIC nint 784 ATOAO(1:NA) = 0 785 DO I = 1, SIZE_SCHUR 786 ATOAO(LISTVAR_SCHUR(I)) = -1 787 ENDDO 788 IAO = 0 789 DO I= 1, NA 790 IF (ATOAO(I).LT.0) CYCLE 791 IAO = IAO +1 792 ATOAO(I) = IAO 793 AOTOA(IAO) = I 794 ENDDO 795 MP = ICNTL(2) 796 MPG= ICNTL(3) 797 NIORM = 3*N 798 NDIAGA = 0 799 IERROR = 0 800 IPE(1:N+1) = 0 801 DO K=1,NZ 802 I = IRN(K) 803 J = ICN(K) 804 IF ((I.GT.NA).OR.(J.GT.NA).OR.(I.LT.1) 805 & .OR.(J.LT.1)) THEN 806 IERROR = IERROR + 1 807 ELSE 808 I = ATOAO(I) 809 J = ATOAO(J) 810 IF ((I.LT.0).OR.(J.LT.0)) CYCLE 811 IF (I.NE.J) THEN 812 IPE(I) = IPE(I) + 1 813 IPE(J) = IPE(J) + 1 814 NIORM = NIORM + 1 815 ELSE 816 NDIAGA = NDIAGA + 1 817 ENDIF 818 ENDIF 819 ENDDO 820 NZOFFA = NIORM - 3*N 821 IF (IERROR.GE.1) THEN 822 NBERR = 0 823 IF (mod(IFLAG,2).EQ.0) IFLAG = IFLAG+1 824 IF ((MP.GT.0).AND.(ICNTL(4).GE.2)) THEN 825 WRITE (MP,99999) 826 DO 70 K=1,NZ 827 I = IRN(K) 828 J = ICN(K) 829 IF ((I.GT.NA).OR.(J.GT.NA).OR.(I.LT.1) 830 & .OR.(J.LT.1)) THEN 831 NBERR = NBERR + 1 832 IF (NBERR.LE.10) THEN 833 IF (mod(K,10).GT.3 .OR. mod(K,10).EQ.0 .OR. 834 & (10.LE.K .AND. K.LE.20)) THEN 835 WRITE (MP,'(I8,A,I8,A,I8,A)') 836 & K,'th entry (in row',I,' and column',J,') ignored' 837 ELSE 838 IF (mod(K,10).EQ.1) WRITE(MP,'(I8,A,I8,A,I8,A)') 839 & K,'st entry (in row',I,' and column',J,') ignored' 840 IF (mod(K,10).EQ.2) WRITE(MP,'(I8,A,I8,A,I8,A)') 841 & K,'nd entry (in row',I,' and column',J,') ignored' 842 IF (mod(K,10).EQ.3) WRITE(MP,'(I8,A,I8,A,I8,A)') 843 & K,'rd entry (in row',I,' and column',J,') ignored' 844 ENDIF 845 ELSE 846 GO TO 100 847 ENDIF 848 ENDIF 849 70 CONTINUE 850 ENDIF 851 ENDIF 852 100 NRORM = NIORM - 2*N 853 IQ(1) = 1 854 N1 = N - 1 855 IF (N1.GT.0) THEN 856 DO 110 I=1,N1 857 IQ(I+1) = IPE(I) + IQ(I) 858 110 CONTINUE 859 ENDIF 860 LAST = max(IPE(N)+IQ(N)-1,IQ(N)) 861 FLAG(1:N) = 0 862 IPE(1:N) = IQ(1:N) 863 IW(1:LAST) = 0 864 IWFR = LAST + 1 865 DO 200 K=1,NZ 866 I = IRN(K) 867 J = ICN(K) 868 IF ((I.GT.NA).OR.(J.GT.NA).OR.(I.LT.1) 869 & .OR.(J.LT.1)) CYCLE 870 I = ATOAO(I) 871 J = ATOAO(J) 872 IF ((I.LT.0).OR.(J.LT.0)) CYCLE 873 IF (I.NE.J) THEN 874 IF (I.LT.J) THEN 875 IW(IQ(I)) = -J 876 IQ(I) = IQ(I) + 1 877 ELSE 878 IW(IQ(J)) = -I 879 IQ(J) = IQ(J) + 1 880 ENDIF 881 ENDIF 882 200 CONTINUE 883 NDUP = 0 884 DO 260 I=1,N 885 K1 = IPE(I) 886 K2 = IQ(I) -1 887 IF (K1.GT.K2) THEN 888 LEN(I) = 0 889 IQ(I) = 0 890 ELSE 891 DO 240 K=K1,K2 892 J = -IW(K) 893 IF (J.LE.0) GO TO 250 894 L = IQ(J) 895 IQ(J) = L + 1 896 IF (FLAG(J).EQ.I) THEN 897 NDUP = NDUP + 1 898 IW(L) = 0 899 IW(K) = 0 900 ELSE 901 IW(L) = I 902 IW(K) = J 903 FLAG(J) = I 904 ENDIF 905 240 CONTINUE 906 250 IQ(I) = IQ(I) - IPE(I) 907 IF (NDUP.EQ.0) LEN(I) = IQ(I) 908 ENDIF 909 260 CONTINUE 910 IF (NDUP.NE.0) THEN 911 IWFR = 1 912 DO 280 I=1,N 913 IF (IQ(I).EQ.0) THEN 914 LEN(I) = 0 915 IPE(I) = IWFR 916 GOTO 280 917 ENDIF 918 K1 = IPE(I) 919 K2 = K1 + IQ(I) - 1 920 L = IWFR 921 IPE(I) = IWFR 922 DO 270 K=K1,K2 923 IF (IW(K).NE.0) THEN 924 IW(IWFR) = IW(K) 925 IWFR = IWFR + 1 926 ENDIF 927 270 CONTINUE 928 LEN(I) = IWFR - L 929 280 CONTINUE 930 ENDIF 931 IPE(N+1) = IPE(N) + LEN(N) 932 IWFR = IPE(N+1) 933 IF (SYM.EQ.0) THEN 934 RSYM = real(NDIAGA+2*NZOFFA - (IWFR-1))/ 935 & real(NZOFFA+NDIAGA) 936 symmetry = nint (100.0E0*RSYM) 937 IF (MPG .GT. 0) 938 & write(MPG,'(A,I5)') 939 & ' ... Structural symmetry (in percent)=', symmetry 940 IF (MP.GT.0 .AND. MPG.NE.MP) 941 & write(MP,'(A,I5)') 942 & ' ... Structural symmetry (in percent)=', symmetry 943 ELSE 944 symmetry = 100 945 ENDIF 946 AvgDens = nint(real(IWFR-1)/real(N)) 947 THRESH = AvgDens*50 - AvgDens/10 + 1 948 NBQD = 0 949 IF (N.GT.2) THEN 950 IQ(1:N) = 0 951 DO I= 1, N 952 K = max(LEN(I),1) 953 IQ(K) = IQ(K) + 1 954 IF (K.GT.THRESH) NBQD = NBQD+1 955 ENDDO 956 K = 0 957 MedDens = 0 958 DO WHILE (K .LT. (N/2)) 959 MedDens = MedDens + 1 960 K = K+IQ(MedDens) 961 ENDDO 962 ELSE 963 MedDens = AvgDens 964 ENDIF 965 IF (MPG .GT. 0) 966 & write(MPG,'(A,3I5)') 967 & ' Density: NBdense, Average, Median =', 968 & NBQD, AvgDens, MedDens 969 IF (MP.GT.0 .AND. MPG.NE.MP) 970 & write(MP,'(A,3I5)') 971 & ' Density: NBdense, Average, Median =', 972 & NBQD, AvgDens, MedDens 973 RETURN 97499999 FORMAT (/'*** Warning message from analysis routine ***') 975 END SUBROUTINE SMUMPS_623 976 SUBROUTINE SMUMPS_549(N,PE,INVPERM,NFILS,WORK) 977 IMPLICIT NONE 978 INTEGER N 979 INTEGER PE(N),INVPERM(N),NFILS(N),WORK(N) 980 INTEGER I,FATHER,STKLEN,STKPOS,PERMPOS,CURVAR 981 NFILS = 0 982 DO I=1,N 983 FATHER = -PE(I) 984 IF(FATHER .NE. 0) NFILS(FATHER) = NFILS(FATHER) + 1 985 ENDDO 986 STKLEN = 0 987 PERMPOS = 1 988 DO I=1,N 989 IF(NFILS(I) .EQ. 0) THEN 990 STKLEN = STKLEN + 1 991 WORK(STKLEN) = I 992 INVPERM(I) = PERMPOS 993 PERMPOS = PERMPOS + 1 994 ENDIF 995 ENDDO 996 DO STKPOS = 1,STKLEN 997 CURVAR = WORK(STKPOS) 998 FATHER = -PE(CURVAR) 999 DO 1000 IF(FATHER .EQ. 0) EXIT 1001 IF(NFILS(FATHER) .EQ. 1) THEN 1002 INVPERM(FATHER) = PERMPOS 1003 FATHER = -PE(FATHER) 1004 PERMPOS = PERMPOS + 1 1005 ELSE 1006 NFILS(FATHER) = NFILS(FATHER) - 1 1007 EXIT 1008 ENDIF 1009 ENDDO 1010 ENDDO 1011 RETURN 1012 END SUBROUTINE SMUMPS_549 1013 SUBROUTINE SMUMPS_548(N,PE,NV,WORK) 1014 IMPLICIT NONE 1015 INTEGER N 1016 INTEGER PE(N),NV(N),WORK(N) 1017 INTEGER I,FATHER,LEN,NEWSON,NEWFATHER 1018 DO I=1,N 1019 IF(NV(I) .GT. 0) CYCLE 1020 LEN = 1 1021 WORK(LEN) = I 1022 FATHER = -PE(I) 1023 DO 1024 IF(NV(FATHER) .GT. 0) THEN 1025 NEWSON = FATHER 1026 EXIT 1027 ENDIF 1028 LEN = LEN + 1 1029 WORK(LEN) = FATHER 1030 NV(FATHER) = 1 1031 FATHER = -PE(FATHER) 1032 ENDDO 1033 NEWFATHER = -PE(FATHER) 1034 PE(WORK(LEN)) = -NEWFATHER 1035 PE(NEWSON) = -WORK(1) 1036 ENDDO 1037 END SUBROUTINE SMUMPS_548 1038