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

Go to the source code of this file.

Functions/Subroutines

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

Function/Subroutine Documentation

◆ pslaebz()

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

Definition at line 873 of file psstebz.f.

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

◆ pslaecv()

subroutine pslaecv ( integer ijob,
integer kf,
integer kl,
real, dimension( * ) intvl,
integer, dimension( * ) intvlct,
integer, dimension( * ) nval,
real abstol,
real reltol )

Definition at line 1199 of file psstebz.f.

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

◆ pslapdct()

subroutine pslapdct ( real sigma,
integer n,
real, dimension( * ) d,
real pivmin,
integer count )

Definition at line 1348 of file psstebz.f.

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

◆ psstebz()

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

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