1102
1104 IMPLICIT NONE
1105 INTEGER(8) :: POSELT, LA
1106 INTEGER NFRONT,NASS,LIW,INODE,IFLAG,INOPV,
1107 & IOLDPS
1108 INTEGER, intent(inout) :: NNEGW, NB22T1W, NBTINYW
1109 INTEGER, intent(inout) :: DET_EXPW, DET_SIGNW
1110 REAL, intent(inout) :: DET_MANTW
1111 INTEGER, intent(in) :: IBEG_BLOCK, IEND_BLOCK
1112 INTEGER, intent(in) :: PIVOT_OPTION,IEND_BLR
1113 INTEGER, intent(inout) :: Inextpiv
1114 LOGICAL, intent(in) :: OOC_EFFECTIVE_ON_FRONT
1115 INTEGER PIVSIZ,LPIV, XSIZE
1116 REAL A(LA)
1117 REAL UU, UULOC, SEUIL
1118 INTEGER IW(LIW)
1119 INTEGER KEEP(500)
1120 INTEGER(8) KEEP8(150)
1121 INTEGER LPN_LIST
1122 INTEGER PIVNUL_LIST(LPN_LIST)
1123 REAL DKEEP(230)
1124 INTEGER PP_FIRST2SWAP_L, PP_LastPanelonDisk
1125 INTEGER PP_LastPIVRPTRIndexFilled
1126 REAL, intent(in) :: MAXFROMM
1127 LOGICAL, intent(inout) :: IS_MAXFROMM_AVAIL
1128 INTEGER, intent(in) :: NVSCHUR
1129 INTEGER, intent(in) :: PARPIV_T1
1130 LOGICAL, intent(in) :: LR_ACTIVATED
1131 include 'mpif.h'
1132 INTEGER (8) :: POSPV1,POSPV2,OFFDAG,APOSJ
1133 INTEGER JMAX, LIM, LIM_SWAP
1134 REAL RMAX,AMAX,TMAX, MAX_PREV_in_PARPIV, ABS_PIVOT
1135 REAL RMAX_NORELAX, TMAX_NORELAX, UULOCM1
1136 INTEGER(8) :: APOSMAX, APOSROW
1137 REAL MAXPIV
1138 REAL PIVNUL
1139 REAL MAXFROMM_UPDATED
1140 REAL FIXA, CSEUIL
1141 REAL PIVOT,DETPIV
1142 REAL ABSDETPIV
1143 include 'mumps_headers.h'
1144 INTEGER :: HF, IPIVNUL
1145 INTEGER :: J
1146 INTEGER(8) :: APOS, J1, J2, JJ, NFRONT8, KK, J1_ini, JJ_ini
1147 INTEGER :: LDA
1148 INTEGER(8) :: LDA8
1149 INTEGER NPIV,IPIV
1150 INTEGER NPIVP1,K
1151 INTEGER :: ISHIFT, K206, IPIV_SHIFT, IPIV_END
1153 REAL ZERO, ONE
1154 parameter( zero = 0.0e0 )
1155 parameter( one = 1.0e0 )
1156 REAL RZERO,RONE
1157 parameter(rzero=0.0e0, rone=1.0e0)
1158#if defined(_OPENMP)
1159 LOGICAL :: OMP_FLAG
1160 INTEGER :: NOMP, CHUNK, J1_end
1161#endif
1162 INTEGER I_PIVRPTR, I_PIVR, NBPANELS_L
1163
1164 pivnul = dkeep(1)
1165 fixa = dkeep(2)
1166 cseuil = seuil
1167 lda = nfront
1168 lda8 = int(lda,8)
1169 nfront8 = int(nfront,8)
1170 k206 = keep(206)
1171 uuloc = uu
1172 IF (uuloc.GT.rzero) THEN
1173 uulocm1 = rone/uuloc
1174 ELSE
1175 uulocm1 = rone
1176 ENDIF
1177 hf = 6 + xsize
1178 IF (keep(50).NE.1 .AND. ooc_effective_on_front) THEN
1180 & i_pivrptr, i_pivr, ioldps+2*nfront+6+keep(ixsz),
1181 & iw, liw)
1182 ENDIF
1183 pivsiz = 1
1184 npiv = iw(ioldps+1+xsize)
1185 npivp1 = npiv + 1
1186 aposmax = poselt+lda8*lda8-1_8
1187 IF(inopv .EQ. -1) THEN
1188 apos = poselt + (lda8+1_8) * int(npiv,8)
1189 pospv1 = apos
1190 CALL smumps_update_minmax_pivot
1191 & ( abs(a(apos)), dkeep, keep, .true.)
1192 IF(abs(a(apos)).LT.seuil) THEN
1193 IF(real(a(apos)) .GE. rzero) THEN
1194 a(apos) = cseuil
1195 ELSE
1196 a(apos) = -cseuil
1197 nnegw = nnegw+1
1198 ENDIF
1199 nbtinyw = nbtinyw + 1
1200 ELSE
1201 IF (keep(258) .NE. 0) THEN
1203 ENDIF
1204 ENDIF
1205 IF (keep(50).NE.1 .AND. ooc_effective_on_front) THEN
1206 CALL smumps_store_perminfo( iw(i_pivrptr), nbpanels_l,
1207 & iw(i_pivr), nass, npivp1, npivp1,
1208 & pp_lastpanelondisk,
1209 & pp_lastpivrptrindexfilled)
1210 ENDIF
1211 GO TO 420
1212 ENDIF
1213 inopv = 0
1214 ishift = 0
1215 ipiv_end = iend_block
1216 IF (k206.GE.1) THEN
1217 IF (inextpiv.GT.npivp1.AND.inextpiv.LE.iend_block) THEN
1218 ishift = inextpiv - npivp1
1219 ENDIF
1220 IF ( k206.EQ.1
1221 & .OR. (k206 .GT.1 .AND. iend_block.EQ.iend_blr) ) THEN
1222 ipiv_end = iend_block + ishift
1223 ENDIF
1224 IF (ishift.GT.0.AND.is_maxfromm_avail) THEN
1225 ipiv = npivp1
1226 apos = poselt + lda8*int(ipiv-1,8) + int(npiv,8)
1227 pospv1 = apos + int(ipiv - npivp1,8)
1228 pivot = a(pospv1)
1229 IF ( maxfromm .GT. pivnul ) THEN
1230 IF (parpiv_t1.NE.0) THEN
1231 maxfromm_updated =
max
1232 & ( maxfromm,
1233 & abs(real(a(aposmax+int(ipiv,8))))
1234 & )
1235 ELSE
1236 maxfromm_updated = maxfromm
1237 ENDIF
1238 IF ( (abs(pivot) .GE. uuloc*maxfromm_updated).AND.
1239 & abs(pivot) .GT.
max(seuil,tiny(maxfromm_updated))
1240 & ) THEN
1241 ishift = 0
1242 ENDIF
1243 ENDIF
1244 ENDIF
1245 IF ( ishift .GT. 0) THEN
1246 is_maxfromm_avail = .false.
1247 ENDIF
1248 ENDIF
1249 DO 460 ipiv_shift = npivp1+ishift, ipiv_end
1250 IF (ipiv_shift .LE. iend_block) THEN
1251 ipiv=ipiv_shift
1252 ELSE
1253 ipiv = ipiv_shift-iend_block-1+npivp1
1254 IF (ibeg_block.EQ.npivp1) THEN
1255 EXIT
1256 ENDIF
1257 ENDIF
1258 apos = poselt + lda8*int(ipiv-1,8) + int(npiv,8)
1259 pospv1 = apos + int(ipiv - npivp1,8)
1260 pivot = a(pospv1)
1261 abs_pivot = abs(pivot)
1262 IF (uuloc.EQ.rzero.OR.pivot_option.EQ.0) THEN
1263 IF(abs_pivot.LT.seuil) THEN
1264 CALL smumps_update_minmax_pivot
1265 & ( abs_pivot, dkeep, keep, .true.)
1266 IF(real(pivot) .GE. rzero) THEN
1267 a(pospv1) = cseuil
1268 ELSE
1269 a(pospv1) = -cseuil
1270 nnegw = nnegw+1
1271 ENDIF
1272 nbtinyw = nbtinyw + 1
1273 ELSE IF (abs_pivot.EQ.rzero) THEN
1274 GO TO 630
1275 ELSE
1276 IF (pivot.LT.rzero) nnegw = nnegw+1
1277 CALL smumps_update_minmax_pivot
1278 & ( abs_pivot, dkeep, keep, .false.)
1279 IF (keep(258) .NE. 0) THEN
1281 ENDIF
1282 ENDIF
1283 GO TO 420
1284 ENDIF
1285 IF ( is_maxfromm_avail ) THEN
1286 IF ( maxfromm .GT. pivnul ) THEN
1287 IF (parpiv_t1.NE.0) THEN
1288 maxfromm_updated =
max
1289 & ( maxfromm,
1290 & abs(real(a(aposmax+int(ipiv,8))))
1291 & )
1292 ELSE
1293 maxfromm_updated = maxfromm
1294 ENDIF
1295 IF ( (abs_pivot .GE. uuloc*maxfromm_updated).AND.
1296 & (abs_pivot .GT.
max(seuil,tiny(maxfromm_updated)))
1297 & ) THEN
1298 IF (pivot .LT. rzero) nnegw = nnegw+1
1299 CALL smumps_update_minmax_pivot
1300 & ( abs_pivot,
1301 & dkeep, keep, .false.)
1302 IF (keep(258) .NE. 0) THEN
1304 ENDIF
1305 GOTO 415
1306 ENDIF
1307 ENDIF
1308 is_maxfromm_avail = .false.
1309 ENDIF
1310 amax = -rone
1311 jmax = 0
1312 IF (pivot_option.EQ.3
1313 & ) THEN
1314 lim = nfront - keep(253)-nvschur
1315 ELSEIF (pivot_option.GE.2
1316 & ) THEN
1317 lim = nass
1318 ELSEIF (pivot_option.GE.1) THEN
1319 lim = iend_blr
1320 ELSE
1321 write(*,*) 'Internal error in FAC_I_LDLT 1x1:',
1322 & pivot_option
1324 ENDIF
1325 j1 = apos
1326 j2 = pospv1 - 1_8
1327 DO jj=j1,j2
1328 IF(abs(a(jj)) .GT. amax) THEN
1329 amax = abs(a(jj))
1330 jmax = ipiv - int(pospv1-jj)
1331 ENDIF
1332 ENDDO
1333 j1 = pospv1 + lda8
1334 DO j=1, iend_block - ipiv
1335 IF(abs(a(j1)) .GT. amax) THEN
1336 amax = abs(a(j1))
1337 jmax = ipiv + j
1338 ENDIF
1339 j1 = j1 + lda8
1340 ENDDO
1341 rmax = rzero
1342 j1_ini = j1
1343#if defined(_OPENMP)
1344 j1_end = lim - iend_block
1345 chunk =
max(j1_end,1)
1346 IF ( j1_end.GE.keep(360)) THEN
1347 omp_flag = .true.
1348 chunk =
max(keep(360)/2,(j1_end+nomp-1)/nomp)
1349 ELSE
1350 omp_flag = .false.
1351 ENDIF
1352#endif
1353
1354
1355 DO j=1, lim - iend_block
1356 j1 = j1_ini + int(j-1,8) * lda8
1357 rmax =
max(abs(a(j1)),rmax)
1358 ENDDO
1359
1360 IF (parpiv_t1.NE.0) THEN
1361 rmax_norelax = real(a(aposmax+int(ipiv,8)))
1362 ELSE
1363 rmax_norelax = rzero
1364 ENDIF
1365 rmax =
max(rmax,rmax_norelax)
1366 IF (
max(amax,rmax,abs(pivot)).LE.pivnul)
THEN
1367 IF ((parpiv_t1.NE.0)
1368 & .AND.(parpiv_t1.NE.-1)
1369 & .AND.(rmax_norelax.LT.0)
1370 & .AND.(ipiv.GT.1)) THEN
1371 max_prev_in_parpiv = rzero
1372 DO jj=1,ipiv-1
1373 max_prev_in_parpiv=
max( max_prev_in_parpiv,
1374 & real(a(aposmax+int(jj,8))) )
1375 ENDDO
1376 IF (max_prev_in_parpiv.GT.pivnul) THEN
1377 aposrow = poselt + nfront8*int(ipiv-1,8)
1378 DO jj=1,ipiv-1
1379 IF (abs(a(aposrow+jj-1)).GT.pivnul) THEN
1380 GOTO 460
1381 ENDIF
1382 ENDDO
1383 ENDIF
1384 ENDIF
1385 CALL smumps_update_minmax_pivot
1386 & ( abs(a(pospv1)), dkeep, keep, .true.)
1387
1388 keep(109) = keep(109) + 1
1389 ipivnul = keep(109)
1390
1391 pivnul_list(ipivnul) = iw( ioldps+hf+npiv+ipiv-npivp1 )
1392 IF(real(fixa).GT.rzero) THEN
1393 IF(real(pivot) .GE. rzero) THEN
1394 a(pospv1) = fixa
1395 ELSE
1396 a(pospv1) = -fixa
1397 ENDIF
1398 ELSE
1399 j1 = apos
1400 j2 = pospv1 - 1_8
1401 DO jj=j1,j2
1402 a(jj) = zero
1403 ENDDO
1404 j1 = pospv1 + lda8
1405 DO j=1, iend_block - ipiv
1406 a(j1) = zero
1407 j1 = j1 + lda8
1408 ENDDO
1409 DO j=1,lim - iend_block
1410 a(j1) = zero
1411 j1 = j1 + lda8
1412 ENDDO
1413 a(pospv1) = one
1414 ENDIF
1415 pivot = a(pospv1)
1416 GO TO 415
1417 ENDIF
1418 rmax =
max(rmax,abs(rmax_norelax))
1419 IF ( abs(pivot).GE.uuloc*
max(rmax,amax)
1420 & .AND. abs(pivot) .GT.
max(seuil,tiny(rmax)) )
THEN
1421 IF (pivot .LT. zero) nnegw = nnegw+1
1422 CALL smumps_update_minmax_pivot
1423 & ( abs(pivot),
1424 & dkeep, keep, .false.)
1425 IF (keep(258) .NE.0 ) THEN
1427 ENDIF
1428 GO TO 415
1429 END IF
1430 IF (npivp1.EQ.iend_block) THEN
1431 GOTO 460
1432 ELSE IF (jmax.EQ.0) THEN
1433 GOTO 460
1434 ENDIF
1435 IF (
max(abs(pivot),rmax,amax).LE.tiny(rmax))
THEN
1436 GOTO 460
1437 ENDIF
1438 IF (
1439 & (keep(19).NE.0).AND.(
max(amax,rmax,abs(pivot)).LE.seuil)
1440 & )
1441 & THEN
1442 GO TO 460
1443 ENDIF
1444 IF (rmax.LT.amax) THEN
1445 j1 = apos
1446 j2 = pospv1 - 1_8
1447 DO jj=j1,j2
1448 IF(int(pospv1-jj) .NE. ipiv-jmax) THEN
1449 rmax =
max(rmax,abs(a(jj)))
1450 ENDIF
1451 ENDDO
1452 j1 = pospv1 + lda8
1453 DO j=1,iend_block-ipiv
1454 IF(ipiv+j .NE. jmax) THEN
1455 rmax =
max(abs(a(j1)),rmax)
1456 ENDIF
1457 j1 = j1 + lda8
1458 ENDDO
1459 ENDIF
1460 aposj = poselt + int(jmax-1,8)*lda8 + int(npiv,8)
1461 pospv2 = aposj + int(jmax - npivp1,8)
1462 IF (ipiv.LT.jmax) THEN
1463 offdag = aposj + int(ipiv - npivp1,8)
1464 ELSE
1465 offdag = apos + int(jmax - npivp1,8)
1466 END IF
1467 tmax = rzero
1468#if defined(_OPENMP)
1469 j1_end = lim-jmax
1470 chunk =
max(j1_end,1)
1471 IF (j1_end.GE.keep(360)) THEN
1472 omp_flag = .true.
1473 chunk =
max(keep(360)/2,(j1_end+nomp-1)/nomp)
1474 ELSE
1475 omp_flag = .false.
1476 ENDIF
1477#endif
1478 IF (jmax .LT. ipiv) THEN
1479 jj_ini = pospv2
1480
1481
1482 DO k = 1, lim - jmax
1483 jj = jj_ini+ int(k,8)*nfront8
1484 IF (jmax+k.NE.ipiv) THEN
1485 tmax=
max(tmax,abs(a(jj)))
1486 ENDIF
1487 ENDDO
1488
1489 DO kk = aposj, pospv2-1_8
1490 tmax =
max(tmax,abs(a(kk)))
1491 ENDDO
1492 ELSE
1493 jj_ini = pospv2
1494
1495
1496 DO k = 1, lim-jmax
1497 jj = jj_ini + int(k,8)*nfront8
1498 tmax=
max(tmax,abs(a(jj)))
1499 ENDDO
1500
1501 DO kk = aposj, pospv2 - 1_8
1502 IF (kk.NE.offdag) THEN
1503 tmax =
max(tmax,abs(a(kk)))
1504 ENDIF
1505 ENDDO
1506 ENDIF
1507 IF (parpiv_t1.NE.0) THEN
1508 tmax_norelax =
max(seuil*uulocm1,
1509 & abs(real(a(aposmax+int(jmax,8))))
1510 & )
1511 ELSE
1512 tmax_norelax = seuil*uulocm1
1513 ENDIF
1514 tmax =
max(tmax,tmax_norelax)
1515 detpiv = a(pospv1)*a(pospv2) - a(offdag)**2
1516 absdetpiv = abs(detpiv)
1517 IF (seuil.GT.rzero) THEN
1518 IF (sqrt(absdetpiv) .LE. seuil ) THEN
1519 GOTO 460
1520 ENDIF
1521 ENDIF
1522 maxpiv =
max(abs(a(pospv1)),abs(a(pospv2)))
1523 IF (maxpiv.EQ.rzero) maxpiv = rone
1524 IF ((abs(a(pospv2))*rmax+amax*tmax)*uuloc.GT.
1525 & absdetpiv .OR. (absdetpiv .EQ. rzero) ) THEN
1526 GO TO 460
1527 ENDIF
1528 IF ((abs(a(pospv1))*tmax+amax*rmax)*uuloc.GT.
1529 & absdetpiv .OR. (absdetpiv.EQ. rzero) ) THEN
1530 GO TO 460
1531 ENDIF
1532 CALL smumps_update_minmax_pivot
1533 & ( sqrt(absdetpiv),
1534 & dkeep, keep, .false.)
1535 IF (keep(258) .NE.0 ) THEN
1537 ENDIF
1538 pivsiz = 2
1539 nb22t1w = nb22t1w + 1
1540 IF(detpiv .LT. rzero) THEN
1541 nnegw = nnegw+1
1542 ELSE IF(a(pospv2) .LT. rzero) THEN
1543 nnegw = nnegw+2
1544 ENDIF
1545 415 CONTINUE
1546 IF (k206.GE.1) THEN
1547 inextpiv =
max(npivp1+pivsiz, ipiv+1)
1548 ENDIF
1549 DO k=1,pivsiz
1550 IF (pivsiz .EQ. 2) THEN
1551 IF (k==1) THEN
1552 lpiv =
min(ipiv,jmax)
1553 ELSE
1554 lpiv =
max(ipiv,jmax)
1555 ENDIF
1556 ELSE
1557 lpiv = ipiv
1558 ENDIF
1559 IF (lpiv.EQ.npivp1) GOTO 416
1560 IF (keep(405) .EQ. 0) THEN
1561 keep8(80) = keep8(80)+1
1562 ELSE
1563
1564 keep8(80) = keep8(80)+1
1565
1566 ENDIF
1567 lim_swap = nfront
1568 CALL smumps_swap_ldlt( a, la, iw, liw,
1569 & ioldps, npivp1, lpiv, poselt, lim_swap,
1570 & lda, nfront, 1, parpiv_t1, keep(50),
1571 & keep(ixsz), -9999)
1572 416 CONTINUE
1573 IF (keep(50).NE.1 .AND. ooc_effective_on_front) THEN
1574 CALL smumps_store_perminfo( iw(i_pivrptr), nbpanels_l,
1575 & iw(i_pivr), nass, npivp1, lpiv, pp_lastpanelondisk,
1576 & pp_lastpivrptrindexfilled)
1577 ENDIF
1578 npivp1 = npivp1 + 1
1579 ENDDO
1580 IF(pivsiz .EQ. 2) THEN
1581 a(poselt+(lda8+1_8)*int(npiv,8)+1_8) = detpiv
1582 ENDIF
1583 GOTO 420
1584 460 CONTINUE
1585 IF (k206 .GE. 1) THEN
1586 inextpiv=iend_block+1
1587 ENDIF
1588 IF (iend_block.EQ.nass) THEN
1589 inopv = 1
1590 ELSE
1591 inopv = 2
1592 ENDIF
1593 GO TO 420
1594 630 CONTINUE
1595 pivsiz = 0
1596 iflag = -10
1597 420 CONTINUE
1598 is_maxfromm_avail = .false.
1599 RETURN