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 MODULE ZMUMPS_FAC_FRONT_AUX_M 14 CONTAINS 15 SUBROUTINE ZMUMPS_FAC_H(NFRONT,NASS,IW,LIW,A,LA, 16 & INOPV,NOFFW,IOLDPS,POSELT,UU,SEUIL,KEEP, DKEEP, 17 & PP_FIRST2SWAP_L, PP_LastPanelonDisk_L, 18 & PP_LastPIVRPTRFilled_L, 19 & PP_FIRST2SWAP_U, PP_LastPanelonDisk_U, 20 & PP_LastPIVRPTRFilled_U,MAXFROMN,IS_MAXFROMN_AVAIL, 21 & Inextpiv 22 &) 23!$ USE OMP_LIB 24 USE MUMPS_OOC_COMMON 25 IMPLICIT NONE 26 INTEGER NFRONT,NASS,LIW,INOPV 27 INTEGER(8) :: LA 28 INTEGER KEEP(500) 29 DOUBLE PRECISION DKEEP(230) 30 DOUBLE PRECISION UU, SEUIL 31 COMPLEX(kind=8) A(LA) 32 INTEGER IW(LIW) 33 DOUBLE PRECISION, intent(in) :: MAXFROMN 34 LOGICAL, intent(inout) :: IS_MAXFROMN_AVAIL 35 INTEGER, intent(inout) :: Inextpiv 36 DOUBLE PRECISION AMROW 37 DOUBLE PRECISION RMAX 38 COMPLEX(kind=8) SWOP 39 INTEGER(8) :: APOS, POSELT 40 INTEGER(8) :: J1, J2, J3_8, JJ, IDIAG 41 INTEGER(8) :: J1_ini 42 INTEGER(8) :: NFRONT8 43 INTEGER IOLDPS 44 INTEGER NOFFW,NPIV,IPIV,IPIV_SHIFT 45 INTEGER J, J3 46 INTEGER NPIVP1,JMAX,ISW,ISWPS1 47 INTEGER ISWPS2,KSW,XSIZE 48 INTEGER I_PIVRPTR_L, I_PIVR_L, NBPANELS_L 49 INTEGER I_PIVRPTR_U, I_PIVR_U, NBPANELS_U 50 INTEGER PP_FIRST2SWAP_L, PP_LastPanelonDisk_L, 51 & PP_LastPIVRPTRFilled_L, 52 & PP_FIRST2SWAP_U, PP_LastPanelonDisk_U, 53 & PP_LastPIVRPTRFilled_U 54 INTEGER ISHIFT, K206 55 INTEGER ZMUMPS_IXAMAX 56 INCLUDE 'mumps_headers.h' 57 INTRINSIC max 58 DOUBLE PRECISION, PARAMETER :: RZERO = 0.0D0 59!$ INTEGER :: JJMAX 60!$ DOUBLE PRECISION :: RRMAX, VALABS 61!$ INTEGER :: NOMP, CHUNK, K360 62!$ K360 = KEEP(360) 63!$ NOMP = OMP_GET_MAX_THREADS() 64 NFRONT8 = int(NFRONT,8) 65 INOPV = 0 66 XSIZE = KEEP(IXSZ) 67 NPIV = IW(IOLDPS+1+XSIZE) 68 NPIVP1 = NPIV + 1 69 K206 = KEEP(206) 70 IF (KEEP(201).EQ.1 .AND. KEEP(50).NE.1) THEN 71 CALL ZMUMPS_GET_OOC_PERM_PTR(TYPEF_L, NBPANELS_L, 72 & I_PIVRPTR_L, I_PIVR_L, 73 & IOLDPS+2*NFRONT+6+IW(IOLDPS+5+XSIZE) 74 & +KEEP(IXSZ), 75 & IW, LIW) 76 CALL ZMUMPS_GET_OOC_PERM_PTR(TYPEF_U, NBPANELS_U, 77 & I_PIVRPTR_U, I_PIVR_U, 78 & IOLDPS+2*NFRONT+6+IW(IOLDPS+5+XSIZE)+XSIZE, 79 & IW, LIW) 80 ENDIF 81 ISHIFT = 0 82 IF (K206.GE.1) THEN 83 IF (Inextpiv.GT.NPIVP1.AND.Inextpiv.LE.NASS) THEN 84 ISHIFT = Inextpiv - NPIVP1 85 ENDIF 86 IF (ISHIFT.GT.0.AND.IS_MAXFROMN_AVAIL) THEN 87 IPIV = NPIVP1 88 APOS = POSELT + NFRONT8*int(NPIV,8) + int(IPIV-1,8) 89 IDIAG = APOS + int(IPIV - NPIVP1,8)*NFRONT8 90 IF (abs(A(IDIAG)).GT.max(UU*MAXFROMN,SEUIL, 91 & tiny(MAXFROMN))) THEN 92 ISHIFT = 0 93 ENDIF 94 ENDIF 95 IF ( ISHIFT .GT. 0) THEN 96 IS_MAXFROMN_AVAIL = .FALSE. 97 ENDIF 98 ENDIF 99 DO 460 IPIV_SHIFT=NPIVP1+ISHIFT,NASS+ISHIFT 100 IF (IPIV_SHIFT .LE. NASS) THEN 101 IPIV=IPIV_SHIFT 102 ELSE 103 IPIV=IPIV_SHIFT-NASS-1+NPIVP1 104 ENDIF 105 APOS = POSELT + NFRONT8*int(NPIV,8) + int(IPIV-1,8) 106 JMAX = 1 107 AMROW = RZERO 108 J1 = APOS 109 J3 = NASS -NPIV 110 IF (KEEP(351).EQ.1) THEN 111!$ IF (NOMP.GT.1 .AND. J3.GE.K360) THEN 112!$ JMAX = 1 113!$ RMAX = RZERO 114!$ CHUNK = max(K360/2,J3/NOMP) 115!$OMP PARALLEL PRIVATE(JJ,VALABS,JJMAX,RRMAX) 116!$OMP& FIRSTPRIVATE(J1,NFRONT8,J3) 117!$ RRMAX = RZERO 118!$OMP DO schedule(static, CHUNK) 119!$ DO J = 1, J3 120!$ JJ = J1 + int(J-1,8)*NFRONT8 121!$ VALABS = abs(A(JJ)) 122!$ IF (VALABS.GT.RRMAX) THEN 123!$ RRMAX = VALABS 124!$ JJMAX = J 125!$ ENDIF 126!$ ENDDO 127!$OMP END DO 128!$ IF (RRMAX.GT.0.0) THEN 129!$OMP CRITICAL 130!$ IF (RRMAX.GT.RMAX) THEN 131!$ RMAX = RRMAX 132!$ JMAX = JJMAX 133!$ ENDIF 134!$OMP END CRITICAL 135!$ ENDIF 136!$OMP END PARALLEL 137!$ ELSE 138 JMAX = ZMUMPS_IXAMAX(J3,A(J1),NFRONT) 139!$ ENDIF 140 ELSE 141 JMAX = ZMUMPS_IXAMAX(J3,A(J1),NFRONT) 142 ENDIF 143 JJ = J1 + int(JMAX-1,8)*NFRONT8 144 AMROW = abs(A(JJ)) 145 RMAX = AMROW 146 J1 = APOS + int(NASS-NPIV,8) * NFRONT8 147 J3 = NFRONT - NASS - KEEP(253) 148 IF (IS_MAXFROMN_AVAIL) THEN 149 RMAX = max(MAXFROMN,RMAX) 150 IS_MAXFROMN_AVAIL = .FALSE. 151 ELSE 152 IF (J3.EQ.0) GOTO 370 153 IF (KEEP(351).EQ.1) THEN 154 J1_ini = J1 155!$ CHUNK = max(K360/2,J3/NOMP) 156!$OMP PARALLEL DO schedule(static, CHUNK) PRIVATE(J1) 157!$OMP& FIRSTPRIVATE(J1_ini,NFRONT8,J3) 158!$OMP& REDUCTION(max:RMAX) IF (J3.GE.K360) 159 DO J=1,J3 160 J1 = J1_ini + int(J-1,8) * NFRONT8 161 RMAX = max(abs(A(J1)),RMAX) 162 END DO 163!$OMP END PARALLEL DO 164 ELSE 165 DO J=1,J3 166 RMAX = max(abs(A(J1)), RMAX) 167 J1 = J1 + NFRONT8 168 END DO 169 ENDIF 170 END IF 171 370 IF (RMAX.LE.tiny(RMAX)) GO TO 460 172 IDIAG = APOS + int(IPIV - NPIVP1,8)*NFRONT8 173 IF (abs(A(IDIAG)).GT.max(UU*RMAX,SEUIL,tiny(RMAX))) THEN 174 JMAX = IPIV - NPIV 175 GO TO 380 176 ENDIF 177 IF (AMROW.LE.max(UU*RMAX,SEUIL,tiny(RMAX))) GO TO 460 178 NOFFW = NOFFW + 1 179 380 CONTINUE 180 IF (K206.GE.1) THEN 181 Inextpiv = IPIV + 1 182 ENDIF 183 IF (KEEP(258) .NE. 0) THEN 184 CALL ZMUMPS_UPDATEDETER( 185 & A(APOS + int(JMAX - 1,8) * NFRONT8 ), 186 & DKEEP(6), 187 & KEEP(259) ) 188 ENDIF 189 IF (IPIV.EQ.NPIVP1) GO TO 400 190 KEEP(260)=-KEEP(260) 191 J1 = POSELT + int(NPIV,8) 192 J3_8 = POSELT + int(IPIV-1,8) 193 DO J= 1,NFRONT 194 SWOP = A(J1) 195 A(J1) = A(J3_8) 196 A(J3_8) = SWOP 197 J1 = J1 + NFRONT8 198 J3_8 = J3_8 + NFRONT8 199 END DO 200 ISWPS1 = IOLDPS + 5 + NPIVP1 + NFRONT + XSIZE 201 ISWPS2 = IOLDPS + 5 + IPIV + NFRONT + XSIZE 202 ISW = IW(ISWPS1) 203 IW(ISWPS1) = IW(ISWPS2) 204 IW(ISWPS2) = ISW 205 400 IF (JMAX.EQ.1) GO TO 420 206 KEEP(260)=-KEEP(260) 207 J1 = POSELT + int(NPIV,8) * NFRONT8 208 J2 = POSELT + int(NPIV + JMAX - 1,8) * NFRONT8 209 DO KSW=1,NFRONT 210 SWOP = A(J1) 211 A(J1) = A(J2) 212 A(J2) = SWOP 213 J1 = J1 + 1_8 214 J2 = J2 + 1_8 215 END DO 216 ISWPS1 = IOLDPS + 5 + NPIV + 1 + XSIZE 217 ISWPS2 = IOLDPS + 5 + NPIV + JMAX + XSIZE 218 ISW = IW(ISWPS1) 219 IW(ISWPS1) = IW(ISWPS2) 220 IW(ISWPS2) = ISW 221 GO TO 420 222 460 CONTINUE 223 INOPV = 1 224 GOTO 430 225 420 CONTINUE 226 IF (KEEP(201).EQ.1) THEN 227 IF (KEEP(251).EQ.0) THEN 228 CALL ZMUMPS_STORE_PERMINFO( IW(I_PIVRPTR_L), 229 & NBPANELS_L, 230 & IW(I_PIVR_L), NASS, NPIVP1, NPIV+JMAX, 231 & PP_LastPanelonDisk_L, 232 & PP_LastPIVRPTRFilled_L) 233 ENDIF 234 CALL ZMUMPS_STORE_PERMINFO( IW(I_PIVRPTR_U), 235 & NBPANELS_U, 236 & IW(I_PIVR_U), NASS, NPIVP1, IPIV, 237 & PP_LastPanelonDisk_U, 238 & PP_LastPIVRPTRFilled_U) 239 ENDIF 240 430 CONTINUE 241 IS_MAXFROMN_AVAIL = .FALSE. 242 RETURN 243 END SUBROUTINE ZMUMPS_FAC_H 244 SUBROUTINE ZMUMPS_FAC_M(IBEG_BLOCK, 245 & NFRONT,NASS,N,INODE,IW,LIW,A,LA, 246 & IOLDPS,POSELT,IFINB,LKJIB,LKJIT,XSIZE) 247 IMPLICIT NONE 248 INTEGER NFRONT,NASS,N,LIW,INODE,IFINB,LKJIB,IBEG_BLOCK 249 INTEGER(8) :: LA 250 COMPLEX(kind=8) A(LA) 251 INTEGER IW(LIW) 252 COMPLEX(kind=8) VALPIV 253 INTEGER(8) :: APOS, POSELT, UUPOS, LPOS 254 INTEGER(8) :: NFRONT8 255 INTEGER IOLDPS 256 INTEGER LKJIT, XSIZE 257 COMPLEX(kind=8) ONE, ALPHA 258 INTEGER NPIV,JROW2 259 INTEGER NEL2,NPIVP1,KROW,NEL 260 INCLUDE 'mumps_headers.h' 261 PARAMETER (ONE=(1.0D0,0.0D0), ALPHA=(-1.0D0,0.0D0)) 262 NFRONT8= int(NFRONT,8) 263 NPIV = IW(IOLDPS+1+XSIZE) 264 NPIVP1 = NPIV + 1 265 NEL = NFRONT - NPIVP1 266 IFINB = 0 267 IF (IW(IOLDPS+3+XSIZE).LE.0) THEN 268 IF (NASS.LT.LKJIT) THEN 269 IW(IOLDPS+3+XSIZE) = NASS 270 ELSE 271 IW(IOLDPS+3+XSIZE) = min0(NASS,LKJIB) 272 ENDIF 273 ENDIF 274 JROW2 = IW(IOLDPS+3+XSIZE) 275 NEL2 = JROW2 - NPIVP1 276 IF (NEL2.EQ.0) THEN 277 IF (JROW2.EQ.NASS) THEN 278 IFINB = -1 279 ELSE 280 IFINB = 1 281 IW(IOLDPS+3+XSIZE) = min0(JROW2+LKJIB,NASS) 282 IBEG_BLOCK = NPIVP1+1 283 ENDIF 284 ELSE 285 APOS = POSELT + int(NPIV,8)*(NFRONT8 + 1_8) 286 VALPIV = ONE/A(APOS) 287 LPOS = APOS + NFRONT8 288 DO 541 KROW = 1,NEL2 289 A(LPOS) = A(LPOS)*VALPIV 290 LPOS = LPOS + NFRONT8 291 541 CONTINUE 292 LPOS = APOS + NFRONT8 293 UUPOS = APOS + 1_8 294 CALL zgeru(NEL,NEL2,ALPHA,A(UUPOS),1,A(LPOS),NFRONT, 295 & A(LPOS+1_8),NFRONT) 296 ENDIF 297 RETURN 298 END SUBROUTINE ZMUMPS_FAC_M 299 SUBROUTINE ZMUMPS_FAC_N(NFRONT,NASS,IW,LIW,A,LA, 300 & IOLDPS,POSELT,IFINB,XSIZE, 301 & KEEP,MAXFROMN,IS_MAXFROMN_AVAIL) 302!$ USE OMP_LIB 303 IMPLICIT NONE 304 INCLUDE 'mumps_headers.h' 305 INTEGER NFRONT,NASS,LIW,IFINB 306 INTEGER(8) :: LA 307 COMPLEX(kind=8) A(LA) 308 INTEGER IW(LIW) 309 COMPLEX(kind=8) ALPHA,VALPIV 310 INTEGER(8) :: APOS, POSELT, UUPOS, LPOS, IRWPOS 311 INTEGER(8) :: NFRONT8 312 INTEGER IOLDPS,NPIV,XSIZE 313 INTEGER, intent(in) :: KEEP(500) 314 DOUBLE PRECISION, intent(inout) :: MAXFROMN 315 LOGICAL, intent(inout) :: IS_MAXFROMN_AVAIL 316 INTEGER NEL,IROW,NEL2,JCOL, NCB 317 INTEGER NPIVP1 318 COMPLEX(kind=8), PARAMETER :: ONE=(1.0D0,0.0D0) 319!$ LOGICAL:: OMP_FLAG 320!$ INTEGER:: NOMP, K360, CHUNK 321!$ NOMP = OMP_GET_MAX_THREADS() 322!$ K360 = KEEP(360) 323 NFRONT8=int(NFRONT,8) 324 NPIV = IW(IOLDPS+1+XSIZE) 325 NPIVP1 = NPIV + 1 326 NEL = NFRONT - NPIVP1 327 NEL2 = NASS - NPIVP1 328 NCB = NFRONT - NASS - KEEP(253) 329 IFINB = 0 330 IF (NPIVP1.EQ.NASS) IFINB = 1 331 APOS = POSELT + int(NPIV,8)*(NFRONT8 + 1_8) 332 VALPIV = ONE/A(APOS) 333!$ OMP_FLAG = .FALSE. 334!$ CHUNK = NEL 335!$ IF (NOMP.GT.1) THEN 336!$ IF (NEL.LT.K360) THEN 337!$ IF (NEL*NEL2.GE.KEEP(361)) THEN 338!$ OMP_FLAG = .TRUE. 339!$ CHUNK = max(20, NEL/NOMP) 340!$ ENDIF 341!$ ELSE 342!$ OMP_FLAG = .TRUE. 343!$ CHUNK = max(K360/2,NEL/NOMP) 344!$ ENDIF 345!$ ENDIF 346 IF (KEEP(351).EQ.2) THEN 347 MAXFROMN = 0.0D0 348 IF (NEL2 > 0) THEN 349 IS_MAXFROMN_AVAIL = .TRUE. 350 ENDIF 351!$OMP PARALLEL DO schedule(static, CHUNK) 352!$OMP& PRIVATE(LPOS, UUPOS, IRWPOS, ALPHA, JCOL) 353!$OMP& FIRSTPRIVATE(APOS,NFRONT8,VALPIV,NEL,NEL2) 354!$OMP& REDUCTION(max:MAXFROMN) IF(OMP_FLAG) 355 DO IROW = 1, NEL 356 LPOS = APOS + NFRONT8*int(IROW,8) 357 A(LPOS) = A(LPOS)*VALPIV 358 ALPHA = -A(LPOS) 359 IRWPOS = LPOS + 1_8 360 UUPOS = APOS + 1_8 361 IF (NEL2 > 0) THEN 362 A(IRWPOS) = A(IRWPOS) + ALPHA*A(UUPOS) 363 MAXFROMN=max(MAXFROMN, abs(A(IRWPOS))) 364 IRWPOS = IRWPOS+1_8 365 UUPOS = UUPOS+1_8 366 DO JCOL = 2, NEL2 367 A(IRWPOS) = A(IRWPOS) + ALPHA*A(UUPOS) 368 IRWPOS = IRWPOS+1_8 369 UUPOS = UUPOS+1_8 370 ENDDO 371 ENDIF 372 END DO 373!$OMP END PARALLEL DO 374 ELSE 375!$OMP PARALLEL DO schedule(static, CHUNK) 376!$OMP& FIRSTPRIVATE(APOS,NFRONT8,VALPIV,NEL,NEL2) 377!$OMP& PRIVATE(LPOS, UUPOS, IRWPOS, ALPHA, JCOL) IF(OMP_FLAG) 378 DO IROW = 1, NEL 379 LPOS = APOS + NFRONT8*int(IROW,8) 380 A(LPOS) = A(LPOS)*VALPIV 381 ALPHA = -A(LPOS) 382 IRWPOS = LPOS + 1_8 383 UUPOS = APOS + 1_8 384 DO JCOL = 1, NEL2 385 A(IRWPOS) = A(IRWPOS) + ALPHA*A(UUPOS) 386 IRWPOS = IRWPOS+1_8 387 UUPOS = UUPOS+1_8 388 ENDDO 389 ENDDO 390!$OMP END PARALLEL DO 391 ENDIF 392 RETURN 393 END SUBROUTINE ZMUMPS_FAC_N 394 SUBROUTINE ZMUMPS_FAC_P(A,LA,NFRONT, 395 & NPIV,NASS,POSELT,CALL_UTRSM 396 & ) 397 IMPLICIT NONE 398 INTEGER(8) :: LA,POSELT 399 COMPLEX(kind=8) A(LA) 400 INTEGER NFRONT, NPIV, NASS 401 LOGICAL CALL_UTRSM 402 INTEGER(8) :: LPOS, LPOS1, LPOS2, UPOS 403 INTEGER NEL1,NEL11 404 COMPLEX(kind=8) ALPHA, ONE 405 PARAMETER (ONE=(1.0D0,0.0D0), ALPHA=(-1.0D0,0.0D0)) 406 NEL1 = NFRONT - NASS 407 NEL11 = NFRONT - NPIV 408 LPOS2 = POSELT + int(NASS,8)*int(NFRONT,8) 409 CALL ztrsm('L','L','N','N',NPIV,NEL1,ONE,A(POSELT),NFRONT, 410 & A(LPOS2),NFRONT) 411 IF (CALL_UTRSM) THEN 412 UPOS = POSELT + int(NASS,8) 413 CALL ztrsm('R', 'U', 'N', 'U', NEL1, NPIV, ONE, 414 & A(POSELT), NFRONT, A(UPOS), NFRONT) 415 ENDIF 416 LPOS = LPOS2 + int(NPIV,8) 417 LPOS1 = POSELT + int(NPIV,8) 418 CALL zgemm('N','N',NEL11,NEL1,NPIV,ALPHA,A(LPOS1), 419 & NFRONT,A(LPOS2),NFRONT,ONE,A(LPOS),NFRONT) 420 RETURN 421 END SUBROUTINE ZMUMPS_FAC_P 422 SUBROUTINE ZMUMPS_FAC_P_PANEL(A,LAFAC,NFRONT, 423 & NPIV,NASS, IW, LIWFAC, 424 & MonBloc, TYPEFile, MYID, KEEP8, 425 & STRAT, IFLAG_OOC, 426 & LNextPiv2beWritten, UNextPiv2beWritten) 427 USE ZMUMPS_OOC 428 IMPLICIT NONE 429 INTEGER NFRONT, NPIV, NASS 430 INTEGER(8) :: LAFAC 431 INTEGER LIWFAC, TYPEFile, MYID, IFLAG_OOC, 432 & LNextPiv2beWritten, UNextPiv2beWritten, STRAT 433 COMPLEX(kind=8) A(LAFAC) 434 INTEGER IW(LIWFAC) 435 INTEGER(8) KEEP8(150) 436 TYPE(IO_BLOCK) :: MonBloc 437 INTEGER(8) :: LPOS2,LPOS1,LPOS 438 INTEGER NEL1,NEL11 439 COMPLEX(kind=8) ALPHA, ONE 440 LOGICAL LAST_CALL 441 PARAMETER (ONE=(1.0D0,0.0D0), ALPHA=(-1.0D0,0.0D0)) 442 NEL1 = NFRONT - NASS 443 NEL11 = NFRONT - NPIV 444 LPOS2 = 1_8 + int(NASS,8) * int(NFRONT,8) 445 CALL ztrsm('L','L','N','N',NPIV,NEL1,ONE,A(1),NFRONT, 446 & A(LPOS2),NFRONT) 447 LAST_CALL=.FALSE. 448 CALL ZMUMPS_OOC_IO_LU_PANEL 449 & ( STRAT, TYPEFile, 450 & A, LAFAC, MonBloc, 451 & LNextPiv2beWritten, UNextPiv2beWritten, 452 & IW, LIWFAC, 453 & MYID, KEEP8(31), IFLAG_OOC,LAST_CALL ) 454 LPOS = LPOS2 + int(NPIV,8) 455 LPOS1 = int(1 + NPIV,8) 456 CALL zgemm('N','N',NEL11,NEL1,NPIV,ALPHA,A(LPOS1), 457 & NFRONT,A(LPOS2),NFRONT,ONE,A(LPOS),NFRONT) 458 RETURN 459 END SUBROUTINE ZMUMPS_FAC_P_PANEL 460 SUBROUTINE ZMUMPS_FAC_T(A,LA,NPIVB,NFRONT, 461 & NPIV,NASS,POSELT) 462 IMPLICIT NONE 463 INTEGER NPIVB,NASS 464 INTEGER(8) :: LA 465 COMPLEX(kind=8) A(LA) 466 INTEGER(8) :: APOS, POSELT 467 INTEGER NFRONT, NPIV, NASSL 468 INTEGER(8) :: LPOS, LPOS1, LPOS2 469 INTEGER NEL1, NEL11, NPIVE 470 COMPLEX(kind=8) ALPHA, ONE 471 PARAMETER (ONE=(1.0D0,0.0D0), ALPHA=(-1.0D0,0.0D0)) 472 NEL1 = NFRONT - NASS 473 NEL11 = NFRONT - NPIV 474 NPIVE = NPIV - NPIVB 475 NASSL = NASS - NPIVB 476 APOS = POSELT + int(NPIVB,8)*int(NFRONT,8) 477 & + int(NPIVB,8) 478 LPOS2 = APOS + int(NASSL,8) 479 CALL ztrsm('R','U','N','U',NEL1,NPIVE,ONE,A(APOS),NFRONT, 480 & A(LPOS2),NFRONT) 481 LPOS = LPOS2 + int(NFRONT,8)*int(NPIVE,8) 482 LPOS1 = APOS + int(NFRONT,8)*int(NPIVE,8) 483 CALL zgemm('N','N',NEL1,NEL11,NPIVE,ALPHA,A(LPOS2), 484 & NFRONT,A(LPOS1),NFRONT,ONE,A(LPOS),NFRONT) 485 RETURN 486 END SUBROUTINE ZMUMPS_FAC_T 487 SUBROUTINE ZMUMPS_FAC_SQ(IBEG_BLOCK, IEND_BLOCK, NPIV, 488 & NFRONT, LAST_ROW, LAST_COL, A, LA, POSELT, 489 & CALL_UTRSM, CALL_GEMM, WITH_COMM_THREAD ) 490 IMPLICIT NONE 491 INTEGER, intent(in) :: IBEG_BLOCK, IEND_BLOCK 492 INTEGER, intent(in) :: NPIV, NFRONT, LAST_ROW, LAST_COL 493 INTEGER(8), intent(in) :: LA 494 COMPLEX(kind=8), intent(inout) :: A(LA) 495 INTEGER(8), intent(in) :: POSELT 496 LOGICAL, intent(in) :: CALL_UTRSM, CALL_GEMM 497 LOGICAL, intent(in) :: WITH_COMM_THREAD 498 INTEGER(8) :: NFRONT8 499 INTEGER(8) :: LPOS, LPOS1, LPOS2, UPOS, POSELT_LOCAL 500 INTEGER NELIM, LKJIW, NEL1, NEL11 501 COMPLEX(kind=8) ALPHA, ONE 502 PARAMETER (ONE=(1.0D0,0.0D0), ALPHA=(-1.0D0,0.0D0)) 503 NFRONT8= int(NFRONT,8) 504 NELIM = IEND_BLOCK - NPIV 505 NEL1 = LAST_ROW - IEND_BLOCK 506 IF ( NEL1 < 0 ) THEN 507 WRITE(*,*) 508 & "Internal error 1 in ZMUMPS_FAC_SQ,IEND_BLOCK>LAST_ROW", 509 & IEND_BLOCK, LAST_ROW 510 CALL MUMPS_ABORT() 511 ENDIF 512 LKJIW = NPIV - IBEG_BLOCK + 1 513 NEL11 = LAST_COL - NPIV 514 IF ((NEL1.NE.0).AND.(LKJIW.NE.0)) THEN 515 LPOS2 = POSELT + int(IEND_BLOCK,8)*NFRONT8 + 516 & int(IBEG_BLOCK - 1,8) 517 UPOS = POSELT + int(IBEG_BLOCK-1,8)*NFRONT8 + 518 & int(IEND_BLOCK,8) 519 POSELT_LOCAL = POSELT + 520 & int(IBEG_BLOCK-1,8)*NFRONT8 + int(IBEG_BLOCK - 1,8) 521 CALL ztrsm('L','L','N','N',LKJIW,NEL1,ONE, 522 & A(POSELT_LOCAL),NFRONT, 523 & A(LPOS2),NFRONT) 524 IF (CALL_UTRSM) THEN 525 CALL ztrsm('R','U','N','U',NEL1,LKJIW,ONE, 526 & A(POSELT_LOCAL),NFRONT, 527 & A(UPOS),NFRONT) 528 ENDIF 529 IF (CALL_GEMM) THEN 530 LPOS = LPOS2 + int(LKJIW,8) 531 LPOS1 = POSELT_LOCAL + int(LKJIW,8) 532 CALL zgemm('N','N',NEL11,NEL1,LKJIW,ALPHA,A(LPOS1), 533 & NFRONT,A(LPOS2),NFRONT,ONE,A(LPOS),NFRONT) 534 END IF 535 ENDIF 536 RETURN 537 END SUBROUTINE ZMUMPS_FAC_SQ 538 SUBROUTINE ZMUMPS_FAC_MQ(IBEG_BLOCK,IEND_BLOCK, 539 & NFRONT, NASS, NPIV, LAST_COL, A, LA, POSELT, IFINB) 540 IMPLICIT NONE 541 INTEGER, intent(in) :: IBEG_BLOCK, IEND_BLOCK, NFRONT, 542 & NASS, NPIV, LAST_COL 543 INTEGER, intent(out) :: IFINB 544 INTEGER(8), intent(in) :: LA, POSELT 545 COMPLEX(kind=8), intent(inout) :: A(LA) 546 COMPLEX(kind=8) :: VALPIV 547 INTEGER(8) :: APOS, UUPOS, LPOS 548 INTEGER(8) :: NFRONT8 549 COMPLEX(kind=8) :: ONE, ALPHA 550 INTEGER :: NEL2,NPIVP1,KROW,NEL 551 PARAMETER (ONE=(1.0D0,0.0D0), ALPHA=(-1.0D0,0.0D0)) 552 NFRONT8= int(NFRONT,8) 553 NPIVP1 = NPIV + 1 554 NEL = LAST_COL - NPIVP1 555 IFINB = 0 556 NEL2 = IEND_BLOCK - NPIVP1 557 IF (NEL2.EQ.0) THEN 558 IF (IEND_BLOCK.EQ.NASS) THEN 559 IFINB = -1 560 ELSE 561 IFINB = 1 562 ENDIF 563 ELSE 564 APOS = POSELT + int(NPIV,8)*(NFRONT8 + 1_8) 565 VALPIV = ONE/A(APOS) 566 LPOS = APOS + NFRONT8 567 DO 541 KROW = 1,NEL2 568 A(LPOS) = A(LPOS)*VALPIV 569 LPOS = LPOS + NFRONT8 570 541 CONTINUE 571 LPOS = APOS + NFRONT8 572 UUPOS = APOS + 1_8 573#if defined(MUMPS_USE_BLAS2) 574 CALL zgeru(NEL,NEL2,ALPHA,A(UUPOS),1,A(LPOS),NFRONT, 575 & A(LPOS+1_8),NFRONT) 576#else 577 CALL zgemm('N','N',NEL,NEL2,1,ALPHA,A(UUPOS),NEL, 578 & A(LPOS),NFRONT,ONE,A(LPOS+1_8),NFRONT) 579#endif 580 ENDIF 581 RETURN 582 END SUBROUTINE ZMUMPS_FAC_MQ 583 SUBROUTINE ZMUMPS_FAC_FR_UPDATE_CBROWS( INODE, NFRONT, NASS, 584 & CALL_UTRSM, A, LA, LAFAC, POSELT, IW, LIW, IOLDPS, 585 & MonBloc, MYID, NOFFW, LIWFAC, 586 & PP_FIRST2SWAP_L, PP_FIRST2SWAP_U, 587 & LNextPiv2beWritten, UNextPiv2beWritten, 588 & PP_LastPIVRPTRFilled_L, PP_LastPIVRPTRFilled_U, 589 & 590 & XSIZE, SEUIL, UU, DKEEP, KEEP8, KEEP, IFLAG) 591 USE ZMUMPS_OOC 592 IMPLICIT NONE 593 INTEGER, intent(in) :: INODE, NFRONT, NASS, 594 & LIW, MYID, XSIZE, IOLDPS, LIWFAC 595 INTEGER(8), intent(in) :: LA, POSELT 596 INTEGER, intent(inout) :: NOFFW, 597 & PP_FIRST2SWAP_L, PP_FIRST2SWAP_U, 598 & LNextPiv2beWritten, UNextPiv2beWritten, 599 & PP_LastPIVRPTRFilled_L, PP_LastPIVRPTRFilled_U, 600 & IFLAG 601 LOGICAL, intent(in) :: CALL_UTRSM 602 INTEGER, intent(inout) :: IW(LIW) 603 COMPLEX(kind=8), intent(inout) :: A(LA) 604 DOUBLE PRECISION, intent(in) :: SEUIL, UU, DKEEP(230) 605 INTEGER, intent(in) :: KEEP( 500 ) 606 INTEGER(8), intent(inout) :: LAFAC 607 INTEGER(8) :: KEEP8(150) 608 TYPE(IO_BLOCK), intent(inout) :: MonBloc 609 INTEGER :: NPIV, NEL1, STRAT, TYPEFile, IFLAG_OOC, 610 & IBEG_BLOCK, IFINB, INOPV 611 INTEGER Inextpiv 612 DOUBLE PRECISION :: MAXFROMN 613 LOGICAL :: IS_MAXFROMN_AVAIL 614 NPIV = IW(IOLDPS+1+XSIZE) 615 NEL1 = NFRONT - NASS 616 IF (KEEP(206).GE.1) THEN 617 Inextpiv = 1 618 ELSE 619 Inextpiv = 0 620 ENDIF 621 IF ((NPIV.GT.0).AND.(NEL1.GT.0)) THEN 622 IF (KEEP(201).EQ.1) THEN 623 STRAT = STRAT_TRY_WRITE 624 TYPEFile = TYPEF_BOTH_LU 625 MonBloc%LastPiv= NPIV 626 CALL ZMUMPS_FAC_P_PANEL(A(POSELT), LAFAC, NFRONT, 627 & NPIV, NASS, IW(IOLDPS), LIWFAC, 628 & MonBloc, TYPEFile, MYID, KEEP8, 629 & STRAT, IFLAG_OOC, 630 & LNextPiv2beWritten, UNextPiv2beWritten) 631 IF (IFLAG_OOC < 0 ) IFLAG=IFLAG_OOC 632 ELSE 633 CALL ZMUMPS_FAC_P(A,LA,NFRONT, NPIV, NASS, POSELT, 634 & CALL_UTRSM 635 & ) 636 ENDIF 637 ENDIF 638 NPIV = IW(IOLDPS+1+XSIZE) 639 IBEG_BLOCK = NPIV 640 IF (NASS.EQ.NPIV) GOTO 500 641 IS_MAXFROMN_AVAIL = .FALSE. 642 120 CALL ZMUMPS_FAC_H(NFRONT,NASS,IW,LIW,A,LA, 643 & INOPV,NOFFW,IOLDPS,POSELT,UU,SEUIL, 644 & KEEP, DKEEP, 645 & PP_FIRST2SWAP_L, MonBloc%LastPanelWritten_L, 646 & PP_LastPIVRPTRFilled_L, 647 & PP_FIRST2SWAP_U, MonBloc%LastPanelWritten_U, 648 & PP_LastPIVRPTRFilled_U, MAXFROMN, IS_MAXFROMN_AVAIL, 649 & Inextpiv 650 & ) 651 IF (INOPV.NE.1) THEN 652 CALL ZMUMPS_FAC_N(NFRONT,NASS,IW,LIW,A,LA, 653 & IOLDPS,POSELT,IFINB,XSIZE, 654 & KEEP, MAXFROMN, IS_MAXFROMN_AVAIL) 655 IW(IOLDPS+1+XSIZE) = IW(IOLDPS+1+XSIZE) + 1 656 IF (IFINB.EQ.0) GOTO 120 657 ENDIF 658 NPIV = IW(IOLDPS+1+XSIZE) 659 NEL1 = NFRONT - NASS 660 IF ((NPIV.LE.IBEG_BLOCK).OR.(NEL1.EQ.0)) GO TO 500 661 CALL ZMUMPS_FAC_T(A,LA,IBEG_BLOCK, 662 & NFRONT,NPIV,NASS,POSELT) 663 500 CONTINUE 664 RETURN 665 END SUBROUTINE ZMUMPS_FAC_FR_UPDATE_CBROWS 666 SUBROUTINE ZMUMPS_FAC_I(NFRONT,NASS,LAST_ROW, 667 & IBEG_BLOCK, IEND_BLOCK, 668 & N,INODE,IW,LIW,A,LA, 669 & INOPV,NOFFW,IFLAG,IOLDPS,POSELT,UU,SEUIL,KEEP,KEEP8, 670 & DKEEP,PIVNUL_LIST,LPN_LIST, 671 & 672 & PP_FIRST2SWAP_L, PP_LastPanelonDisk_L, 673 & PP_LastPIVRPTRFilled_L, 674 & PP_FIRST2SWAP_U, PP_LastPanelonDisk_U, 675 & PP_LastPIVRPTRFilled_U, 676 & PIVOT_OPTION, LR_ACTIVATED, IEND_BLR, Inextpiv, 677 & TIPIV 678 & ) 679!$ USE OMP_LIB 680 USE MUMPS_OOC_COMMON 681 IMPLICIT NONE 682 INTEGER, intent(in) :: IBEG_BLOCK, IEND_BLOCK 683 INTEGER, intent(inout), OPTIONAL :: TIPIV(:) 684 INTEGER(8), intent(in) :: LA 685 COMPLEX(kind=8), intent(inout) :: A(LA) 686 INTEGER, intent(in) :: NFRONT,NASS,N,LIW,INODE,LAST_ROW 687 INTEGER, intent(inout) :: IFLAG,INOPV,NOFFW 688 DOUBLE PRECISION, intent(in) :: UU, SEUIL 689 INTEGER, intent(inout) :: IW(LIW) 690 INTEGER, intent(in) :: IOLDPS 691 INTEGER(8), intent(in) :: POSELT 692 INTEGER KEEP(500) 693 INTEGER(8) KEEP8(150) 694 INTEGER, intent(in) :: LPN_LIST 695 INTEGER, intent(inout) :: PIVNUL_LIST(LPN_LIST) 696 DOUBLE PRECISION DKEEP(230) 697 INTEGER PP_FIRST2SWAP_L, PP_LastPanelonDisk_L, 698 & PP_LastPIVRPTRFilled_L, 699 & PP_FIRST2SWAP_U, PP_LastPanelonDisk_U, 700 & PP_LastPIVRPTRFilled_U 701 INTEGER, intent(in) :: PIVOT_OPTION, IEND_BLR 702 LOGICAL, intent(in) :: LR_ACTIVATED 703 INTEGER, intent(inout) :: Inextpiv 704 INCLUDE 'mumps_headers.h' 705 COMPLEX(kind=8) SWOP 706 INTEGER XSIZE 707 INTEGER(8) :: APOS, IDIAG 708 INTEGER(8) :: J1, J2, JJ, J3 709 INTEGER(8) :: NFRONT8 710 INTEGER ILOC 711 COMPLEX(kind=8) ZERO 712 PARAMETER( ZERO = (0.0D0,0.0D0) ) 713 DOUBLE PRECISION RZERO, RMAX, AMROW 714 DOUBLE PRECISION PIVNUL 715 COMPLEX(kind=8) FIXA, CSEUIL 716 INTEGER NPIV,IPIV 717 INTEGER NPIVP1,JMAX,J,ISW,ISWPS1 718 INTEGER ISWPS2,KSW, HF 719 INTEGER ZMUMPS_IXAMAX 720 INTEGER :: ISHIFT, K206 721 INTEGER :: IPIV_SHIFT,IPIV_END 722 INTRINSIC max 723 DATA RZERO /0.0D0/ 724!$ INTEGER :: J4,JJMAX,NOMP,CHUNK,K361 725!$ DOUBLE PRECISION :: RRMAX,VALABS 726 INTEGER I_PIVRPTR_L, I_PIVR_L, NBPANELS_L 727 INTEGER I_PIVRPTR_U, I_PIVR_U, NBPANELS_U 728!$ NOMP = OMP_GET_MAX_THREADS() 729!$ K361 = KEEP(361) 730 PIVNUL = DKEEP(1) 731 FIXA = cmplx(DKEEP(2),kind=kind(FIXA)) 732 CSEUIL = cmplx(SEUIL,kind=kind(CSEUIL)) 733 NFRONT8 = int(NFRONT,8) 734 K206 = KEEP(206) 735 XSIZE = KEEP(IXSZ) 736 NPIV = IW(IOLDPS+1+XSIZE) 737 HF = 6 + IW(IOLDPS+5+XSIZE)+XSIZE 738 NPIVP1 = NPIV + 1 739 IF (KEEP(201).EQ.1) THEN 740 CALL ZMUMPS_GET_OOC_PERM_PTR(TYPEF_L, NBPANELS_L, 741 & I_PIVRPTR_L, I_PIVR_L, 742 & IOLDPS+2*NFRONT+6+IW(IOLDPS+5+XSIZE)+XSIZE, 743 & IW, LIW) 744 CALL ZMUMPS_GET_OOC_PERM_PTR(TYPEF_U, NBPANELS_U, 745 & I_PIVRPTR_U, I_PIVR_U, 746 & IOLDPS+2*NFRONT+6+IW(IOLDPS+5+XSIZE)+XSIZE, 747 & IW, LIW) 748 ENDIF 749 IF ( present(TIPIV) ) THEN 750 ILOC = NPIVP1 - IBEG_BLOCK + 1 751 TIPIV(ILOC) = ILOC 752 ENDIF 753 IF (INOPV .EQ. -1) THEN 754 APOS = POSELT + NFRONT8*int(NPIVP1-1,8) + int(NPIV,8) 755 IDIAG = APOS 756 IF(abs(A(APOS)).LT.SEUIL) THEN 757 IF (dble(A(APOS)) .GE. RZERO) THEN 758 A(APOS) = CSEUIL 759 ELSE 760 A(APOS) = -CSEUIL 761 ENDIF 762 KEEP(98) = KEEP(98)+1 763 ELSE IF (KEEP(258) .NE. 0) THEN 764 CALL ZMUMPS_UPDATEDETER(A(APOS), DKEEP(6), KEEP(259)) 765 ENDIF 766 IF (KEEP(201).EQ.1) THEN 767 IF (KEEP(251).EQ.0) THEN 768 CALL ZMUMPS_STORE_PERMINFO( IW(I_PIVRPTR_L), 769 & NBPANELS_L, 770 & IW(I_PIVR_L), NASS, NPIVP1, NPIVP1, 771 & PP_LastPanelonDisk_L, 772 & PP_LastPIVRPTRFilled_L) 773 ENDIF 774 CALL ZMUMPS_STORE_PERMINFO( IW(I_PIVRPTR_U), 775 & NBPANELS_U, 776 & IW(I_PIVR_U), NASS, NPIVP1, NPIVP1, 777 & PP_LastPanelonDisk_U, 778 & PP_LastPIVRPTRFilled_U) 779 ENDIF 780 GO TO 420 781 ENDIF 782 INOPV = 0 783 ISHIFT = 0 784 IPIV_END = IEND_BLOCK 785 IF (K206.GE.1) THEN 786 IF (Inextpiv.GT.NPIVP1.AND.Inextpiv.LE.IEND_BLOCK) THEN 787 ISHIFT = Inextpiv - NPIVP1 788 ENDIF 789 IF ( K206.EQ.1 790 & .OR. (K206 .GT.1 .AND. IEND_BLOCK.EQ.IEND_BLR) ) THEN 791 IPIV_END = IEND_BLOCK + ISHIFT 792 ENDIF 793 ENDIF 794 DO 460 IPIV_SHIFT = NPIVP1+ISHIFT, IPIV_END 795 IF (IPIV_SHIFT .LE. IEND_BLOCK) THEN 796 IPIV=IPIV_SHIFT 797 ELSE 798 IPIV = IPIV_SHIFT-IEND_BLOCK-1+NPIVP1 799 IF (IBEG_BLOCK.EQ.NPIVP1) THEN 800 EXIT 801 ENDIF 802 ENDIF 803 APOS = POSELT + NFRONT8*int(IPIV-1,8) + int(NPIV,8) 804 JMAX = 1 805 IF (PIVOT_OPTION.GT.0.AND.UU.GT.RZERO) GO TO 340 806 IF (A(APOS).EQ.ZERO) GO TO 630 807 GO TO 380 808 340 CONTINUE 809 AMROW = RZERO 810 J1 = APOS 811 IF (PIVOT_OPTION.EQ.1 .OR. 812 & (LR_ACTIVATED.AND.KEEP(480).GE.2)) THEN 813 J = IEND_BLR - NPIV 814 ELSE 815 J = NASS - NPIV 816 ENDIF 817 J2 = J1 + J - 1_8 818 IF (KEEP(351).EQ.1) THEN 819!$ IF (NOMP.GT.1 .AND. J.GE.K361) THEN 820!$ JMAX = 1 821!$ RMAX = RZERO 822!$ CHUNK = max(K361/2,J/NOMP) 823!$OMP PARALLEL PRIVATE(J3,VALABS,JJMAX,RRMAX) 824!$OMP& FIRSTPRIVATE(J1,J) 825!$ RRMAX = RZERO 826!$OMP DO schedule(static, CHUNK) 827!$ DO J4 = 1, J 828!$ J3 = J1 + int(J4-1,8) 829!$ VALABS = abs(A(J3)) 830!$ IF(VALABS.GT.RRMAX) THEN 831!$ RRMAX = VALABS 832!$ JJMAX = J4 833!$ ENDIF 834!$ ENDDO 835!$OMP END DO 836!$ IF (RRMAX.GT.0.0) THEN 837!$OMP CRITICAL 838!$ IF (RRMAX.GT.RMAX) THEN 839!$ RMAX = RRMAX 840!$ JMAX = JJMAX 841!$ ENDIF 842!$OMP END CRITICAL 843!$ ENDIF 844!$OMP END PARALLEL 845!$ ELSE 846 JMAX = ZMUMPS_IXAMAX(J,A(J1),1) 847!$ ENDIF 848 ELSE 849 JMAX = ZMUMPS_IXAMAX(J,A(J1),1) 850 ENDIF 851 JJ = J1 + int(JMAX - 1,8) 852 AMROW = abs(A(JJ)) 853 RMAX = AMROW 854 IF (PIVOT_OPTION.GE.2) THEN 855 J1 = J2 + 1_8 856 IF (PIVOT_OPTION.GE.3) THEN 857 J2 = APOS +int(- NPIV + NFRONT - 1 - KEEP(253),8) 858 ELSE 859 J2 = APOS +int(- NPIV + NASS - 1 - KEEP(253),8) 860 ENDIF 861 IF (J2.LT.J1) GO TO 370 862 IF (KEEP(351).EQ.1) THEN 863!$ CHUNK = max(K361/2,int(J2-J1)/NOMP) 864!$OMP PARALLEL DO schedule(static, CHUNK) PRIVATE(JJ) 865!$OMP& FIRSTPRIVATE(J1,J2) 866!$OMP& REDUCTION(max:RMAX) IF ((J2-J1).GE.K361) 867 DO JJ=J1,J2 868 RMAX = max(abs(A(JJ)),RMAX) 869 ENDDO 870!$OMP END PARALLEL DO 871 ELSE 872 DO 360 JJ=J1,J2 873 RMAX = max(abs(A(JJ)),RMAX) 874 360 CONTINUE 875 ENDIF 876 ENDIF 877 370 IDIAG = APOS + int(IPIV - NPIVP1,8) 878 IF ( RMAX .LE. PIVNUL ) THEN 879 IF (NFRONT - KEEP(253) .EQ. NASS) THEN 880 IF (IEND_BLOCK.NE.NASS ) THEN 881 GOTO 460 882 ENDIF 883 J1=POSELT+int(IPIV-1,8)+int(NPIV,8)*NFRONT8 884 J2=POSELT+int(IPIV-1,8)+int(LAST_ROW-1,8)*NFRONT8 885 ELSE 886 J1=POSELT+int(IPIV-1,8) 887 J2=POSELT+int(IPIV-1,8)+int(LAST_ROW-1,8)*NFRONT8 888 ENDIF 889 DO JJ=J1, J2, NFRONT8 890 IF ( abs(A(JJ)) .GT. PIVNUL ) THEN 891 GOTO 460 892 END IF 893 ENDDO 894 KEEP(109) = KEEP(109)+1 895 ISW = IOLDPS+HF+ 896 & IW(IOLDPS+1+XSIZE)+IPIV-NPIVP1 897 PIVNUL_LIST(KEEP(109)) = IW(ISW) 898 IF(dble(FIXA).GT.RZERO) THEN 899 IF(dble(A(IDIAG)) .GE. RZERO) THEN 900 A(IDIAG) = FIXA 901 ELSE 902 A(IDIAG) = -FIXA 903 ENDIF 904 ELSE 905 J1 = APOS 906 J2 = APOS +int(- NPIV + NFRONT - 1 - KEEP(253),8) 907 DO JJ=J1,J2 908 A(JJ) = ZERO 909 ENDDO 910 A(IDIAG) = -FIXA 911 ENDIF 912 JMAX = IPIV - NPIV 913 GOTO 385 914 ENDIF 915 IF (abs(A(IDIAG)) .GE. UU*RMAX .AND. 916 & abs(A(IDIAG)) .GT. max(SEUIL,tiny(RMAX))) THEN 917 JMAX = IPIV - NPIV 918 GO TO 380 919 ENDIF 920 IF ( .NOT. (AMROW .GE. UU*RMAX .AND. 921 & AMROW .GT. max(SEUIL,tiny(RMAX))) ) GO TO 460 922 NOFFW = NOFFW + 1 923 380 CONTINUE 924 IF (K206.GE.1) THEN 925 Inextpiv = IPIV + 1 926 ENDIF 927 IF (KEEP(258) .NE. 0) THEN 928 CALL ZMUMPS_UPDATEDETER( A(APOS+int(JMAX-1,8)), 929 & DKEEP(6), 930 & KEEP(259)) 931 ENDIF 932 385 CONTINUE 933 IF (IPIV.EQ.NPIVP1) GO TO 400 934 KEEP(260)=-KEEP(260) 935 J1 = POSELT + int(NPIV,8)*NFRONT8 936 J2 = J1 + NFRONT8 - 1_8 937 J3 = POSELT + int(IPIV-1,8)*NFRONT8 938 DO 390 JJ=J1,J2 939 SWOP = A(JJ) 940 A(JJ) = A(J3) 941 A(J3) = SWOP 942 J3 = J3 + 1_8 943 390 CONTINUE 944 ISWPS1 = IOLDPS + HF - 1 + NPIVP1 945 ISWPS2 = IOLDPS + HF - 1 + IPIV 946 ISW = IW(ISWPS1) 947 IW(ISWPS1) = IW(ISWPS2) 948 IW(ISWPS2) = ISW 949 400 IF (JMAX.EQ.1) GO TO 420 950 KEEP(260)=-KEEP(260) 951 IF ( present(TIPIV) ) THEN 952 TIPIV(ILOC) = ILOC + JMAX - 1 953 ENDIF 954 J1 = POSELT + int(NPIV,8) 955 J2 = POSELT + int(NPIV + JMAX - 1,8) 956 DO 410 KSW=1,LAST_ROW 957 SWOP = A(J1) 958 A(J1) = A(J2) 959 A(J2) = SWOP 960 J1 = J1 + NFRONT8 961 J2 = J2 + NFRONT8 962 410 CONTINUE 963 ISWPS1 = IOLDPS + HF - 1 + NFRONT + NPIV + 1 964 ISWPS2 = IOLDPS + HF - 1 + NFRONT + NPIV + JMAX 965 ISW = IW(ISWPS1) 966 IW(ISWPS1) = IW(ISWPS2) 967 IW(ISWPS2) = ISW 968 GO TO 420 969 460 CONTINUE 970 IF (K206 .GE. 1) THEN 971 Inextpiv=IEND_BLOCK+1 972 ENDIF 973 IF (IEND_BLOCK.EQ.NASS) THEN 974 INOPV = 1 975 ELSE 976 INOPV = 2 977 ENDIF 978 GO TO 430 979 630 CONTINUE 980 IFLAG = -10 981 GOTO 430 982 420 CONTINUE 983 IF (KEEP(201).EQ.1) THEN 984 IF (KEEP(251).EQ.0) THEN 985 CALL ZMUMPS_STORE_PERMINFO( IW(I_PIVRPTR_L), 986 & NBPANELS_L, 987 & IW(I_PIVR_L), NASS, NPIVP1, IPIV, 988 & PP_LastPanelonDisk_L, 989 & PP_LastPIVRPTRFilled_L) 990 ENDIF 991 CALL ZMUMPS_STORE_PERMINFO( IW(I_PIVRPTR_U), 992 & NBPANELS_U, 993 & IW(I_PIVR_U), NASS, NPIVP1, NPIV+JMAX, 994 & PP_LastPanelonDisk_U, 995 & PP_LastPIVRPTRFilled_U) 996 ENDIF 997 430 CONTINUE 998 RETURN 999 END SUBROUTINE ZMUMPS_FAC_I 1000 SUBROUTINE ZMUMPS_FAC_I_LDLT 1001 & (NFRONT,NASS,INODE,IBEG_BLOCK,IEND_BLOCK, 1002 & IW,LIW, A,LA, INOPV, 1003 & NNEG, 1004 & IFLAG,IOLDPS,POSELT,UU, SEUIL,KEEP,KEEP8,PIVSIZ, 1005 & DKEEP,PIVNUL_LIST,LPN_LIST, XSIZE, 1006 & PP_FIRST2SWAP_L, PP_LastPanelonDisk, 1007 & PP_LastPIVRPTRIndexFilled,MAXFROMM,IS_MAXFROMM_AVAIL, 1008 & PIVOT_OPTION, IEND_BLR, Inextpiv) 1009 USE MUMPS_OOC_COMMON 1010 IMPLICIT NONE 1011 INTEGER(8) :: POSELT, LA 1012 INTEGER NFRONT,NASS,LIW,INODE,IFLAG,INOPV, 1013 & IOLDPS, NNEG 1014 INTEGER, intent(in) :: IBEG_BLOCK, IEND_BLOCK 1015 INTEGER, intent(in) :: PIVOT_OPTION,IEND_BLR 1016 INTEGER, intent(inout) :: Inextpiv 1017 INTEGER PIVSIZ,LPIV, XSIZE 1018 COMPLEX(kind=8) A(LA) 1019 DOUBLE PRECISION UU, UULOC, SEUIL 1020 INTEGER IW(LIW) 1021 INTEGER KEEP(500) 1022 INTEGER(8) KEEP8(150) 1023 INTEGER LPN_LIST 1024 INTEGER PIVNUL_LIST(LPN_LIST) 1025 DOUBLE PRECISION DKEEP(230) 1026 INTEGER PP_FIRST2SWAP_L, PP_LastPanelonDisk 1027 INTEGER PP_LastPIVRPTRIndexFilled 1028 DOUBLE PRECISION, intent(in) :: MAXFROMM 1029 LOGICAL, intent(inout) :: IS_MAXFROMM_AVAIL 1030 include 'mpif.h' 1031 INTEGER (8) :: POSPV1,POSPV2,OFFDAG,APOSJ 1032 INTEGER JMAX, LIM 1033 DOUBLE PRECISION RMAX,AMAX,TMAX 1034 DOUBLE PRECISION MAXPIV 1035 DOUBLE PRECISION PIVNUL 1036 COMPLEX(kind=8) FIXA, CSEUIL 1037 COMPLEX(kind=8) PIVOT,DETPIV 1038 INCLUDE 'mumps_headers.h' 1039 INTEGER :: J 1040 INTEGER(8) :: APOS, J1, J2, JJ, NFRONT8, KK, J1_ini, JJ_ini 1041 INTEGER :: LDA 1042 INTEGER(8) :: LDA8 1043 INTEGER NPIV,IPIV 1044 INTEGER NPIVP1,K 1045 INTEGER :: ISHIFT, K206, IPIV_SHIFT, IPIV_END 1046 INTRINSIC max 1047 COMPLEX(kind=8) ZERO, ONE 1048 PARAMETER( ZERO = (0.0D0,0.0D0) ) 1049 PARAMETER( ONE = (1.0D0,1.0D0) ) 1050 DOUBLE PRECISION RZERO,RONE 1051 PARAMETER(RZERO=0.0D0, RONE=1.0D0) 1052 LOGICAL OMP_FLAG 1053 INTEGER I_PIVRPTR, I_PIVR, NBPANELS_L 1054 PIVNUL = DKEEP(1) 1055 FIXA = cmplx(DKEEP(2),kind=kind(FIXA)) 1056 CSEUIL = cmplx(SEUIL,kind=kind(CSEUIL)) 1057 LDA = NFRONT 1058 LDA8 = int(LDA,8) 1059 NFRONT8 = int(NFRONT,8) 1060 K206 = KEEP(206) 1061 IF (KEEP(201).EQ.1 .AND. KEEP(50).NE.1) THEN 1062 CALL ZMUMPS_GET_OOC_PERM_PTR(TYPEF_L, NBPANELS_L, 1063 & I_PIVRPTR, I_PIVR, IOLDPS+2*NFRONT+6+KEEP(IXSZ), 1064 & IW, LIW) 1065 ENDIF 1066 UULOC = UU 1067 PIVSIZ = 1 1068 NPIV = IW(IOLDPS+1+XSIZE) 1069 NPIVP1 = NPIV + 1 1070 IF(INOPV .EQ. -1) THEN 1071 APOS = POSELT + (LDA8+1_8) * int(NPIV,8) 1072 POSPV1 = APOS 1073 IF(abs(A(APOS)).LT.SEUIL) THEN 1074 IF(dble(A(APOS)) .GE. RZERO) THEN 1075 A(APOS) = CSEUIL 1076 ELSE 1077 A(APOS) = -CSEUIL 1078 ENDIF 1079 KEEP(98) = KEEP(98)+1 1080 ELSE IF (KEEP(258) .NE. 0) THEN 1081 CALL ZMUMPS_UPDATEDETER( A(APOS), DKEEP(6), KEEP(259) ) 1082 ENDIF 1083 IF (KEEP(201).EQ.1.AND.KEEP(50).NE.1) THEN 1084 CALL ZMUMPS_STORE_PERMINFO( IW(I_PIVRPTR), NBPANELS_L, 1085 & IW(I_PIVR), NASS, NPIVP1, NPIVP1, 1086 & PP_LastPanelonDisk, 1087 & PP_LastPIVRPTRIndexFilled) 1088 ENDIF 1089 GO TO 420 1090 ENDIF 1091 INOPV = 0 1092 ISHIFT = 0 1093 IPIV_END = IEND_BLOCK 1094 IF (K206.GE.1) THEN 1095 IF (Inextpiv.GT.NPIVP1.AND.Inextpiv.LE.IEND_BLOCK) THEN 1096 ISHIFT = Inextpiv - NPIVP1 1097 ENDIF 1098 IF ( K206.EQ.1 1099 & .OR. (K206 .GT.1 .AND. IEND_BLOCK.EQ.IEND_BLR) ) THEN 1100 IPIV_END = IEND_BLOCK + ISHIFT 1101 ENDIF 1102 IF (ISHIFT.GT.0.AND.IS_MAXFROMM_AVAIL) THEN 1103 IPIV = NPIVP1 1104 APOS = POSELT + LDA8*int(IPIV-1,8) + int(NPIV,8) 1105 POSPV1 = APOS + int(IPIV - NPIVP1,8) 1106 PIVOT = A(POSPV1) 1107 IF ( MAXFROMM .GT. PIVNUL ) THEN 1108 IF ( abs(PIVOT) .GE. UULOC*MAXFROMM 1109 & .AND. abs(PIVOT) .GT. max(SEUIL,tiny(MAXFROMM)) ) THEN 1110 ISHIFT = 0 1111 ENDIF 1112 ENDIF 1113 ENDIF 1114 IF ( ISHIFT .GT. 0) THEN 1115 IS_MAXFROMM_AVAIL = .FALSE. 1116 ENDIF 1117 ENDIF 1118 DO 460 IPIV_SHIFT = NPIVP1+ISHIFT, IPIV_END 1119 IF (IPIV_SHIFT .LE. IEND_BLOCK) THEN 1120 IPIV=IPIV_SHIFT 1121 ELSE 1122 IPIV = IPIV_SHIFT-IEND_BLOCK-1+NPIVP1 1123 IF (IBEG_BLOCK.EQ.NPIVP1) THEN 1124 EXIT 1125 ENDIF 1126 ENDIF 1127 APOS = POSELT + LDA8*int(IPIV-1,8) + int(NPIV,8) 1128 POSPV1 = APOS + int(IPIV - NPIVP1,8) 1129 PIVOT = A(POSPV1) 1130 IF (UULOC.EQ.RZERO.OR.PIVOT_OPTION.EQ.0) THEN 1131 IF (abs(A(APOS)).EQ.RZERO) GO TO 630 1132 IF (KEEP(258) .NE. 0) THEN 1133 CALL ZMUMPS_UPDATEDETER(A(APOS), DKEEP(6), KEEP(259)) 1134 ENDIF 1135 GO TO 420 1136 ENDIF 1137 IF ( IS_MAXFROMM_AVAIL ) THEN 1138 IF ( MAXFROMM .GT. PIVNUL ) THEN 1139 IF ( abs(PIVOT) .GE. UULOC*MAXFROMM 1140 & .AND. abs(PIVOT) .GT. max(SEUIL,tiny(MAXFROMM)) ) THEN 1141 IF (KEEP(258) .NE. 0) THEN 1142 CALL ZMUMPS_UPDATEDETER(PIVOT, DKEEP(6), KEEP(259)) 1143 ENDIF 1144 GOTO 415 1145 ENDIF 1146 ENDIF 1147 IS_MAXFROMM_AVAIL = .FALSE. 1148 ENDIF 1149 AMAX = -RONE 1150 JMAX = 0 1151 J1 = APOS 1152 J2 = POSPV1 - 1_8 1153 DO JJ=J1,J2 1154 IF(abs(A(JJ)) .GT. AMAX) THEN 1155 AMAX = abs(A(JJ)) 1156 JMAX = IPIV - int(POSPV1-JJ) 1157 ENDIF 1158 ENDDO 1159 J1 = POSPV1 + LDA8 1160 DO J=1, IEND_BLOCK - IPIV 1161 IF(abs(A(J1)) .GT. AMAX) THEN 1162 AMAX = abs(A(J1)) 1163 JMAX = IPIV + J 1164 ENDIF 1165 J1 = J1 + LDA8 1166 ENDDO 1167 RMAX = RZERO 1168 IF (PIVOT_OPTION.EQ.3) THEN 1169 LIM = NFRONT 1170 ELSEIF (PIVOT_OPTION.EQ.2) THEN 1171 LIM = NASS 1172 ELSEIF (PIVOT_OPTION.EQ.1) THEN 1173 LIM = IEND_BLR 1174 ELSE 1175 write(*,*) 'Internal error in FAC_I_LDLT: PIVOT_OPTION=', 1176 & PIVOT_OPTION 1177 ENDIF 1178 J1_ini = J1 1179 IF ( (LIM - KEEP(253) - IEND_BLOCK).GE.300 ) THEN 1180 OMP_FLAG = .TRUE. 1181 ELSE 1182 OMP_FLAG = .FALSE. 1183 ENDIF 1184!$OMP PARALLEL DO PRIVATE(J1) REDUCTION(max:RMAX) IF(OMP_FLAG) 1185 DO J=1, LIM - KEEP(253) - IEND_BLOCK 1186 J1 = J1_ini + int(J-1,8) * LDA8 1187 RMAX = max(abs(A(J1)),RMAX) 1188 ENDDO 1189!$OMP END PARALLEL DO 1190 IF (max(AMAX,RMAX,abs(PIVOT)).LE.PIVNUL) THEN 1191 KEEP(109) = KEEP(109)+1 1192 PIVNUL_LIST(KEEP(109)) = -1 1193 IF(dble(FIXA).GT.RZERO) THEN 1194 IF(dble(PIVOT) .GE. RZERO) THEN 1195 A(POSPV1) = FIXA 1196 ELSE 1197 A(POSPV1) = -FIXA 1198 ENDIF 1199 ELSE 1200 J1 = APOS 1201 J2 = POSPV1 - 1_8 1202 DO JJ=J1,J2 1203 A(JJ) = ZERO 1204 ENDDO 1205 J1 = POSPV1 + LDA8 1206 DO J=1, IEND_BLOCK - IPIV 1207 A(J1) = ZERO 1208 J1 = J1 + LDA8 1209 ENDDO 1210 DO J=1,NFRONT - IEND_BLOCK 1211 A(J1) = ZERO 1212 J1 = J1 + LDA8 1213 ENDDO 1214 A(POSPV1) = ONE 1215 ENDIF 1216 PIVOT = A(POSPV1) 1217 GO TO 415 1218 ENDIF 1219 IF ( abs(PIVOT).GE.UULOC*max(RMAX,AMAX) 1220 & .AND. abs(PIVOT) .GT. max(SEUIL,tiny(RMAX)) ) THEN 1221 IF (KEEP(258) .NE.0 ) THEN 1222 CALL ZMUMPS_UPDATEDETER(PIVOT, DKEEP(6), KEEP(259)) 1223 ENDIF 1224 GO TO 415 1225 END IF 1226 IF (NPIVP1.EQ.IEND_BLOCK) THEN 1227 GOTO 460 1228 ENDIF 1229 IF (max(abs(PIVOT),RMAX,AMAX).LE.tiny(RMAX)) THEN 1230 GOTO 460 1231 ENDIF 1232 IF ( 1233 & (KEEP(19).NE.0).AND.(max(AMAX,RMAX,abs(PIVOT)).LE.SEUIL) 1234 & ) 1235 & THEN 1236 GO TO 460 1237 ENDIF 1238 IF (RMAX.LT.AMAX) THEN 1239 J1 = APOS 1240 J2 = POSPV1 - 1_8 1241 DO JJ=J1,J2 1242 IF(int(POSPV1-JJ) .NE. IPIV-JMAX) THEN 1243 RMAX = max(RMAX,abs(A(JJ))) 1244 ENDIF 1245 ENDDO 1246 J1 = POSPV1 + LDA8 1247 DO J=1,NASS-IPIV 1248 IF(IPIV+J .NE. JMAX) THEN 1249 RMAX = max(abs(A(J1)),RMAX) 1250 ENDIF 1251 J1 = J1 + LDA8 1252 ENDDO 1253 ENDIF 1254 IF (PIVOT_OPTION.EQ.3) THEN 1255 LIM = NFRONT 1256 ELSEIF (PIVOT_OPTION.EQ.2) THEN 1257 LIM = NASS 1258 ELSEIF (PIVOT_OPTION.EQ.1) THEN 1259 LIM = IEND_BLR 1260 ELSE 1261 write(*,*) 'Internal error in FAC_I_LDLT: PIVOT_OPTION=', 1262 & PIVOT_OPTION 1263 ENDIF 1264 APOSJ = POSELT + int(JMAX-1,8)*LDA8 + int(NPIV,8) 1265 POSPV2 = APOSJ + int(JMAX - NPIVP1,8) 1266 IF (IPIV.LT.JMAX) THEN 1267 OFFDAG = APOSJ + int(IPIV - NPIVP1,8) 1268 ELSE 1269 OFFDAG = APOS + int(JMAX - NPIVP1,8) 1270 END IF 1271 TMAX = RZERO 1272 IF (JMAX .LT. IPIV) THEN 1273 JJ_ini = POSPV2 1274 OMP_FLAG = (LIM-JMAX-KEEP(253). GE. 300) 1275!$OMP PARALLEL DO IF (OMP_FLAG) 1276!$OMP& PRIVATE(JJ) REDUCTION(max:TMAX) 1277 DO K = 1, LIM - JMAX - KEEP(253) 1278 JJ = JJ_ini+ int(K,8)*NFRONT8 1279 IF (JMAX+K.NE.IPIV) THEN 1280 TMAX=max(TMAX,abs(A(JJ))) 1281 ENDIF 1282 ENDDO 1283!$OMP END PARALLEL DO 1284 DO KK = APOSJ, POSPV2-1_8 1285 TMAX = max(TMAX,abs(A(KK))) 1286 ENDDO 1287 ELSE 1288 JJ_ini = POSPV2 1289 OMP_FLAG = (LIM-JMAX-KEEP(253). GE. 300) 1290!$OMP PARALLEL DO PRIVATE(JJ) REDUCTION(max:TMAX) IF(OMP_FLAG) 1291 DO K = 1, LIM-JMAX-KEEP(253) 1292 JJ = JJ_ini + int(K,8)*NFRONT8 1293 TMAX=max(TMAX,abs(A(JJ))) 1294 ENDDO 1295!$OMP END PARALLEL DO 1296 DO KK = APOSJ, POSPV2 - 1_8 1297 IF (KK.NE.OFFDAG) THEN 1298 TMAX = max(TMAX,abs(A(KK))) 1299 ENDIF 1300 ENDDO 1301 ENDIF 1302 DETPIV = A(POSPV1)*A(POSPV2) - A(OFFDAG)**2 1303 IF (SEUIL.GT.RZERO) THEN 1304 IF (sqrt(abs(DETPIV)) .LE. SEUIL ) THEN 1305 GOTO 460 1306 ENDIF 1307 ENDIF 1308 MAXPIV = max(abs(A(POSPV1)),abs(A(POSPV2))) 1309 IF (MAXPIV.EQ.RZERO) MAXPIV = RONE 1310 IF ((abs(A(POSPV2))*RMAX+AMAX*TMAX)*UULOC.GT. 1311 & abs(DETPIV) .OR. abs(DETPIV) .EQ. RZERO) THEN 1312 GO TO 460 1313 ENDIF 1314 IF ((abs(A(POSPV1))*TMAX+AMAX*RMAX)*UULOC.GT. 1315 & abs(DETPIV) .OR. abs(DETPIV) .EQ. RZERO) THEN 1316 GO TO 460 1317 ENDIF 1318 IF (KEEP(258) .NE.0 ) THEN 1319 CALL ZMUMPS_UPDATEDETER(DETPIV, DKEEP(6), KEEP(259)) 1320 ENDIF 1321 PIVSIZ = 2 1322 KEEP(103) = KEEP(103)+1 1323 415 CONTINUE 1324 IF (K206.GE.1) THEN 1325 Inextpiv = max(NPIVP1+PIVSIZ, IPIV+1) 1326 ENDIF 1327 DO K=1,PIVSIZ 1328 IF (PIVSIZ .EQ. 2) THEN 1329 IF (K==1) THEN 1330 LPIV = min(IPIV,JMAX) 1331 ELSE 1332 LPIV = max(IPIV,JMAX) 1333 ENDIF 1334 ELSE 1335 LPIV = IPIV 1336 ENDIF 1337 IF (LPIV.EQ.NPIVP1) GOTO 416 1338 CALL ZMUMPS_SWAP_LDLT( A, LA, IW, LIW, 1339 & IOLDPS, NPIVP1, LPIV, POSELT, NASS, 1340 & LDA, NFRONT, 1, KEEP(219), KEEP(50), 1341 & KEEP(IXSZ), -9999) 1342 416 CONTINUE 1343 IF (KEEP(201).EQ.1.AND.KEEP(50).NE.1) THEN 1344 CALL ZMUMPS_STORE_PERMINFO( IW(I_PIVRPTR), NBPANELS_L, 1345 & IW(I_PIVR), NASS, NPIVP1, LPIV, PP_LastPanelonDisk, 1346 & PP_LastPIVRPTRIndexFilled) 1347 ENDIF 1348 NPIVP1 = NPIVP1 + 1 1349 ENDDO 1350 IF(PIVSIZ .EQ. 2) THEN 1351 A(POSELT+(LDA8+1_8)*int(NPIV,8)+1_8) = DETPIV 1352 ENDIF 1353 GOTO 420 1354 460 CONTINUE 1355 IF (K206 .GE. 1) THEN 1356 Inextpiv=IEND_BLOCK+1 1357 ENDIF 1358 IF (IEND_BLOCK.EQ.NASS) THEN 1359 INOPV = 1 1360 ELSE 1361 INOPV = 2 1362 ENDIF 1363 GO TO 420 1364 630 CONTINUE 1365 PIVSIZ = 0 1366 IFLAG = -10 1367 420 CONTINUE 1368 IS_MAXFROMM_AVAIL = .FALSE. 1369 RETURN 1370 END SUBROUTINE ZMUMPS_FAC_I_LDLT 1371 SUBROUTINE ZMUMPS_FAC_MQ_LDLT(IEND_BLOCK, 1372 & NFRONT,NASS,NPIV,INODE, 1373 & A,LA,LDA, 1374 & POSELT,IFINB,PIVSIZ, 1375 & MAXFROMM, IS_MAXFROMM_AVAIL, IS_MAX_USEFUL, 1376 & KEEP253, PIVOT_OPTION, IEND_BLR 1377 & ) 1378 IMPLICIT NONE 1379 INTEGER, intent(out):: IFINB 1380 INTEGER, intent(in) :: INODE, NFRONT, NASS, NPIV 1381 INTEGER, intent(in) :: IEND_BLOCK 1382 INTEGER, intent(in) :: LDA 1383 INTEGER(8), intent(in) :: LA 1384 COMPLEX(kind=8), intent(inout) :: A(LA) 1385 INTEGER, intent(in) :: PIVOT_OPTION, IEND_BLR 1386 INTEGER(8) :: POSELT 1387 DOUBLE PRECISION, intent(out) :: MAXFROMM 1388 LOGICAL, intent(out) :: IS_MAXFROMM_AVAIL 1389 LOGICAL, intent(in) :: IS_MAX_USEFUL 1390 INTEGER, INTENT(in) :: KEEP253 1391 COMPLEX(kind=8) VALPIV 1392 DOUBLE PRECISION :: MAXFROMMTMP 1393 INTEGER NCB1 1394 INTEGER(8) :: NFRONT8 1395 INTEGER(8) :: LDA8 1396 INTEGER(8) :: K1POS 1397 INTEGER NEL2, NEL, LIM 1398 COMPLEX(kind=8) ONE, ZERO 1399 COMPLEX(kind=8) A11,A22,A12 1400 INTEGER(8) :: APOS, LPOS, LPOS1, LPOS2 1401 INTEGER(8) :: POSPV1, POSPV2 1402 INTEGER PIVSIZ,NPIV_NEW,J2,I 1403 INTEGER(8) :: OFFDAG, OFFDAG_OLD, IBEG, IEND 1404 INTEGER(8) :: JJ, K1, K2, IROW 1405 COMPLEX(kind=8) SWOP,DETPIV,MULT1,MULT2 1406 INCLUDE 'mumps_headers.h' 1407 PARAMETER(ONE = (1.0D0,0.0D0), 1408 & ZERO = (0.0D0,0.0D0)) 1409 LDA8 = int(LDA,8) 1410 NFRONT8 = int(NFRONT,8) 1411 NPIV_NEW = NPIV + PIVSIZ 1412 NEL = NFRONT - NPIV_NEW 1413 IFINB = 0 1414 IS_MAXFROMM_AVAIL = .FALSE. 1415 NEL2 = IEND_BLOCK - NPIV_NEW 1416 IF (NEL2.EQ.0) THEN 1417 IF (IEND_BLOCK.EQ.NASS) THEN 1418 IFINB = -1 1419 ELSE 1420 IFINB = 1 1421 ENDIF 1422 ENDIF 1423 IF(PIVSIZ .EQ. 1) THEN 1424 APOS = POSELT + int(NPIV,8)*(NFRONT8 + 1_8) 1425 VALPIV = ONE/A(APOS) 1426 LPOS = APOS + LDA8 1427 MAXFROMM = 0.0D00 1428 IF (NEL2 > 0) THEN 1429 IF (.NOT. IS_MAX_USEFUL) THEN 1430 DO I=1, NEL2 1431 K1POS = LPOS + int(I-1,8)*LDA8 1432 A(APOS+int(I,8))=A(K1POS) 1433 A(K1POS) = A(K1POS) * VALPIV 1434 DO JJ=1_8, int(I,8) 1435 A(K1POS+JJ)=A(K1POS+JJ) - A(K1POS) * A(APOS+JJ) 1436 ENDDO 1437 ENDDO 1438 ELSE 1439 IS_MAXFROMM_AVAIL = .TRUE. 1440 DO I=1, NEL2 1441 K1POS = LPOS + int(I-1,8)*LDA8 1442 A(APOS+int(I,8))=A(K1POS) 1443 A(K1POS) = A(K1POS) * VALPIV 1444 A(K1POS+1_8)=A(K1POS+1_8) - A(K1POS) * A(APOS+1_8) 1445 MAXFROMM=max( MAXFROMM,abs(A(K1POS+1_8)) ) 1446 DO JJ = 2_8, int(I,8) 1447 A(K1POS+JJ)=A(K1POS+JJ) - A(K1POS) * A(APOS+JJ) 1448 ENDDO 1449 ENDDO 1450 ENDIF 1451 ENDIF 1452 IF (PIVOT_OPTION.EQ.3) THEN 1453 LIM = NFRONT 1454 ELSEIF (PIVOT_OPTION.EQ.2) THEN 1455 LIM = NASS 1456 ELSE 1457 LIM = IEND_BLR 1458 ENDIF 1459 NCB1 = LIM - IEND_BLOCK 1460 IF (.NOT. IS_MAX_USEFUL) THEN 1461!$OMP PARALLEL DO PRIVATE(JJ,K1POS) IF (NCB1 > 300) 1462 DO I=NEL2+1, NEL2 + NCB1 1463 K1POS = LPOS+ int(I-1,8)*LDA8 1464 A(APOS+int(I,8))=A(K1POS) 1465 A(K1POS) = A(K1POS) * VALPIV 1466 DO JJ = 1_8, int(NEL2,8) 1467 A(K1POS+JJ)=A(K1POS+JJ) - A(K1POS) * A(APOS+JJ) 1468 ENDDO 1469 ENDDO 1470!$OMP END PARALLEL DO 1471 ELSE 1472 MAXFROMMTMP=0.0D0 1473!$OMP PARALLEL DO PRIVATE(JJ,K1POS) 1474!$OMP& REDUCTION(max:MAXFROMMTMP) IF (NCB1 - KEEP253 > 300) 1475 DO I=NEL2+1, NEL2 + NCB1 - KEEP253 1476 K1POS = LPOS+ int(I-1,8)*LDA8 1477 A(APOS+int(I,8))=A(K1POS) 1478 A(K1POS) = A(K1POS) * VALPIV 1479 IF (NEL2 > 0) THEN 1480 A(K1POS+1_8) = A(K1POS+1_8) - A(K1POS) * A(APOS+1_8) 1481 MAXFROMMTMP=max(MAXFROMMTMP, abs(A(K1POS+1_8))) 1482 DO JJ = 2_8, int(NEL2,8) 1483 A(K1POS+JJ)=A(K1POS+JJ) - A(K1POS) * A(APOS+JJ) 1484 ENDDO 1485 ENDIF 1486 ENDDO 1487!$OMP END PARALLEL DO 1488 DO I = NEL2 + NCB1 - KEEP253 + 1, NEL2 + NCB1 1489 K1POS = LPOS+ int(I-1,8)*LDA8 1490 A(APOS+int(I,8))=A(K1POS) 1491 A(K1POS) = A(K1POS) * VALPIV 1492 DO JJ = 1_8, int(NEL2,8) 1493 A(K1POS+JJ)=A(K1POS+JJ) - A(K1POS) * A(APOS+JJ) 1494 ENDDO 1495 ENDDO 1496 MAXFROMM=max(MAXFROMM, MAXFROMMTMP) 1497 ENDIF 1498 ELSE 1499 IF (PIVOT_OPTION.EQ.3) THEN 1500 LIM = NFRONT 1501 ELSEIF (PIVOT_OPTION.EQ.2) THEN 1502 LIM = NASS 1503 ELSE 1504 LIM = IEND_BLR 1505 ENDIF 1506 POSPV1 = POSELT + int(NPIV,8)*(NFRONT8 + 1_8) 1507 POSPV2 = POSPV1 + NFRONT8 + 1_8 1508 OFFDAG_OLD = POSPV2 - 1_8 1509 OFFDAG = POSPV1 + 1_8 1510 SWOP = A(POSPV2) 1511 DETPIV = A(OFFDAG) 1512 A22 = A(POSPV1)/DETPIV 1513 A11 = SWOP/DETPIV 1514 A12 = -A(OFFDAG_OLD)/DETPIV 1515 A(OFFDAG) = A(OFFDAG_OLD) 1516 A(OFFDAG_OLD) = ZERO 1517 LPOS1 = POSPV2 + LDA8 - 1_8 1518 LPOS2 = LPOS1 + 1_8 1519 CALL zcopy(LIM-NPIV_NEW, A(LPOS1), LDA, A(POSPV1+2_8), 1) 1520 CALL zcopy(LIM-NPIV_NEW, A(LPOS2), LDA, A(POSPV2+1_8), 1) 1521 JJ = POSPV2 + NFRONT8-1_8 1522 IBEG = JJ + 2_8 1523 IEND = IBEG 1524 DO J2 = 1,NEL2 1525 K1 = JJ 1526 K2 = JJ+1_8 1527 MULT1 = - (A11*A(K1)+A12*A(K2)) 1528 MULT2 = - (A12*A(K1)+A22*A(K2)) 1529 K1 = POSPV1 + 2_8 1530 K2 = POSPV2 + 1_8 1531 DO IROW = IBEG, IEND 1532 A(IROW) = A(IROW) + MULT1*A(K1) + MULT2*A(K2) 1533 K1 = K1 + 1_8 1534 K2 = K2 + 1_8 1535 ENDDO 1536 A( JJ ) = -MULT1 1537 A( JJ + 1_8 ) = -MULT2 1538 IBEG = IBEG + NFRONT8 1539 IEND = IEND + NFRONT8 + 1_8 1540 JJ = JJ+NFRONT8 1541 ENDDO 1542 IEND = IEND-1_8 1543 DO J2 = IEND_BLOCK+1,LIM 1544 K1 = JJ 1545 K2 = JJ+1_8 1546 MULT1 = - (A11*A(K1)+A12*A(K2)) 1547 MULT2 = - (A12*A(K1)+A22*A(K2)) 1548 K1 = POSPV1 + 2_8 1549 K2 = POSPV2 + 1_8 1550 DO IROW = IBEG, IEND 1551 A(IROW) = A(IROW) + MULT1*A(K1) + MULT2*A(K2) 1552 K1 = K1 + 1_8 1553 K2 = K2 + 1_8 1554 ENDDO 1555 A( JJ ) = -MULT1 1556 A( JJ + 1_8 ) = -MULT2 1557 IBEG = IBEG + NFRONT8 1558 IEND = IEND + NFRONT8 1559 JJ = JJ + NFRONT8 1560 ENDDO 1561 ENDIF 1562 RETURN 1563 END SUBROUTINE ZMUMPS_FAC_MQ_LDLT 1564 SUBROUTINE ZMUMPS_FAC_SQ_LDLT(IBEG_BLOCK,IEND_BLOCK,NPIV, 1565 & NFRONT,NASS,LAST_VAR,INODE,A,LA, 1566 & LDA, 1567 & POSELT, 1568 & KEEP,KEEP8, 1569 & PIVOT_OPTION, CALL_TRSM) 1570 IMPLICIT NONE 1571 INTEGER, intent(in) :: NPIV 1572 INTEGER, intent(in) :: NFRONT, NASS, IBEG_BLOCK, IEND_BLOCK 1573 INTEGER(8), intent(in) :: LA 1574 COMPLEX(kind=8), intent(inout) :: A(LA) 1575 INTEGER, intent(in) :: INODE 1576 INTEGER, intent(in) :: LAST_VAR 1577 INTEGER :: KEEP(500) 1578 INTEGER(8) KEEP8(150) 1579 INTEGER(8), intent(in) :: POSELT 1580 INTEGER, intent(in) :: LDA 1581 INTEGER, intent(in) :: PIVOT_OPTION 1582 LOGICAL, intent(in) :: CALL_TRSM 1583 INTEGER(8) :: LDA8 1584 INTEGER NPIV_BLOCK, NEL1, I, II 1585 INTEGER(8) :: LPOS,UPOS,APOS 1586 INTEGER IROW 1587 INTEGER Block 1588 INTEGER BLSIZE, ELSIZE 1589 COMPLEX(kind=8) ONE, ALPHA, VALPIV 1590 INCLUDE 'mumps_headers.h' 1591 PARAMETER (ONE=(1.0D0,0.0D0), ALPHA=(-1.0D0,0.0D0)) 1592 LDA8 = int(LDA,8) 1593 ELSIZE = IEND_BLOCK - IBEG_BLOCK +1 1594 NEL1 = LAST_VAR - IEND_BLOCK 1595 NPIV_BLOCK = NPIV - IBEG_BLOCK + 1 1596 IF (NPIV_BLOCK.EQ.0) GO TO 500 1597 IF (NEL1.NE.0) THEN 1598 IF (PIVOT_OPTION.LE.1.AND.CALL_TRSM) THEN 1599 APOS = POSELT + LDA8*int(IBEG_BLOCK-1,8) + int(IBEG_BLOCK-1,8) 1600 LPOS = POSELT + LDA8*int(IEND_BLOCK,8) + int(IBEG_BLOCK-1,8) 1601 UPOS = POSELT + LDA8*int(IBEG_BLOCK-1,8) + int(IEND_BLOCK,8) 1602 CALL ztrsm('L', 'U', 'T', 'U', ELSIZE, NEL1, ONE, 1603 & A(APOS), LDA, A(LPOS), LDA) 1604!$OMP PARALLEL PRIVATE(VALPIV,I,II) 1605 DO I = 1, ELSIZE 1606 VALPIV = ONE/A(POSELT+(LDA8+1_8)*int(IBEG_BLOCK+I-2,8)) 1607!$OMP DO 1608 DO II = 1,NEL1 1609 A(UPOS+int(I-1,8)*LDA8 + int(II-1,8)) = 1610 & A(LPOS+int(I-1,8) + int(II-1,8)*LDA8) 1611 A(LPOS+int(I-1,8) + int(II-1,8)*LDA8) = 1612 & A(LPOS+int(I-1,8) + int(II-1,8)*LDA8)*VALPIV 1613 ENDDO 1614!$OMP END DO NOWAIT 1615 ENDDO 1616!$OMP END PARALLEL 1617 ENDIF 1618 IF ( LAST_VAR - IEND_BLOCK > KEEP(7) ) THEN 1619 BLSIZE = KEEP(8) 1620 ELSE 1621 BLSIZE = LAST_VAR - IEND_BLOCK 1622 END IF 1623 IF ( NASS - IEND_BLOCK .GT. 0 ) THEN 1624#if defined(SAK_BYROW) 1625 DO IROW = IEND_BLOCK+1, LAST_VAR, BLSIZE 1626 Block = min( BLSIZE, NASS - IROW + 1 ) 1627 LPOS = POSELT + int(IROW - 1,8) * LDA8 + 1628 & int(IBEG_BLOCK - 1,8) 1629 UPOS = POSELT + int(IBEG_BLOCK - 1,8) * LDA8 + 1630 & int(IROW - 1,8) 1631 APOS = POSELT + int(IROW - 1,8) * LDA8 + 1632 & int(IEND_BLOCK,8) 1633 CALL zgemm( 'N','N', IROW + Block - IEND_BLOCK - 1, 1634 & Block, NPIV_BLOCK, 1635 & ALPHA, A( UPOS ), LDA, 1636 & A( LPOS ), LDA, ONE, A( APOS ), LDA ) 1637 ENDDO 1638#else 1639 DO IROW = IEND_BLOCK+1, LAST_VAR, BLSIZE 1640 Block = min( BLSIZE, LAST_VAR - IROW + 1 ) 1641 LPOS = POSELT + int( IROW - 1,8) * LDA8 + 1642 & int(IBEG_BLOCK - 1,8) 1643 UPOS = POSELT + int(IBEG_BLOCK - 1,8) * LDA8 + 1644 & int( IROW - 1,8) 1645 APOS = POSELT + int( IROW - 1,8) * LDA8 + int( IROW - 1,8) 1646 CALL zgemm( 'N','N', Block, LAST_VAR - IROW + 1, NPIV_BLOCK, 1647 & ALPHA, A( UPOS ), LDA, 1648 & A( LPOS ), LDA, ONE, A( APOS ), LDA ) 1649 END DO 1650#endif 1651 END IF 1652 LPOS = POSELT + int(LAST_VAR,8)*LDA8 + int(IBEG_BLOCK - 1,8) 1653 UPOS = POSELT + int(IBEG_BLOCK-1,8) * LDA8 + 1654 & int(IEND_BLOCK,8) 1655 APOS = POSELT + int(LAST_VAR,8)*LDA8 + int(IEND_BLOCK,8) 1656 IF (PIVOT_OPTION.EQ.3) THEN 1657 CALL zgemm('N', 'N', NEL1, NFRONT-LAST_VAR, NPIV_BLOCK, ALPHA, 1658 & A(UPOS), LDA, A(LPOS), LDA, ONE, 1659 & A(APOS), LDA) 1660 ELSEIF (PIVOT_OPTION.EQ.2.AND.(NASS.GT. LAST_VAR)) THEN 1661 CALL zgemm('N', 'N', NEL1, NASS-LAST_VAR, NPIV_BLOCK, ALPHA, 1662 & A(UPOS), LDA, A(LPOS), LDA, ONE, 1663 & A(APOS), LDA) 1664 ENDIF 1665 ENDIF 1666 500 CONTINUE 1667 RETURN 1668 END SUBROUTINE ZMUMPS_FAC_SQ_LDLT 1669 SUBROUTINE ZMUMPS_SWAP_LDLT( A, LA, IW, LIW, 1670 & IOLDPS, NPIVP1, IPIV, POSELT, NASS, 1671 & LDA, NFRONT, LEVEL, K219, K50, XSIZE, 1672 & IBEG_BLOCK_TO_SEND ) 1673 IMPLICIT NONE 1674 INTEGER(8) :: POSELT, LA 1675 INTEGER LIW, IOLDPS, NPIVP1, IPIV 1676 INTEGER LDA, NFRONT, NASS, LEVEL, K219, K50, XSIZE 1677 COMPLEX(kind=8) A( LA ) 1678 INTEGER IW( LIW ) 1679 INTEGER, INTENT(IN) :: IBEG_BLOCK_TO_SEND 1680 INCLUDE 'mumps_headers.h' 1681 INTEGER :: LASTROW2SWAP, IBEG 1682 INTEGER ISW, ISWPS1, ISWPS2, HF 1683 INTEGER(8) :: IDIAG, APOS 1684 INTEGER(8) :: LDA8 1685 COMPLEX(kind=8) SWOP 1686 LDA8 = int(LDA,8) 1687 APOS = POSELT + LDA8*int(IPIV-1,8) + int(NPIVP1-1,8) 1688 IDIAG = APOS + int(IPIV - NPIVP1,8) 1689 HF = 6 + IW( IOLDPS + 5 + XSIZE) + XSIZE 1690 ISWPS1 = IOLDPS + HF + NPIVP1 - 1 1691 ISWPS2 = IOLDPS + HF + IPIV - 1 1692 ISW = IW(ISWPS1) 1693 IW(ISWPS1) = IW(ISWPS2) 1694 IW(ISWPS2) = ISW 1695 ISW = IW(ISWPS1+NFRONT) 1696 IW(ISWPS1+NFRONT) = IW(ISWPS2+NFRONT) 1697 IW(ISWPS2+NFRONT) = ISW 1698 IF ( LEVEL .eq. 2 ) THEN 1699 IBEG = IBEG_BLOCK_TO_SEND 1700 CALL zswap( NPIVP1 - 1 - IBEG + 1, 1701 & A( POSELT + int(NPIVP1-1,8) + 1702 & int(IBEG-1,8) * LDA8), LDA, 1703 & A( POSELT + int(IPIV-1,8) + 1704 & int(IBEG-1,8) * LDA8), LDA ) 1705 END IF 1706 CALL zswap( NPIVP1-1, 1707 & A( POSELT+int(NPIVP1-1,8) * LDA8 ), 1, 1708 & A( POSELT + int(IPIV-1,8) * LDA8 ), 1 ) 1709 CALL zswap( IPIV - NPIVP1 - 1, 1710 & A( POSELT+int(NPIVP1,8) * LDA8 + int(NPIVP1-1,8) ), 1711 & LDA, A( APOS + 1_8 ), 1 ) 1712 SWOP = A(IDIAG) 1713 A(IDIAG) = A( POSELT+int(NPIVP1-1,8)*LDA8+int(NPIVP1-1,8) ) 1714 A( POSELT + int(NPIVP1-1,8)*LDA8 + int(NPIVP1-1,8) ) = SWOP 1715 IF (LEVEL .EQ. 1) THEN 1716 LASTROW2SWAP = NFRONT 1717 ELSE 1718 LASTROW2SWAP = NASS 1719 ENDIF 1720 CALL zswap( LASTROW2SWAP - IPIV, 1721 & A( APOS + LDA8 ), LDA, 1722 & A( IDIAG + LDA8 ), LDA ) 1723 IF (K219.NE.0 .AND.K50.EQ.2) THEN 1724 IF ( LEVEL .eq. 2) THEN 1725 APOS = POSELT+LDA8*LDA8-1_8 1726 SWOP = A(APOS+int(NPIVP1,8)) 1727 A(APOS+int(NPIVP1,8))= A(APOS+int(IPIV,8)) 1728 A(APOS+int(IPIV,8)) = SWOP 1729 ENDIF 1730 ENDIF 1731 RETURN 1732 END SUBROUTINE ZMUMPS_SWAP_LDLT 1733 SUBROUTINE ZMUMPS_FAC_T_LDLT(NFRONT,NASS, 1734 & IW,LIW,A,LA, 1735 & LDA, 1736 & IOLDPS,POSELT,KEEP,KEEP8, 1737 & POSTPONE_COL_UPDATE, ETATASS, 1738 & TYPEFile, LAFAC, MonBloc, NextPiv2beWritten, 1739 & LIWFAC, MYID, IFLAG, OFFSET_IW) 1740 USE ZMUMPS_OOC 1741 IMPLICIT NONE 1742 INTEGER NFRONT, NASS,LIW 1743 INTEGER(8) :: LA 1744 COMPLEX(kind=8) A(LA) 1745 INTEGER IW(LIW) 1746 INTEGER KEEP(500) 1747 INTEGER(8) KEEP8(150) 1748 INTEGER(8) :: POSELT 1749 INTEGER LDA 1750 INTEGER IOLDPS, ETATASS 1751 LOGICAL POSTPONE_COL_UPDATE 1752 INTEGER(8) :: LAFAC 1753 INTEGER TYPEFile, NextPiv2beWritten 1754 INTEGER LIWFAC, MYID, IFLAG 1755 TYPE(IO_BLOCK):: MonBloc 1756 INTEGER IDUMMY 1757 LOGICAL LAST_CALL 1758 INTEGER :: OFFSET_IW 1759 INCLUDE 'mumps_headers.h' 1760 INTEGER(8) :: UPOS, APOS, LPOS 1761 INTEGER(8) :: DPOS, POSPV1, POSPV2, OFFDAG 1762 INTEGER(8) :: LDA8 1763 INTEGER BLSIZE, BLSIZE2, Block, IROW, NPIV, I, J, IROWEND 1764 INTEGER I2, I2END, Block2 1765 COMPLEX(kind=8) ONE, ALPHA, BETA, ZERO 1766 COMPLEX(kind=8) :: MULT1, MULT2, A11, DETPIV, A22, A12 1767 PARAMETER (ONE=(1.0D0,0.0D0), ALPHA=(-1.0D0,0.0D0)) 1768 PARAMETER (ZERO=(0.0D0,0.0D0)) 1769 LDA8 = int(LDA,8) 1770 IF (ETATASS.EQ.1) THEN 1771 BETA = ZERO 1772 ELSE 1773 BETA = ONE 1774 ENDIF 1775 IF ( NFRONT - NASS > KEEP(57) ) THEN 1776 BLSIZE = KEEP(58) 1777 ELSE 1778 BLSIZE = NFRONT - NASS 1779 END IF 1780 BLSIZE2 = KEEP(218) 1781 NPIV = IW( IOLDPS + 1 + KEEP(IXSZ)) 1782 IF ( NFRONT - NASS .GT. 0 ) THEN 1783 IF ( POSTPONE_COL_UPDATE ) THEN 1784 CALL ztrsm( 'L', 'U', 'T', 'U', 1785 & NPIV, NFRONT-NPIV, ONE, 1786 & A( POSELT ), LDA, 1787 & A( POSELT + LDA8 * int(NPIV,8) ), LDA ) 1788 ENDIF 1789 DO IROWEND = NFRONT - NASS, 1, -BLSIZE 1790 Block = min( BLSIZE, IROWEND ) 1791 IROW = IROWEND - Block + 1 1792 LPOS = POSELT + int(NASS,8)*LDA8 + int(IROW-1,8) * LDA8 1793 APOS = POSELT + int(NASS,8)*LDA8 + int(IROW-1,8) * LDA8 + 1794 & int(NASS + IROW - 1,8) 1795 UPOS = POSELT + int(NASS,8) 1796 IF (.NOT. POSTPONE_COL_UPDATE) THEN 1797 UPOS = POSELT + int(NASS + IROW - 1,8) 1798 ENDIF 1799 IF (POSTPONE_COL_UPDATE) THEN 1800 DPOS = POSELT 1801 I = 1 1802 DO 1803 IF(I .GT. NPIV) EXIT 1804 IF(IW(OFFSET_IW+I-1) .GT. 0) THEN 1805 A11 = ONE/A(DPOS) 1806 CALL zcopy(Block, A(LPOS+int(I-1,8)), LDA, 1807 & A(UPOS+int(I-1,8)*LDA8), 1) 1808 CALL zscal(Block, A11, A(LPOS+int(I-1,8)), LDA) 1809 DPOS = DPOS + int(LDA+1,8) 1810 I = I+1 1811 ELSE 1812 CALL zcopy(Block, A(LPOS+int(I-1,8)), LDA, 1813 & A(UPOS+int(I-1,8)*LDA8), 1) 1814 CALL zcopy(Block, A(LPOS+int(I,8)), LDA, 1815 & A(UPOS+int(I,8)*LDA8), 1) 1816 POSPV1 = DPOS 1817 POSPV2 = DPOS + int(LDA+1,8) 1818 OFFDAG = POSPV1+1_8 1819 A11 = A(POSPV1) 1820 A22 = A(POSPV2) 1821 A12 = A(OFFDAG) 1822 DETPIV = A11*A22 - A12**2 1823 A22 = A11/DETPIV 1824 A11 = A(POSPV2)/DETPIV 1825 A12 = -A12/DETPIV 1826 DO J = 1,Block 1827 MULT1 = A11*A(LPOS+int(J-1,8)*LDA8+int(I-1,8)) 1828 & + A12*A(LPOS+int(J-1,8)*LDA8+int(I,8)) 1829 MULT2 = A12*A(LPOS+int(J-1,8)*LDA8+int(I-1,8)) 1830 & + A22*A(LPOS+int(J-1,8)*LDA8+int(I,8)) 1831 A(LPOS+int(J-1,8)*LDA8+int(I-1,8)) = MULT1 1832 A(LPOS+int(J-1,8)*LDA8+int(I,8)) = MULT2 1833 ENDDO 1834 DPOS = POSPV2 + int(LDA+1,8) 1835 I = I+2 1836 ENDIF 1837 ENDDO 1838 ENDIF 1839 DO I2END = Block, 1, -BLSIZE2 1840 Block2 = min(BLSIZE2, I2END) 1841 I2 = I2END - Block2+1 1842 CALL zgemm('N', 'N', Block2, Block-I2+1, NPIV, ALPHA, 1843 & A(UPOS+int(I2-1,8)), LDA, 1844 & A(LPOS+int(I2-1,8)*LDA8), LDA, 1845 & BETA, 1846 & A(APOS + int(I2-1,8) + int(I2-1,8)*LDA8), LDA) 1847 IF (KEEP(201).EQ.1) THEN 1848 IF (NextPiv2beWritten.LE.NPIV) THEN 1849 LAST_CALL=.FALSE. 1850 CALL ZMUMPS_OOC_IO_LU_PANEL( 1851 & STRAT_TRY_WRITE, TYPEFile, 1852 & A(POSELT), LAFAC, MonBloc, 1853 & NextPiv2beWritten, IDUMMY, 1854 & IW(IOLDPS), LIWFAC, MYID, 1855 & KEEP8(31), 1856 & IFLAG,LAST_CALL ) 1857 IF (IFLAG .LT. 0 ) RETURN 1858 ENDIF 1859 ENDIF 1860 ENDDO 1861 IF ( NFRONT - NASS - IROW + 1 - Block > 0 ) THEN 1862 CALL zgemm( 'N', 'N', Block, NFRONT-NASS-Block-IROW+1, NPIV, 1863 & ALPHA, A( UPOS ), LDA, 1864 & A( LPOS + LDA8 * int(Block,8) ), LDA, 1865 & BETA, 1866 & A( APOS + LDA8 * int(Block,8) ), LDA ) 1867 ENDIF 1868 END DO 1869 END IF 1870 RETURN 1871 END SUBROUTINE ZMUMPS_FAC_T_LDLT 1872 SUBROUTINE ZMUMPS_STORE_PERMINFO( PIVRPTR, NBPANELS, PIVR, NASS, 1873 & K, P, LastPanelonDisk, 1874 & LastPIVRPTRIndexFilled ) 1875 IMPLICIT NONE 1876 INTEGER, intent(in) :: NBPANELS, NASS, K, P 1877 INTEGER, intent(inout) :: PIVRPTR(NBPANELS), PIVR(NASS) 1878 INTEGER LastPanelonDisk, LastPIVRPTRIndexFilled 1879 INTEGER I 1880 IF ( LastPanelonDisk+1 > NBPANELS ) THEN 1881 WRITE(*,*) "INTERNAL ERROR IN ZMUMPS_STORE_PERMINFO!" 1882 WRITE(*,*) "NASS=",NASS,"PIVRPTR=",PIVRPTR(1:NBPANELS) 1883 WRITE(*,*) "K=",K, "P=",P, "LastPanelonDisk=",LastPanelonDisk 1884 WRITE(*,*) "LastPIVRPTRIndexFilled=", LastPIVRPTRIndexFilled 1885 CALL MUMPS_ABORT() 1886 ENDIF 1887 PIVRPTR(LastPanelonDisk+1) = K + 1 1888 IF (LastPanelonDisk.NE.0) THEN 1889 PIVR(K - PIVRPTR(1) + 1) = P 1890 DO I = LastPIVRPTRIndexFilled + 1, LastPanelonDisk 1891 PIVRPTR(I)=PIVRPTR(LastPIVRPTRIndexFilled) 1892 ENDDO 1893 ENDIF 1894 LastPIVRPTRIndexFilled = LastPanelonDisk + 1 1895 RETURN 1896 END SUBROUTINE ZMUMPS_STORE_PERMINFO 1897 END MODULE ZMUMPS_FAC_FRONT_AUX_M 1898