OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
pdstebz.f File Reference

Go to the source code of this file.

Functions/Subroutines

subroutine pdstebz (ictxt, range, order, n, vl, vu, il, iu, abstol, d, e, m, nsplit, w, iblock, isplit, work, lwork, iwork, liwork, info)
subroutine pdlaebz (ijob, n, mmax, minp, abstol, reltol, pivmin, d, nval, intvl, intvlct, mout, lsave, ieflag, info)
subroutine pdlaecv (ijob, kf, kl, intvl, intvlct, nval, abstol, reltol)
subroutine pdlapdct (sigma, n, d, pivmin, count)

Function/Subroutine Documentation

◆ pdlaebz()

subroutine pdlaebz ( integer ijob,
integer n,
integer mmax,
integer minp,
double precision abstol,
double precision reltol,
double precision pivmin,
double precision, dimension( * ) d,
integer, dimension( * ) nval,
double precision, dimension( * ) intvl,
integer, dimension( * ) intvlct,
integer mout,
double precision lsave,
integer ieflag,
integer info )

Definition at line 883 of file pdstebz.f.

886*
887* -- ScaLAPACK routine (version 1.7) --
888* University of Tennessee, Knoxville, Oak Ridge National Laboratory,
889* and University of California, Berkeley.
890* November 15, 1997
891*
892*
893* .. Scalar Arguments ..
894 INTEGER IEFLAG, IJOB, INFO, MINP, MMAX, MOUT, N
895 DOUBLE PRECISION ABSTOL, LSAVE, PIVMIN, RELTOL
896* ..
897* .. Array Arguments ..
898 INTEGER INTVLCT( * ), NVAL( * )
899 DOUBLE PRECISION D( * ), INTVL( * )
900* ..
901*
902* Purpose
903* =======
904*
905* PDLAEBZ contains the iteration loop which computes the eigenvalues
906* contained in the input intervals [ INTVL(2*j-1), INTVL(2*j) ] where
907* j = 1,...,MINP. It uses and computes the function N(w), which is
908* the count of eigenvalues of a symmetric tridiagonal matrix less than
909* or equal to its argument w.
910*
911* This is a ScaLAPACK internal subroutine and arguments are not
912* checked for unreasonable values.
913*
914* Arguments
915* =========
916*
917* IJOB (input) INTEGER
918* Specifies the computation done by PDLAEBZ
919* = 0 : Find an interval with desired values of N(w) at the
920* endpoints of the interval.
921* = 1 : Find a floating point number contained in the initial
922* interval with a desired value of N(w).
923* = 2 : Perform bisection iteration to find eigenvalues of T.
924*
925* N (input) INTEGER
926* The order of the tridiagonal matrix T. N >= 1.
927*
928* MMAX (input) INTEGER
929* The maximum number of intervals that may be generated. If
930* more than MMAX intervals are generated, then PDLAEBZ will
931* quit with INFO = MMAX+1.
932*
933* MINP (input) INTEGER
934* The initial number of intervals. MINP <= MMAX.
935*
936* ABSTOL (input) DOUBLE PRECISION
937* The minimum (absolute) width of an interval. When an interval
938* is narrower than ABSTOL, or than RELTOL times the larger (in
939* magnitude) endpoint, then it is considered to be sufficiently
940* small, i.e., converged.
941* This must be at least zero.
942*
943* RELTOL (input) DOUBLE PRECISION
944* The minimum relative width of an interval. When an interval
945* is narrower than ABSTOL, or than RELTOL times the larger (in
946* magnitude) endpoint, then it is considered to be sufficiently
947* small, i.e., converged.
948* Note : This should be at least radix*machine epsilon.
949*
950* PIVMIN (input) DOUBLE PRECISION
951* The minimum absolute of a "pivot" in the "paranoid"
952* implementation of the Sturm sequence loop. This must be at
953* least max_j |e(j)^2| *safe_min, and at least safe_min, where
954* safe_min is at least the smallest number that can divide 1.0
955* without overflow.
956* See PDLAPDCT for the "paranoid" implementation of the Sturm
957* sequence loop.
958*
959* D (input) DOUBLE PRECISION array, dimension (2*N - 1)
960* Contains the diagonals and the squares of the off-diagonal
961* elements of the tridiagonal matrix T. These elements are
962* assumed to be interleaved in memory for better cache
963* performance. The diagonal entries of T are in the entries
964* D(1),D(3),...,D(2*N-1), while the squares of the off-diagonal
965* entries are D(2),D(4),...,D(2*N-2). To avoid overflow, the
966* matrix must be scaled so that its largest entry is no greater
967* than overflow**(1/2) * underflow**(1/4) in absolute value,
968* and for greatest accuracy, it should not be much smaller
969* than that.
970*
971* NVAL (input/output) INTEGER array, dimension (4)
972* If IJOB = 0, the desired values of N(w) are in NVAL(1) and
973* NVAL(2).
974* If IJOB = 1, NVAL(2) is the desired value of N(w).
975* If IJOB = 2, not referenced.
976* This array will, in general, be reordered on output.
977*
978* INTVL (input/output) DOUBLE PRECISION array, dimension (2*MMAX)
979* The endpoints of the intervals. INTVL(2*j-1) is the left
980* endpoint of the j-th interval, and INTVL(2*j) is the right
981* endpoint of the j-th interval. The input intervals will,
982* in general, be modified, split and reordered by the
983* calculation.
984* On input, INTVL contains the MINP input intervals.
985* On output, INTVL contains the converged intervals.
986*
987* INTVLCT (input/output) INTEGER array, dimension (2*MMAX)
988* The counts at the endpoints of the intervals. INTVLCT(2*j-1)
989* is the count at the left endpoint of the j-th interval, i.e.,
990* the function value N(INTVL(2*j-1)), and INTVLCT(2*j) is the
991* count at the right endpoint of the j-th interval.
992* On input, INTVLCT contains the counts at the endpoints of
993* the MINP input intervals.
994* On output, INTVLCT contains the counts at the endpoints of
995* the converged intervals.
996*
997* MOUT (output) INTEGER
998* The number of intervals output.
999*
1000* LSAVE (output) DOUBLE PRECISION
1001* If IJOB = 0 or 2, not referenced.
1002* If IJOB = 1, this is the largest floating point number
1003* encountered which has count N(w) = NVAL(1).
1004*
1005* IEFLAG (input) INTEGER
1006* A flag which indicates whether N(w) should be speeded up by
1007* exploiting IEEE Arithmetic.
1008*
1009* INFO (output) INTEGER
1010* = 0 : All intervals converged.
1011* = 1 - MMAX : The last INFO intervals did not converge.
1012* = MMAX + 1 : More than MMAX intervals were generated.
1013*
1014* =====================================================================
1015*
1016* .. Intrinsic Functions ..
1017 INTRINSIC abs, int, log, max, min
1018* ..
1019* .. External Subroutines ..
1020 EXTERNAL pdlaecv, pdlaiectb, pdlaiectl, pdlapdct
1021* ..
1022* .. Parameters ..
1023 DOUBLE PRECISION ZERO, TWO, HALF
1024 parameter( zero = 0.0d+0, two = 2.0d+0,
1025 $ half = 1.0d+0 / two )
1026* ..
1027* .. Local Scalars ..
1028 INTEGER I, ITMAX, J, K, KF, KL, KLNEW, L, LCNT, LREQ,
1029 $ NALPHA, NBETA, NMID, RCNT, RREQ
1030 DOUBLE PRECISION ALPHA, BETA, MID
1031* ..
1032* .. Executable Statements ..
1033*
1034 kf = 1
1035 kl = minp + 1
1036 info = 0
1037 IF( intvl( 2 )-intvl( 1 ).LE.zero ) THEN
1038 info = minp
1039 mout = kf
1040 RETURN
1041 END IF
1042 IF( ijob.EQ.0 ) THEN
1043*
1044* Check if some input intervals have "converged"
1045*
1046 CALL pdlaecv( 0, kf, kl, intvl, intvlct, nval,
1047 $ max( abstol, pivmin ), reltol )
1048 IF( kf.GE.kl )
1049 $ GO TO 60
1050*
1051* Compute upper bound on number of iterations needed
1052*
1053 itmax = int( ( log( intvl( 2 )-intvl( 1 )+pivmin )-
1054 $ log( pivmin ) ) / log( two ) ) + 2
1055*
1056* Iteration Loop
1057*
1058 DO 20 i = 1, itmax
1059 klnew = kl
1060 DO 10 j = kf, kl - 1
1061 k = 2*j
1062*
1063* Bisect the interval and find the count at that point
1064*
1065 mid = half*( intvl( k-1 )+intvl( k ) )
1066 IF( ieflag.EQ.0 ) THEN
1067 CALL pdlapdct( mid, n, d, pivmin, nmid )
1068 ELSE IF( ieflag.EQ.1 ) THEN
1069 CALL pdlaiectb( mid, n, d, nmid )
1070 ELSE
1071 CALL pdlaiectl( mid, n, d, nmid )
1072 END IF
1073 lreq = nval( k-1 )
1074 rreq = nval( k )
1075 IF( kl.EQ.1 )
1076 $ nmid = min( intvlct( k ),
1077 $ max( intvlct( k-1 ), nmid ) )
1078 IF( nmid.LE.nval( k-1 ) ) THEN
1079 intvl( k-1 ) = mid
1080 intvlct( k-1 ) = nmid
1081 END IF
1082 IF( nmid.GE.nval( k ) ) THEN
1083 intvl( k ) = mid
1084 intvlct( k ) = nmid
1085 END IF
1086 IF( nmid.GT.lreq .AND. nmid.LT.rreq ) THEN
1087 l = 2*klnew
1088 intvl( l-1 ) = mid
1089 intvl( l ) = intvl( k )
1090 intvlct( l-1 ) = nval( k )
1091 intvlct( l ) = intvlct( k )
1092 intvl( k ) = mid
1093 intvlct( k ) = nval( k-1 )
1094 nval( l-1 ) = nval( k )
1095 nval( l ) = nval( l-1 )
1096 nval( k ) = nval( k-1 )
1097 klnew = klnew + 1
1098 END IF
1099 10 CONTINUE
1100 kl = klnew
1101 CALL pdlaecv( 0, kf, kl, intvl, intvlct, nval,
1102 $ max( abstol, pivmin ), reltol )
1103 IF( kf.GE.kl )
1104 $ GO TO 60
1105 20 CONTINUE
1106 ELSE IF( ijob.EQ.1 ) THEN
1107 alpha = intvl( 1 )
1108 beta = intvl( 2 )
1109 nalpha = intvlct( 1 )
1110 nbeta = intvlct( 2 )
1111 lsave = alpha
1112 lreq = nval( 1 )
1113 rreq = nval( 2 )
1114 30 CONTINUE
1115 IF( nbeta.NE.rreq .AND. beta-alpha.GT.
1116 $ max( abstol, reltol*max( abs( alpha ), abs( beta ) ) ) )
1117 $ THEN
1118*
1119* Bisect the interval and find the count at that point
1120*
1121 mid = half*( alpha+beta )
1122 IF( ieflag.EQ.0 ) THEN
1123 CALL pdlapdct( mid, n, d, pivmin, nmid )
1124 ELSE IF( ieflag.EQ.1 ) THEN
1125 CALL pdlaiectb( mid, n, d, nmid )
1126 ELSE
1127 CALL pdlaiectl( mid, n, d, nmid )
1128 END IF
1129 nmid = min( nbeta, max( nalpha, nmid ) )
1130 IF( nmid.GE.rreq ) THEN
1131 beta = mid
1132 nbeta = nmid
1133 ELSE
1134 alpha = mid
1135 nalpha = nmid
1136 IF( nmid.EQ.lreq )
1137 $ lsave = alpha
1138 END IF
1139 GO TO 30
1140 END IF
1141 kl = kf
1142 intvl( 1 ) = alpha
1143 intvl( 2 ) = beta
1144 intvlct( 1 ) = nalpha
1145 intvlct( 2 ) = nbeta
1146 ELSE IF( ijob.EQ.2 ) THEN
1147*
1148* Check if some input intervals have "converged"
1149*
1150 CALL pdlaecv( 1, kf, kl, intvl, intvlct, nval,
1151 $ max( abstol, pivmin ), reltol )
1152 IF( kf.GE.kl )
1153 $ GO TO 60
1154*
1155* Compute upper bound on number of iterations needed
1156*
1157 itmax = int( ( log( intvl( 2 )-intvl( 1 )+pivmin )-
1158 $ log( pivmin ) ) / log( two ) ) + 2
1159*
1160* Iteration Loop
1161*
1162 DO 50 i = 1, itmax
1163 klnew = kl
1164 DO 40 j = kf, kl - 1
1165 k = 2*j
1166 mid = half*( intvl( k-1 )+intvl( k ) )
1167 IF( ieflag.EQ.0 ) THEN
1168 CALL pdlapdct( mid, n, d, pivmin, nmid )
1169 ELSE IF( ieflag.EQ.1 ) THEN
1170 CALL pdlaiectb( mid, n, d, nmid )
1171 ELSE
1172 CALL pdlaiectl( mid, n, d, nmid )
1173 END IF
1174 lcnt = intvlct( k-1 )
1175 rcnt = intvlct( k )
1176 nmid = min( rcnt, max( lcnt, nmid ) )
1177*
1178* Form New Interval(s)
1179*
1180 IF( nmid.EQ.lcnt ) THEN
1181 intvl( k-1 ) = mid
1182 ELSE IF( nmid.EQ.rcnt ) THEN
1183 intvl( k ) = mid
1184 ELSE IF( klnew.LT.mmax+1 ) THEN
1185 l = 2*klnew
1186 intvl( l-1 ) = mid
1187 intvl( l ) = intvl( k )
1188 intvlct( l-1 ) = nmid
1189 intvlct( l ) = intvlct( k )
1190 intvl( k ) = mid
1191 intvlct( k ) = nmid
1192 klnew = klnew + 1
1193 ELSE
1194 info = mmax + 1
1195 RETURN
1196 END IF
1197 40 CONTINUE
1198 kl = klnew
1199 CALL pdlaecv( 1, kf, kl, intvl, intvlct, nval,
1200 $ max( abstol, pivmin ), reltol )
1201 IF( kf.GE.kl )
1202 $ GO TO 60
1203 50 CONTINUE
1204 END IF
1205 60 CONTINUE
1206 info = max( kl-kf, 0 )
1207 mout = kl - 1
1208 RETURN
1209*
1210* End of PDLAEBZ
1211*
alpha
#define alpha
Definition eval.h:35
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
pdlaecv
subroutine pdlaecv(ijob, kf, kl, intvl, intvlct, nval, abstol, reltol)
Definition pdstebz.f:1217
pdlapdct
subroutine pdlapdct(sigma, n, d, pivmin, count)
Definition pdstebz.f:1365

◆ pdlaecv()

subroutine pdlaecv ( integer ijob,
integer kf,
integer kl,
double precision, dimension( * ) intvl,
integer, dimension( * ) intvlct,
integer, dimension( * ) nval,
double precision abstol,
double precision reltol )

Definition at line 1215 of file pdstebz.f.

1217*
1218* -- ScaLAPACK routine (version 1.7) --
1219* University of Tennessee, Knoxville, Oak Ridge National Laboratory,
1220* and University of California, Berkeley.
1221* November 15, 1997
1222*
1223*
1224* .. Scalar Arguments ..
1225 INTEGER IJOB, KF, KL
1226 DOUBLE PRECISION ABSTOL, RELTOL
1227* ..
1228* .. Array Arguments ..
1229 INTEGER INTVLCT( * ), NVAL( * )
1230 DOUBLE PRECISION INTVL( * )
1231* ..
1232*
1233* Purpose
1234* =======
1235*
1236* PDLAECV checks if the input intervals [ INTVL(2*i-1), INTVL(2*i) ],
1237* i = KF, ... , KL-1, have "converged".
1238* PDLAECV modifies KF to be the index of the last converged interval,
1239* i.e., on output, all intervals [ INTVL(2*i-1), INTVL(2*i) ], i < KF,
1240* have converged. Note that the input intervals may be reordered by
1241* PDLAECV.
1242*
1243* This is a SCALAPACK internal procedure and arguments are not checked
1244* for unreasonable values.
1245*
1246* Arguments
1247* =========
1248*
1249* IJOB (input) INTEGER
1250* Specifies the criterion for "convergence" of an interval.
1251* = 0 : When an interval is narrower than ABSTOL, or than
1252* RELTOL times the larger (in magnitude) endpoint, then
1253* it is considered to have "converged".
1254* = 1 : When an interval is narrower than ABSTOL, or than
1255* RELTOL times the larger (in magnitude) endpoint, or if
1256* the counts at the endpoints are identical to the counts
1257* specified by NVAL ( see NVAL ) then the interval is
1258* considered to have "converged".
1259*
1260* KF (input/output) INTEGER
1261* On input, the index of the first input interval is 2*KF-1.
1262* On output, the index of the last converged interval
1263* is 2*KF-3.
1264*
1265* KL (input) INTEGER
1266* The index of the last input interval is 2*KL-3.
1267*
1268* INTVL (input/output) DOUBLE PRECISION array, dimension (2*(KL-KF))
1269* The endpoints of the intervals. INTVL(2*j-1) is the left
1270* oendpoint f the j-th interval, and INTVL(2*j) is the right
1271* endpoint of the j-th interval. The input intervals will,
1272* in general, be reordered on output.
1273* On input, INTVL contains the KL-KF input intervals.
1274* On output, INTVL contains the converged intervals, 1 thru'
1275* KF-1, and the unconverged intervals, KF thru' KL-1.
1276*
1277* INTVLCT (input/output) INTEGER array, dimension (2*(KL-KF))
1278* The counts at the endpoints of the intervals. INTVLCT(2*j-1)
1279* is the count at the left endpoint of the j-th interval, i.e.,
1280* the function value N(INTVL(2*j-1)), and INTVLCT(2*j) is the
1281* count at the right endpoint of the j-th interval. This array
1282* will, in general, be reordered on output.
1283* See the comments in PDLAEBZ for more on the function N(w).
1284*
1285* NVAL (input/output) INTEGER array, dimension (2*(KL-KF))
1286* The desired counts, N(w), at the endpoints of the
1287* corresponding intervals. This array will, in general,
1288* be reordered on output.
1289*
1290* ABSTOL (input) DOUBLE PRECISION
1291* The minimum (absolute) width of an interval. When an interval
1292* is narrower than ABSTOL, or than RELTOL times the larger (in
1293* magnitude) endpoint, then it is considered to be sufficiently
1294* small, i.e., converged.
1295* Note : This must be at least zero.
1296*
1297* RELTOL (input) DOUBLE PRECISION
1298* The minimum relative width of an interval. When an interval
1299* is narrower than ABSTOL, or than RELTOL times the larger (in
1300* magnitude) endpoint, then it is considered to be sufficiently
1301* small, i.e., converged.
1302* Note : This should be at least radix*machine epsilon.
1303*
1304* =====================================================================
1305*
1306* .. Intrinsic Functions ..
1307 INTRINSIC abs, max
1308* ..
1309* .. Local Scalars ..
1310 LOGICAL CONDN
1311 INTEGER I, ITMP1, ITMP2, J, K, KFNEW
1312 DOUBLE PRECISION TMP1, TMP2, TMP3, TMP4
1313* ..
1314* .. Executable Statements ..
1315*
1316 kfnew = kf
1317 DO 10 i = kf, kl - 1
1318 k = 2*i
1319 tmp3 = intvl( k-1 )
1320 tmp4 = intvl( k )
1321 tmp1 = abs( tmp4-tmp3 )
1322 tmp2 = max( abs( tmp3 ), abs( tmp4 ) )
1323 condn = tmp1.LT.max( abstol, reltol*tmp2 )
1324 IF( ijob.EQ.0 )
1325 $ condn = condn .OR. ( ( intvlct( k-1 ).EQ.nval( k-1 ) ) .AND.
1326 $ intvlct( k ).EQ.nval( k ) )
1327 IF( condn ) THEN
1328 IF( i.GT.kfnew ) THEN
1329*
1330* Reorder Intervals
1331*
1332 j = 2*kfnew
1333 tmp1 = intvl( k-1 )
1334 tmp2 = intvl( k )
1335 itmp1 = intvlct( k-1 )
1336 itmp2 = intvlct( k )
1337 intvl( k-1 ) = intvl( j-1 )
1338 intvl( k ) = intvl( j )
1339 intvlct( k-1 ) = intvlct( j-1 )
1340 intvlct( k ) = intvlct( j )
1341 intvl( j-1 ) = tmp1
1342 intvl( j ) = tmp2
1343 intvlct( j-1 ) = itmp1
1344 intvlct( j ) = itmp2
1345 IF( ijob.EQ.0 ) THEN
1346 itmp1 = nval( k-1 )
1347 nval( k-1 ) = nval( j-1 )
1348 nval( j-1 ) = itmp1
1349 itmp1 = nval( k )
1350 nval( k ) = nval( j )
1351 nval( j ) = itmp1
1352 END IF
1353 END IF
1354 kfnew = kfnew + 1
1355 END IF
1356 10 CONTINUE
1357 kf = kfnew
1358 RETURN
1359*
1360* End of PDLAECV
1361*

◆ pdlapdct()

subroutine pdlapdct ( double precision sigma,
integer n,
double precision, dimension( * ) d,
double precision pivmin,
integer count )

Definition at line 1364 of file pdstebz.f.

1365*
1366* -- ScaLAPACK routine (version 1.7) --
1367* University of Tennessee, Knoxville, Oak Ridge National Laboratory,
1368* and University of California, Berkeley.
1369* November 15, 1997
1370*
1371*
1372* .. Scalar Arguments ..
1373 INTEGER COUNT, N
1374 DOUBLE PRECISION PIVMIN, SIGMA
1375* ..
1376* .. Array Arguments ..
1377 DOUBLE PRECISION D( * )
1378* ..
1379*
1380* Purpose
1381* =======
1382*
1383* PDLAPDCT counts the number of negative eigenvalues of (T - SIGMA I).
1384* This implementation of the Sturm Sequence loop has conditionals in
1385* the innermost loop to avoid overflow and determine the sign of a
1386* floating point number. PDLAPDCT will be referred to as the "paranoid"
1387* implementation of the Sturm Sequence loop.
1388*
1389* This is a SCALAPACK internal procedure and arguments are not checked
1390* for unreasonable values.
1391*
1392* Arguments
1393* =========
1394*
1395* SIGMA (input) DOUBLE PRECISION
1396* The shift. PDLAPDCT finds the number of eigenvalues of T less
1397* than or equal to SIGMA.
1398*
1399* N (input) INTEGER
1400* The order of the tridiagonal matrix T. N >= 1.
1401*
1402* D (input) DOUBLE PRECISION array, dimension (2*N - 1)
1403* Contains the diagonals and the squares of the off-diagonal
1404* elements of the tridiagonal matrix T. These elements are
1405* assumed to be interleaved in memory for better cache
1406* performance. The diagonal entries of T are in the entries
1407* D(1),D(3),...,D(2*N-1), while the squares of the off-diagonal
1408* entries are D(2),D(4),...,D(2*N-2). To avoid overflow, the
1409* matrix must be scaled so that its largest entry is no greater
1410* than overflow**(1/2) * underflow**(1/4) in absolute value,
1411* and for greatest accuracy, it should not be much smaller
1412* than that.
1413*
1414* PIVMIN (input) DOUBLE PRECISION
1415* The minimum absolute of a "pivot" in this "paranoid"
1416* implementation of the Sturm sequence loop. This must be at
1417* least max_j |e(j)^2| *safe_min, and at least safe_min, where
1418* safe_min is at least the smallest number that can divide 1.0
1419* without overflow.
1420*
1421* COUNT (output) INTEGER
1422* The count of the number of eigenvalues of T less than or
1423* equal to SIGMA.
1424*
1425* =====================================================================
1426*
1427* .. Intrinsic Functions ..
1428 INTRINSIC abs
1429* ..
1430* .. Parameters ..
1431 DOUBLE PRECISION ZERO
1432 parameter( zero = 0.0d+0 )
1433* ..
1434* .. Local Scalars ..
1435 INTEGER I
1436 DOUBLE PRECISION TMP
1437* ..
1438* .. Executable Statements ..
1439*
1440 tmp = d( 1 ) - sigma
1441 IF( abs( tmp ).LE.pivmin )
1442 $ tmp = -pivmin
1443 count = 0
1444 IF( tmp.LE.zero )
1445 $ count = 1
1446 DO 10 i = 3, 2*n - 1, 2
1447 tmp = d( i ) - d( i-1 ) / tmp - sigma
1448 IF( abs( tmp ).LE.pivmin )
1449 $ tmp = -pivmin
1450 IF( tmp.LE.zero )
1451 $ count = count + 1
1452 10 CONTINUE
1453*
1454 RETURN
1455*
1456* End of PDLAPDCT
1457*

◆ pdstebz()

subroutine pdstebz ( integer ictxt,
character range,
character order,
integer n,
double precision vl,
double precision vu,
integer il,
integer iu,
double precision abstol,
double precision, dimension( * ) d,
double precision, dimension( * ) e,
integer m,
integer nsplit,
double precision, dimension( * ) w,
integer, dimension( * ) iblock,
integer, dimension( * ) isplit,
double precision, dimension( * ) work,
integer lwork,
integer, dimension( * ) iwork,
integer liwork,
integer info )

Definition at line 1 of file pdstebz.f.

4*
5* -- ScaLAPACK routine (version 1.7) --
6* University of Tennessee, Knoxville, Oak Ridge National Laboratory,
7* and University of California, Berkeley.
8* November 15, 1997
9*
10* .. Scalar Arguments ..
11 CHARACTER ORDER, RANGE
12 INTEGER ICTXT, IL, INFO, IU, LIWORK, LWORK, M, N,
13 $ NSPLIT
14 DOUBLE PRECISION ABSTOL, VL, VU
15* ..
16* .. Array Arguments ..
17 INTEGER IBLOCK( * ), ISPLIT( * ), IWORK( * )
18 DOUBLE PRECISION D( * ), E( * ), W( * ), WORK( * )
19* ..
20*
21* Purpose
22* =======
23*
24* PDSTEBZ computes the eigenvalues of a symmetric tridiagonal matrix in
25* parallel. The user may ask for all eigenvalues, all eigenvalues in
26* the interval [VL, VU], or the eigenvalues indexed IL through IU. A
27* static partitioning of work is done at the beginning of PDSTEBZ which
28* results in all processes finding an (almost) equal number of
29* eigenvalues.
30*
31* NOTE : It is assumed that the user is on an IEEE machine. If the user
32* is not on an IEEE mchine, set the compile time flag NO_IEEE
33* to 1 (in SLmake.inc). The features of IEEE arithmetic that
34* are needed for the "fast" Sturm Count are : (a) infinity
35* arithmetic (b) the sign bit of a single precision floating
36* point number is assumed be in the 32nd bit position
37* (c) the sign of negative zero.
38*
39* See W. Kahan "Accurate Eigenvalues of a Symmetric Tridiagonal
40* Matrix", Report CS41, Computer Science Dept., Stanford
41* University, July 21, 1966.
42*
43* Arguments
44* =========
45*
46* ICTXT (global input) INTEGER
47* The BLACS context handle.
48*
49* RANGE (global input) CHARACTER
50* Specifies which eigenvalues are to be found.
51* = 'A': ("All") all eigenvalues will be found.
52* = 'V': ("Value") all eigenvalues in the interval
53* [VL, VU] will be found.
54* = 'I': ("Index") the IL-th through IU-th eigenvalues (of the
55* entire matrix) will be found.
56*
57* ORDER (global input) CHARACTER
58* Specifies the order in which the eigenvalues and their block
59* numbers are stored in W and IBLOCK.
60* = 'B': ("By Block") the eigenvalues will be grouped by
61* split-off block (see IBLOCK, ISPLIT) and
62* ordered from smallest to largest within
63* the block.
64* = 'E': ("Entire matrix")
65* the eigenvalues for the entire matrix
66* will be ordered from smallest to largest.
67*
68* N (global input) INTEGER
69* The order of the tridiagonal matrix T. N >= 0.
70*
71* VL (global input) DOUBLE PRECISION
72* If RANGE='V', the lower bound of the interval to be searched
73* for eigenvalues. Eigenvalues less than VL will not be
74* returned. Not referenced if RANGE='A' or 'I'.
75*
76* VU (global input) DOUBLE PRECISION
77* If RANGE='V', the upper bound of the interval to be searched
78* for eigenvalues. Eigenvalues greater than VU will not be
79* returned. VU must be greater than VL. Not referenced if
80* RANGE='A' or 'I'.
81*
82* IL (global input) INTEGER
83* If RANGE='I', the index (from smallest to largest) of the
84* smallest eigenvalue to be returned. IL must be at least 1.
85* Not referenced if RANGE='A' or 'V'.
86*
87* IU (global input) INTEGER
88* If RANGE='I', the index (from smallest to largest) of the
89* largest eigenvalue to be returned. IU must be at least IL
90* and no greater than N. Not referenced if RANGE='A' or 'V'.
91*
92* ABSTOL (global input) DOUBLE PRECISION
93* The absolute tolerance for the eigenvalues. An eigenvalue
94* (or cluster) is considered to be located if it has been
95* determined to lie in an interval whose width is ABSTOL or
96* less. If ABSTOL is less than or equal to zero, then ULP*|T|
97* will be used, where |T| means the 1-norm of T.
98* Eigenvalues will be computed most accurately when ABSTOL is
99* set to the underflow threshold DLAMCH('U'), not zero.
100* Note : If eigenvectors are desired later by inverse iteration
101* ( PDSTEIN ), ABSTOL should be set to 2*PDLAMCH('S').
102*
103* D (global input) DOUBLE PRECISION array, dimension (N)
104* The n diagonal elements of the tridiagonal matrix T. To
105* avoid overflow, the matrix must be scaled so that its largest
106* entry is no greater than overflow**(1/2) * underflow**(1/4)
107* in absolute value, and for greatest accuracy, it should not
108* be much smaller than that.
109*
110* E (global input) DOUBLE PRECISION array, dimension (N-1)
111* The (n-1) off-diagonal elements of the tridiagonal matrix T.
112* To avoid overflow, the matrix must be scaled so that its
113* largest entry is no greater than overflow**(1/2) *
114* underflow**(1/4) in absolute value, and for greatest
115* accuracy, it should not be much smaller than that.
116*
117* M (global output) INTEGER
118* The actual number of eigenvalues found. 0 <= M <= N.
119* (See also the description of INFO=2)
120*
121* NSPLIT (global output) INTEGER
122* The number of diagonal blocks in the matrix T.
123* 1 <= NSPLIT <= N.
124*
125* W (global output) DOUBLE PRECISION array, dimension (N)
126* On exit, the first M elements of W contain the eigenvalues
127* on all processes.
128*
129* IBLOCK (global output) INTEGER array, dimension (N)
130* At each row/column j where E(j) is zero or small, the
131* matrix T is considered to split into a block diagonal
132* matrix. On exit IBLOCK(i) specifies which block (from 1
133* to the number of blocks) the eigenvalue W(i) belongs to.
134* NOTE: in the (theoretically impossible) event that bisection
135* does not converge for some or all eigenvalues, INFO is set
136* to 1 and the ones for which it did not are identified by a
137* negative block number.
138*
139* ISPLIT (global output) INTEGER array, dimension (N)
140* The splitting points, at which T breaks up into submatrices.
141* The first submatrix consists of rows/columns 1 to ISPLIT(1),
142* the second of rows/columns ISPLIT(1)+1 through ISPLIT(2),
143* etc., and the NSPLIT-th consists of rows/columns
144* ISPLIT(NSPLIT-1)+1 through ISPLIT(NSPLIT)=N.
145* (Only the first NSPLIT elements will actually be used, but
146* since the user cannot know a priori what value NSPLIT will
147* have, N words must be reserved for ISPLIT.)
148*
149* WORK (local workspace) DOUBLE PRECISION array,
150* dimension ( MAX( 5*N, 7 ) )
151*
152* LWORK (local input) INTEGER
153* size of array WORK must be >= MAX( 5*N, 7 )
154* If LWORK = -1, then LWORK is global input and a workspace
155* query is assumed; the routine only calculates the minimum
156* and optimal size for all work arrays. Each of these
157* values is returned in the first entry of the corresponding
158* work array, and no error message is issued by PXERBLA.
159*
160* IWORK (local workspace) INTEGER array, dimension ( MAX( 4*N, 14 ) )
161*
162* LIWORK (local input) INTEGER
163* size of array IWORK must be >= MAX( 4*N, 14, NPROCS )
164* If LIWORK = -1, then LIWORK is global input and a workspace
165* query is assumed; the routine only calculates the minimum
166* and optimal size for all work arrays. Each of these
167* values is returned in the first entry of the corresponding
168* work array, and no error message is issued by PXERBLA.
169*
170* INFO (global output) INTEGER
171* = 0 : successful exit
172* < 0 : if INFO = -i, the i-th argument had an illegal value
173* > 0 : some or all of the eigenvalues failed to converge or
174* were not computed:
175* = 1 : Bisection failed to converge for some eigenvalues;
176* these eigenvalues are flagged by a negative block
177* number. The effect is that the eigenvalues may not
178* be as accurate as the absolute and relative
179* tolerances. This is generally caused by arithmetic
180* which is less accurate than PDLAMCH says.
181* = 2 : There is a mismatch between the number of
182* eigenvalues output and the number desired.
183* = 3 : RANGE='i', and the Gershgorin interval initially
184* used was incorrect. No eigenvalues were computed.
185* Probable cause: your machine has sloppy floating
186* point arithmetic.
187* Cure: Increase the PARAMETER "FUDGE", recompile,
188* and try again.
189*
190* Internal Parameters
191* ===================
192*
193* RELFAC DOUBLE PRECISION, default = 2.0
194* The relative tolerance. An interval [a,b] lies within
195* "relative tolerance" if b-a < RELFAC*ulp*max(|a|,|b|),
196* where "ulp" is the machine precision (distance from 1 to
197* the next larger floating point number.)
198*
199* FUDGE DOUBLE PRECISION, default = 2.0
200* A "fudge factor" to widen the Gershgorin intervals. Ideally,
201* a value of 1 should work, but on machines with sloppy
202* arithmetic, this needs to be larger. The default for
203* publicly released versions should be large enough to handle
204* the worst machine around. Note that this has no effect
205* on the accuracy of the solution.
206*
207* =====================================================================
208*
209* .. Intrinsic Functions ..
210 INTRINSIC abs, dble, ichar, max, min, mod
211* ..
212* .. External Functions ..
213 LOGICAL LSAME
214 INTEGER BLACS_PNUM
215 DOUBLE PRECISION PDLAMCH
216 EXTERNAL lsame, blacs_pnum, pdlamch
217* ..
218* .. External Subroutines ..
219 EXTERNAL blacs_freebuff, blacs_get, blacs_gridexit,
220 $ blacs_gridinfo, blacs_gridmap, dgebr2d,
221 $ dgebs2d, dgerv2d, dgesd2d, dlasrt2, globchk,
222 $ igebr2d, igebs2d, igerv2d, igesd2d, igsum2d,
223 $ pdlaebz, pdlaiectb, pdlaiectl, pdlapdct,
224 $ pdlasnbt, pxerbla
225* ..
226* .. Parameters ..
227 INTEGER BLOCK_CYCLIC_2D, DLEN_, DTYPE_, CTXT_, M_, N_,
228 $ MB_, NB_, RSRC_, CSRC_, LLD_
229 parameter( block_cyclic_2d = 1, dlen_ = 9, dtype_ = 1,
230 $ ctxt_ = 2, m_ = 3, n_ = 4, mb_ = 5, nb_ = 6,
231 $ rsrc_ = 7, csrc_ = 8, lld_ = 9 )
232 INTEGER BIGNUM, DESCMULT
233 parameter( bignum = 10000, descmult = 100 )
234 DOUBLE PRECISION ZERO, ONE, TWO, FIVE, HALF
235 parameter( zero = 0.0d+0, one = 1.0d+0, two = 2.0d+0,
236 $ five = 5.0d+0, half = 1.0d+0 / two )
237 DOUBLE PRECISION FUDGE, RELFAC
238 parameter( fudge = 2.0d+0, relfac = 2.0d+0 )
239* ..
240* .. Local Scalars ..
241 LOGICAL LQUERY
242 INTEGER BLKNO, FOUND, I, IBEGIN, IEFLAG, IEND, IFRST,
243 $ IINFO, ILAST, ILOAD, IM, IMYLOAD, IN, INDRIW1,
244 $ INDRIW2, INDRW1, INDRW2, INXTLOAD, IOFF,
245 $ IORDER, IOUT, IRANGE, IRECV, IREM, ITMP1,
246 $ ITMP2, J, JB, K, LAST, LEXTRA, LREQ, MYCOL,
247 $ MYROW, NALPHA, NBETA, NCMP, NEIGINT, NEXT, NGL,
248 $ NGLOB, NGU, NINT, NPCOL, NPROW, OFFSET,
249 $ ONEDCONTEXT, P, PREV, REXTRA, RREQ, SELF
250 DOUBLE PRECISION ALPHA, ATOLI, BETA, BNORM, DRECV, DSEND, GL,
251 $ GU, INITVL, INITVU, LSAVE, MID, PIVMIN, RELTOL,
252 $ SAFEMN, TMP1, TMP2, TNORM, ULP
253* ..
254* .. Local Arrays ..
255 INTEGER IDUM( 5, 2 )
256 INTEGER TORECV( 1, 1 )
257* ..
258* .. Executable Statements ..
259* This is just to keep ftnchek happy
260 IF( block_cyclic_2d*csrc_*ctxt_*dlen_*dtype_*lld_*mb_*m_*nb_*n_*
261 $ rsrc_.LT.0 )RETURN
262*
263* Set up process grid
264*
265 CALL blacs_gridinfo( ictxt, nprow, npcol, myrow, mycol )
266*
267 info = 0
268 m = 0
269*
270* Decode RANGE
271*
272 IF( lsame( range, 'A' ) ) THEN
273 irange = 1
274 ELSE IF( lsame( range, 'V' ) ) THEN
275 irange = 2
276 ELSE IF( lsame( range, 'I' ) ) THEN
277 irange = 3
278 ELSE
279 irange = 0
280 END IF
281*
282* Decode ORDER
283*
284 IF( lsame( order, 'B' ) ) THEN
285 iorder = 2
286 ELSE IF( lsame( order, 'E' ) .OR. lsame( order, 'A' ) ) THEN
287 iorder = 1
288 ELSE
289 iorder = 0
290 END IF
291*
292* Check for Errors
293*
294 IF( nprow.EQ.-1 ) THEN
295 info = -1
296 ELSE
297*
298* Get machine constants
299*
300 safemn = pdlamch( ictxt, 'S' )
301 ulp = pdlamch( ictxt, 'P' )
302 reltol = ulp*relfac
303 idum( 1, 1 ) = ichar( range )
304 idum( 1, 2 ) = 2
305 idum( 2, 1 ) = ichar( order )
306 idum( 2, 2 ) = 3
307 idum( 3, 1 ) = n
308 idum( 3, 2 ) = 4
309 nglob = 5
310 IF( irange.EQ.3 ) THEN
311 idum( 4, 1 ) = il
312 idum( 4, 2 ) = 7
313 idum( 5, 1 ) = iu
314 idum( 5, 2 ) = 8
315 ELSE
316 idum( 4, 1 ) = 0
317 idum( 4, 2 ) = 0
318 idum( 5, 1 ) = 0
319 idum( 5, 2 ) = 0
320 END IF
321 IF( myrow.EQ.0 .AND. mycol.EQ.0 ) THEN
322 work( 1 ) = abstol
323 IF( irange.EQ.2 ) THEN
324 work( 2 ) = vl
325 work( 3 ) = vu
326 ELSE
327 work( 2 ) = zero
328 work( 3 ) = zero
329 END IF
330 CALL dgebs2d( ictxt, 'ALL', ' ', 3, 1, work, 3 )
331 ELSE
332 CALL dgebr2d( ictxt, 'ALL', ' ', 3, 1, work, 3, 0, 0 )
333 END IF
334 lquery = ( lwork.EQ.-1 .OR. liwork.EQ.-1 )
335 IF( info.EQ.0 ) THEN
336 IF( irange.EQ.0 ) THEN
337 info = -2
338 ELSE IF( iorder.EQ.0 ) THEN
339 info = -3
340 ELSE IF( irange.EQ.2 .AND. vl.GE.vu ) THEN
341 info = -5
342 ELSE IF( irange.EQ.3 .AND. ( il.LT.1 .OR. il.GT.max( 1,
343 $ n ) ) ) THEN
344 info = -6
345 ELSE IF( irange.EQ.3 .AND. ( iu.LT.min( n,
346 $ il ) .OR. iu.GT.n ) ) THEN
347 info = -7
348 ELSE IF( lwork.LT.max( 5*n, 7 ) .AND. .NOT.lquery ) THEN
349 info = -18
350 ELSE IF( liwork.LT.max( 4*n, 14, nprow*npcol ) .AND. .NOT.
351 $ lquery ) THEN
352 info = -20
353 ELSE IF( irange.EQ.2 .AND. ( abs( work( 2 )-vl ).GT.five*
354 $ ulp*abs( vl ) ) ) THEN
355 info = -5
356 ELSE IF( irange.EQ.2 .AND. ( abs( work( 3 )-vu ).GT.five*
357 $ ulp*abs( vu ) ) ) THEN
358 info = -6
359 ELSE IF( abs( work( 1 )-abstol ).GT.five*ulp*abs( abstol ) )
360 $ THEN
361 info = -9
362 END IF
363 END IF
364 IF( info.EQ.0 )
365 $ info = bignum
366 CALL globchk( ictxt, nglob, idum, 5, iwork, info )
367 IF( info.EQ.bignum ) THEN
368 info = 0
369 ELSE IF( mod( info, descmult ).EQ.0 ) THEN
370 info = -info / descmult
371 ELSE
372 info = -info
373 END IF
374 END IF
375 work( 1 ) = dble( max( 5*n, 7 ) )
376 iwork( 1 ) = max( 4*n, 14, nprow*npcol )
377*
378 IF( info.NE.0 ) THEN
379 CALL pxerbla( ictxt, 'PDSTEBZ', -info )
380 RETURN
381 ELSE IF( lwork.EQ.-1 .AND. liwork.EQ.-1 ) THEN
382 RETURN
383 END IF
384*
385*
386* Quick return if possible
387*
388 IF( n.EQ.0 )
389 $ RETURN
390*
391 k = 1
392 DO 20 i = 0, nprow - 1
393 DO 10 j = 0, npcol - 1
394 iwork( k ) = blacs_pnum( ictxt, i, j )
395 k = k + 1
396 10 CONTINUE
397 20 CONTINUE
398*
399 p = nprow*npcol
400 nprow = 1
401 npcol = p
402*
403 CALL blacs_get( ictxt, 10, onedcontext )
404 CALL blacs_gridmap( onedcontext, iwork, nprow, nprow, npcol )
405 CALL blacs_gridinfo( onedcontext, i, j, k, self )
406*
407* Simplifications:
408*
409 IF( irange.EQ.3 .AND. il.EQ.1 .AND. iu.EQ.n )
410 $ irange = 1
411*
412 next = mod( self+1, p )
413 prev = mod( p+self-1, p )
414*
415* Compute squares of off-diagonals, splitting points and pivmin.
416* Interleave diagonals and off-diagonals.
417*
418 indrw1 = max( 2*n, 4 )
419 indrw2 = indrw1 + 2*n
420 indriw1 = max( 2*n, 8 )
421 nsplit = 1
422 work( indrw1+2*n ) = zero
423 pivmin = one
424*
425 DO 30 i = 1, n - 1
426 tmp1 = e( i )**2
427 j = 2*i
428 work( indrw1+j-1 ) = d( i )
429 IF( abs( d( i+1 )*d( i ) )*ulp**2+safemn.GT.tmp1 ) THEN
430 isplit( nsplit ) = i
431 nsplit = nsplit + 1
432 work( indrw1+j ) = zero
433 ELSE
434 work( indrw1+j ) = tmp1
435 pivmin = max( pivmin, tmp1 )
436 END IF
437 30 CONTINUE
438 work( indrw1+2*n-1 ) = d( n )
439 isplit( nsplit ) = n
440 pivmin = pivmin*safemn
441*
442* Compute Gershgorin interval [gl,gu] for entire matrix
443*
444 gu = d( 1 )
445 gl = d( 1 )
446 tmp1 = zero
447*
448 DO 40 i = 1, n - 1
449 tmp2 = abs( e( i ) )
450 gu = max( gu, d( i )+tmp1+tmp2 )
451 gl = min( gl, d( i )-tmp1-tmp2 )
452 tmp1 = tmp2
453 40 CONTINUE
454 gu = max( gu, d( n )+tmp1 )
455 gl = min( gl, d( n )-tmp1 )
456 tnorm = max( abs( gl ), abs( gu ) )
457 gl = gl - fudge*tnorm*ulp*n - fudge*two*pivmin
458 gu = gu + fudge*tnorm*ulp*n + fudge*pivmin
459*
460 IF( abstol.LE.zero ) THEN
461 atoli = ulp*tnorm
462 ELSE
463 atoli = abstol
464 END IF
465*
466* Find out if on an IEEE machine, the sign bit is the
467* 32nd bit (Big Endian) or the 64th bit (Little Endian)
468*
469 IF( irange.EQ.1 .OR. nsplit.EQ.1 ) THEN
470 CALL pdlasnbt( ieflag )
471 ELSE
472 ieflag = 0
473 END IF
474 lextra = 0
475 rextra = 0
476*
477* Form Initial Interval containing desired eigenvalues
478*
479 IF( irange.EQ.1 ) THEN
480 initvl = gl
481 initvu = gu
482 work( 1 ) = gl
483 work( 2 ) = gu
484 iwork( 1 ) = 0
485 iwork( 2 ) = n
486 ifrst = 1
487 ilast = n
488 ELSE IF( irange.EQ.2 ) THEN
489 IF( vl.GT.gl ) THEN
490 IF( ieflag.EQ.0 ) THEN
491 CALL pdlapdct( vl, n, work( indrw1+1 ), pivmin, ifrst )
492 ELSE IF( ieflag.EQ.1 ) THEN
493 CALL pdlaiectb( vl, n, work( indrw1+1 ), ifrst )
494 ELSE
495 CALL pdlaiectl( vl, n, work( indrw1+1 ), ifrst )
496 END IF
497 ifrst = ifrst + 1
498 initvl = vl
499 ELSE
500 initvl = gl
501 ifrst = 1
502 END IF
503 IF( vu.LT.gu ) THEN
504 IF( ieflag.EQ.0 ) THEN
505 CALL pdlapdct( vu, n, work( indrw1+1 ), pivmin, ilast )
506 ELSE IF( ieflag.EQ.1 ) THEN
507 CALL pdlaiectb( vu, n, work( indrw1+1 ), ilast )
508 ELSE
509 CALL pdlaiectl( vu, n, work( indrw1+1 ), ilast )
510 END IF
511 initvu = vu
512 ELSE
513 initvu = gu
514 ilast = n
515 END IF
516 work( 1 ) = initvl
517 work( 2 ) = initvu
518 iwork( 1 ) = ifrst - 1
519 iwork( 2 ) = ilast
520 ELSE IF( irange.EQ.3 ) THEN
521 work( 1 ) = gl
522 work( 2 ) = gu
523 iwork( 1 ) = 0
524 iwork( 2 ) = n
525 iwork( 5 ) = il - 1
526 iwork( 6 ) = iu
527 CALL pdlaebz( 0, n, 2, 1, atoli, reltol, pivmin,
528 $ work( indrw1+1 ), iwork( 5 ), work, iwork, nint,
529 $ lsave, ieflag, iinfo )
530 IF( iinfo.NE.0 ) THEN
531 info = 3
532 GO TO 230
533 END IF
534 IF( nint.GT.1 ) THEN
535 IF( iwork( 5 ).EQ.il-1 ) THEN
536 work( 2 ) = work( 4 )
537 iwork( 2 ) = iwork( 4 )
538 ELSE
539 work( 1 ) = work( 3 )
540 iwork( 1 ) = iwork( 3 )
541 END IF
542 IF( iwork( 1 ).LT.0 .OR. iwork( 1 ).GT.il-1 .OR.
543 $ iwork( 2 ).LE.min( iu-1, iwork( 1 ) ) .OR.
544 $ iwork( 2 ).GT.n ) THEN
545 info = 3
546 GO TO 230
547 END IF
548 END IF
549 lextra = il - 1 - iwork( 1 )
550 rextra = iwork( 2 ) - iu
551 initvl = work( 1 )
552 initvu = work( 2 )
553 ifrst = il
554 ilast = iu
555 END IF
556* NVL = IFRST - 1
557* NVU = ILAST
558 gl = initvl
559 gu = initvu
560 ngl = iwork( 1 )
561 ngu = iwork( 2 )
562 im = 0
563 found = 0
564 indriw2 = indriw1 + ngu - ngl
565 iend = 0
566 IF( ifrst.GT.ilast )
567 $ GO TO 100
568 IF( ifrst.EQ.1 .AND. ilast.EQ.n )
569 $ irange = 1
570*
571* Find Eigenvalues -- Loop Over Blocks
572*
573 DO 90 jb = 1, nsplit
574 ioff = iend
575 ibegin = ioff + 1
576 iend = isplit( jb )
577 in = iend - ioff
578 IF( jb.NE.1 ) THEN
579 IF( irange.NE.1 ) THEN
580 found = im
581*
582* Find total number of eigenvalues found thus far
583*
584 CALL igsum2d( onedcontext, 'All', ' ', 1, 1, found, 1,
585 $ -1, -1 )
586 ELSE
587 found = ioff
588 END IF
589 END IF
590* IF( SELF.GE.P )
591* $ GO TO 30
592 IF( in.NE.n ) THEN
593*
594* Compute Gershgorin interval [gl,gu] for split matrix
595*
596 gu = d( ibegin )
597 gl = d( ibegin )
598 tmp1 = zero
599*
600 DO 50 j = ibegin, iend - 1
601 tmp2 = abs( e( j ) )
602 gu = max( gu, d( j )+tmp1+tmp2 )
603 gl = min( gl, d( j )-tmp1-tmp2 )
604 tmp1 = tmp2
605 50 CONTINUE
606*
607 gu = max( gu, d( iend )+tmp1 )
608 gl = min( gl, d( iend )-tmp1 )
609 bnorm = max( abs( gl ), abs( gu ) )
610 gl = gl - fudge*bnorm*ulp*in - fudge*pivmin
611 gu = gu + fudge*bnorm*ulp*in + fudge*pivmin
612*
613* Compute ATOLI for the current submatrix
614*
615 IF( abstol.LE.zero ) THEN
616 atoli = ulp*bnorm
617 ELSE
618 atoli = abstol
619 END IF
620*
621 IF( gl.LT.initvl ) THEN
622 gl = initvl
623 IF( ieflag.EQ.0 ) THEN
624 CALL pdlapdct( gl, in, work( indrw1+2*ioff+1 ),
625 $ pivmin, ngl )
626 ELSE IF( ieflag.EQ.1 ) THEN
627 CALL pdlaiectb( gl, in, work( indrw1+2*ioff+1 ), ngl )
628 ELSE
629 CALL pdlaiectl( gl, in, work( indrw1+2*ioff+1 ), ngl )
630 END IF
631 ELSE
632 ngl = 0
633 END IF
634 IF( gu.GT.initvu ) THEN
635 gu = initvu
636 IF( ieflag.EQ.0 ) THEN
637 CALL pdlapdct( gu, in, work( indrw1+2*ioff+1 ),
638 $ pivmin, ngu )
639 ELSE IF( ieflag.EQ.1 ) THEN
640 CALL pdlaiectb( gu, in, work( indrw1+2*ioff+1 ), ngu )
641 ELSE
642 CALL pdlaiectl( gu, in, work( indrw1+2*ioff+1 ), ngu )
643 END IF
644 ELSE
645 ngu = in
646 END IF
647 IF( ngl.GE.ngu )
648 $ GO TO 90
649 work( 1 ) = gl
650 work( 2 ) = gu
651 iwork( 1 ) = ngl
652 iwork( 2 ) = ngu
653 END IF
654 offset = found - ngl
655 blkno = jb
656*
657* Do a static partitioning of work so that each process
658* has to find an (almost) equal number of eigenvalues
659*
660 ncmp = ngu - ngl
661 iload = ncmp / p
662 irem = ncmp - iload*p
663 itmp1 = mod( self-found, p )
664 IF( itmp1.LT.0 )
665 $ itmp1 = itmp1 + p
666 IF( itmp1.LT.irem ) THEN
667 imyload = iload + 1
668 ELSE
669 imyload = iload
670 END IF
671 IF( imyload.EQ.0 ) THEN
672 GO TO 90
673 ELSE IF( in.EQ.1 ) THEN
674 work( indrw2+im+1 ) = work( indrw1+2*ioff+1 )
675 iwork( indriw1+im+1 ) = blkno
676 iwork( indriw2+im+1 ) = offset + 1
677 im = im + 1
678 GO TO 90
679 ELSE
680 inxtload = iload
681 itmp2 = mod( self+1-found, p )
682 IF( itmp2.LT.0 )
683 $ itmp2 = itmp2 + p
684 IF( itmp2.LT.irem )
685 $ inxtload = inxtload + 1
686 lreq = ngl + itmp1*iload + min( irem, itmp1 )
687 rreq = lreq + imyload
688 iwork( 5 ) = lreq
689 iwork( 6 ) = rreq
690 tmp1 = work( 1 )
691 itmp1 = iwork( 1 )
692 CALL pdlaebz( 1, in, 1, 1, atoli, reltol, pivmin,
693 $ work( indrw1+2*ioff+1 ), iwork( 5 ), work,
694 $ iwork, nint, lsave, ieflag, iinfo )
695 alpha = work( 1 )
696 beta = work( 2 )
697 nalpha = iwork( 1 )
698 nbeta = iwork( 2 )
699 dsend = beta
700 IF( nbeta.GT.rreq+inxtload ) THEN
701 nbeta = rreq
702 dsend = alpha
703 END IF
704 last = mod( found+min( ngu-ngl, p )-1, p )
705 IF( last.LT.0 )
706 $ last = last + p
707 IF( self.NE.last ) THEN
708 CALL dgesd2d( onedcontext, 1, 1, dsend, 1, 0, next )
709 CALL igesd2d( onedcontext, 1, 1, nbeta, 1, 0, next )
710 END IF
711 IF( self.NE.mod( found, p ) ) THEN
712 CALL dgerv2d( onedcontext, 1, 1, drecv, 1, 0, prev )
713 CALL igerv2d( onedcontext, 1, 1, irecv, 1, 0, prev )
714 ELSE
715 drecv = tmp1
716 irecv = itmp1
717 END IF
718 work( 1 ) = max( lsave, drecv )
719 iwork( 1 ) = irecv
720 alpha = max( alpha, work( 1 ) )
721 nalpha = max( nalpha, irecv )
722 IF( beta-alpha.LE.max( atoli, reltol*max( abs( alpha ),
723 $ abs( beta ) ) ) ) THEN
724 mid = half*( alpha+beta )
725 DO 60 j = offset + nalpha + 1, offset + nbeta
726 work( indrw2+im+1 ) = mid
727 iwork( indriw1+im+1 ) = blkno
728 iwork( indriw2+im+1 ) = j
729 im = im + 1
730 60 CONTINUE
731 work( 2 ) = alpha
732 iwork( 2 ) = nalpha
733 END IF
734 END IF
735 neigint = iwork( 2 ) - iwork( 1 )
736 IF( neigint.LE.0 )
737 $ GO TO 90
738*
739* Call the main computational routine
740*
741 CALL pdlaebz( 2, in, neigint, 1, atoli, reltol, pivmin,
742 $ work( indrw1+2*ioff+1 ), iwork, work, iwork,
743 $ iout, lsave, ieflag, iinfo )
744 IF( iinfo.NE.0 ) THEN
745 info = 1
746 END IF
747 DO 80 i = 1, iout
748 mid = half*( work( 2*i-1 )+work( 2*i ) )
749 IF( i.GT.iout-iinfo )
750 $ blkno = -blkno
751 DO 70 j = offset + iwork( 2*i-1 ) + 1,
752 $ offset + iwork( 2*i )
753 work( indrw2+im+1 ) = mid
754 iwork( indriw1+im+1 ) = blkno
755 iwork( indriw2+im+1 ) = j
756 im = im + 1
757 70 CONTINUE
758 80 CONTINUE
759 90 CONTINUE
760*
761* Find out total number of eigenvalues computed
762*
763 100 CONTINUE
764 m = im
765 CALL igsum2d( onedcontext, 'ALL', ' ', 1, 1, m, 1, -1, -1 )
766*
767* Move the eigenvalues found to their final destinations
768*
769 DO 130 i = 1, p
770 IF( self.EQ.i-1 ) THEN
771 CALL igebs2d( onedcontext, 'ALL', ' ', 1, 1, im, 1 )
772 IF( im.NE.0 ) THEN
773 CALL igebs2d( onedcontext, 'ALL', ' ', im, 1,
774 $ iwork( indriw2+1 ), im )
775 CALL dgebs2d( onedcontext, 'ALL', ' ', im, 1,
776 $ work( indrw2+1 ), im )
777 CALL igebs2d( onedcontext, 'ALL', ' ', im, 1,
778 $ iwork( indriw1+1 ), im )
779 DO 110 j = 1, im
780 w( iwork( indriw2+j ) ) = work( indrw2+j )
781 iblock( iwork( indriw2+j ) ) = iwork( indriw1+j )
782 110 CONTINUE
783 END IF
784 ELSE
785 CALL igebr2d( onedcontext, 'ALL', ' ', 1, 1, torecv, 1, 0,
786 $ i-1 )
787 IF( torecv( 1, 1 ).NE.0 ) THEN
788 CALL igebr2d( onedcontext, 'all', ' ', TORECV( 1, 1 ), 1,
789 $ IWORK, TORECV( 1, 1 ), 0, I-1 )
790 CALL DGEBR2D( ONEDCONTEXT, 'all', ' ', TORECV( 1, 1 ), 1,
791 $ WORK, TORECV( 1, 1 ), 0, I-1 )
792 CALL IGEBR2D( ONEDCONTEXT, 'all', ' ', TORECV( 1, 1 ), 1,
793 $ IWORK( N+1 ), TORECV( 1, 1 ), 0, I-1 )
794 DO 120 J = 1, TORECV( 1, 1 )
795 W( IWORK( J ) ) = WORK( J )
796 IBLOCK( IWORK( J ) ) = IWORK( N+J )
797 120 CONTINUE
798 END IF
799 END IF
800 130 CONTINUE
801.GT..AND..EQ. IF( NSPLIT1 IORDER1 ) THEN
802*
803* Sort the eigenvalues
804*
805*
806 DO 140 I = 1, M
807 IWORK( M+I ) = I
808 140 CONTINUE
809 CALL DLASRT2( 'i', M, W, IWORK( M+1 ), IINFO )
810 DO 150 I = 1, M
811 IWORK( I ) = IBLOCK( I )
812 150 CONTINUE
813 DO 160 I = 1, M
814 IBLOCK( I ) = IWORK( IWORK( M+I ) )
815 160 CONTINUE
816 END IF
817.EQ..AND..GT..OR..GT. IF( IRANGE3 ( LEXTRA0 REXTRA0 ) ) THEN
818*
819* Discard unwanted eigenvalues (occurs only when RANGE = 'I',
820* and eigenvalues IL, and/or IU are in a cluster)
821*
822 DO 170 I = 1, M
823 WORK( I ) = W( I )
824 IWORK( I ) = I
825 IWORK( M+I ) = I
826 170 CONTINUE
827 DO 190 I = 1, LEXTRA
828 ITMP1 = I
829 DO 180 J = I + 1, M
830.LT. IF( WORK( J )WORK( ITMP1 ) ) THEN
831 ITMP1 = J
832 END IF
833 180 CONTINUE
834 TMP1 = WORK( I )
835 WORK( I ) = WORK( ITMP1 )
836 WORK( ITMP1 ) = TMP1
837 IWORK( IWORK( M+ITMP1 ) ) = I
838 IWORK( IWORK( M+I ) ) = ITMP1
839 ITMP2 = IWORK( M+I )
840 IWORK( M+I ) = IWORK( M+ITMP1 )
841 IWORK( M+ITMP1 ) = ITMP2
842 190 CONTINUE
843 DO 210 I = 1, REXTRA
844 ITMP1 = M - I + 1
845 DO 200 J = M - I, LEXTRA + 1, -1
846.GT. IF( WORK( J )WORK( ITMP1 ) ) THEN
847 ITMP1 = J
848 END IF
849 200 CONTINUE
850 TMP1 = WORK( M-I+1 )
851 WORK( M-I+1 ) = WORK( ITMP1 )
852 WORK( ITMP1 ) = TMP1
853 IWORK( IWORK( M+ITMP1 ) ) = M - I + 1
854 IWORK( IWORK( 2*M-I+1 ) ) = ITMP1
855 ITMP2 = IWORK( 2*M-I+1 )
856 IWORK( 2*M-I+1 ) = IWORK( M+ITMP1 )
857 IWORK( M+ITMP1 ) = ITMP2
858* IWORK( ITMP1 ) = 1
859 210 CONTINUE
860 J = 0
861 DO 220 I = 1, M
862.GT..AND..LE. IF( IWORK( I )LEXTRA IWORK( I )M-REXTRA ) THEN
863 J = J + 1
864 W( J ) = WORK( IWORK( I ) )
865 IBLOCK( J ) = IBLOCK( I )
866 END IF
867 220 CONTINUE
868 M = M - LEXTRA - REXTRA
869 END IF
870.NE. IF( MILAST-IFRST+1 ) THEN
871 INFO = 2
872 END IF
873*
874 230 CONTINUE
875 CALL BLACS_FREEBUFF( ONEDCONTEXT, 1 )
876 CALL BLACS_GRIDEXIT( ONEDCONTEXT )
877 RETURN
878*
879* End of PDSTEBZ
880*
dlasrt2
subroutine dlasrt2(id, n, d, key, info)
Definition dlasrt2.f:4
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:53
dgebs2d
subroutine dgebs2d(contxt, scope, top, m, n, a, lda)
Definition mpi.f:1082
pxerbla
subroutine pxerbla(contxt, srname, info)
Definition mpi.f:1600
dgebr2d
subroutine dgebr2d(contxt, scope, top, m, n, a, lda)
Definition mpi.f:1123
blacs_gridexit
subroutine blacs_gridexit(cntxt)
Definition mpi.f:762
blacs_gridinfo
subroutine blacs_gridinfo(cntxt, nprow, npcol, myrow, mycol)
Definition mpi.f:754
globchk
subroutine globchk(ictxt, n, x, ldx, iwork, info)
Definition pchkxmat.f:403
pdlamch
double precision function pdlamch(ictxt, cmach)
Definition pdblastst.f:6769
pdlaebz
subroutine pdlaebz(ijob, n, mmax, minp, abstol, reltol, pivmin, d, nval, intvl, intvlct, mout, lsave, ieflag, info)
Definition pdstebz.f:886