![]() ![]() | ![]() |
Mise à jour de lapack vers la version 3.2.2.
1: SUBROUTINE DLARRJ( N, D, E2, IFIRST, ILAST, 2: $ RTOL, OFFSET, W, WERR, WORK, IWORK, 3: $ PIVMIN, SPDIAM, INFO ) 4: * 5: * -- LAPACK auxiliary routine (version 3.2.2) -- 6: * -- LAPACK is a software package provided by Univ. of Tennessee, -- 7: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 8: * June 2010 9: * 10: * .. Scalar Arguments .. 11: INTEGER IFIRST, ILAST, INFO, N, OFFSET 12: DOUBLE PRECISION PIVMIN, RTOL, SPDIAM 13: * .. 14: * .. Array Arguments .. 15: INTEGER IWORK( * ) 16: DOUBLE PRECISION D( * ), E2( * ), W( * ), 17: $ WERR( * ), WORK( * ) 18: * .. 19: * 20: * Purpose 21: * ======= 22: * 23: * Given the initial eigenvalue approximations of T, DLARRJ 24: * does bisection to refine the eigenvalues of T, 25: * W( IFIRST-OFFSET ) through W( ILAST-OFFSET ), to more accuracy. Initial 26: * guesses for these eigenvalues are input in W, the corresponding estimate 27: * of the error in these guesses in WERR. During bisection, intervals 28: * [left, right] are maintained by storing their mid-points and 29: * semi-widths in the arrays W and WERR respectively. 30: * 31: * Arguments 32: * ========= 33: * 34: * N (input) INTEGER 35: * The order of the matrix. 36: * 37: * D (input) DOUBLE PRECISION array, dimension (N) 38: * The N diagonal elements of T. 39: * 40: * E2 (input) DOUBLE PRECISION array, dimension (N-1) 41: * The Squares of the (N-1) subdiagonal elements of T. 42: * 43: * IFIRST (input) INTEGER 44: * The index of the first eigenvalue to be computed. 45: * 46: * ILAST (input) INTEGER 47: * The index of the last eigenvalue to be computed. 48: * 49: * RTOL (input) DOUBLE PRECISION 50: * Tolerance for the convergence of the bisection intervals. 51: * An interval [LEFT,RIGHT] has converged if 52: * RIGHT-LEFT.LT.RTOL*MAX(|LEFT|,|RIGHT|). 53: * 54: * OFFSET (input) INTEGER 55: * Offset for the arrays W and WERR, i.e., the IFIRST-OFFSET 56: * through ILAST-OFFSET elements of these arrays are to be used. 57: * 58: * W (input/output) DOUBLE PRECISION array, dimension (N) 59: * On input, W( IFIRST-OFFSET ) through W( ILAST-OFFSET ) are 60: * estimates of the eigenvalues of L D L^T indexed IFIRST through 61: * ILAST. 62: * On output, these estimates are refined. 63: * 64: * WERR (input/output) DOUBLE PRECISION array, dimension (N) 65: * On input, WERR( IFIRST-OFFSET ) through WERR( ILAST-OFFSET ) are 66: * the errors in the estimates of the corresponding elements in W. 67: * On output, these errors are refined. 68: * 69: * WORK (workspace) DOUBLE PRECISION array, dimension (2*N) 70: * Workspace. 71: * 72: * IWORK (workspace) INTEGER array, dimension (2*N) 73: * Workspace. 74: * 75: * PIVMIN (input) DOUBLE PRECISION 76: * The minimum pivot in the Sturm sequence for T. 77: * 78: * SPDIAM (input) DOUBLE PRECISION 79: * The spectral diameter of T. 80: * 81: * INFO (output) INTEGER 82: * Error flag. 83: * 84: * Further Details 85: * =============== 86: * 87: * Based on contributions by 88: * Beresford Parlett, University of California, Berkeley, USA 89: * Jim Demmel, University of California, Berkeley, USA 90: * Inderjit Dhillon, University of Texas, Austin, USA 91: * Osni Marques, LBNL/NERSC, USA 92: * Christof Voemel, University of California, Berkeley, USA 93: * 94: * ===================================================================== 95: * 96: * .. Parameters .. 97: DOUBLE PRECISION ZERO, ONE, TWO, HALF 98: PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0, 99: $ HALF = 0.5D0 ) 100: INTEGER MAXITR 101: * .. 102: * .. Local Scalars .. 103: INTEGER CNT, I, I1, I2, II, ITER, J, K, NEXT, NINT, 104: $ OLNINT, P, PREV, SAVI1 105: DOUBLE PRECISION DPLUS, FAC, LEFT, MID, RIGHT, S, TMP, WIDTH 106: * 107: * .. 108: * .. Intrinsic Functions .. 109: INTRINSIC ABS, MAX 110: * .. 111: * .. Executable Statements .. 112: * 113: INFO = 0 114: * 115: MAXITR = INT( ( LOG( SPDIAM+PIVMIN )-LOG( PIVMIN ) ) / 116: $ LOG( TWO ) ) + 2 117: * 118: * Initialize unconverged intervals in [ WORK(2*I-1), WORK(2*I) ]. 119: * The Sturm Count, Count( WORK(2*I-1) ) is arranged to be I-1, while 120: * Count( WORK(2*I) ) is stored in IWORK( 2*I ). The integer IWORK( 2*I-1 ) 121: * for an unconverged interval is set to the index of the next unconverged 122: * interval, and is -1 or 0 for a converged interval. Thus a linked 123: * list of unconverged intervals is set up. 124: * 125: 126: I1 = IFIRST 127: I2 = ILAST 128: * The number of unconverged intervals 129: NINT = 0 130: * The last unconverged interval found 131: PREV = 0 132: DO 75 I = I1, I2 133: K = 2*I 134: II = I - OFFSET 135: LEFT = W( II ) - WERR( II ) 136: MID = W(II) 137: RIGHT = W( II ) + WERR( II ) 138: WIDTH = RIGHT - MID 139: TMP = MAX( ABS( LEFT ), ABS( RIGHT ) ) 140: 141: * The following test prevents the test of converged intervals 142: IF( WIDTH.LT.RTOL*TMP ) THEN 143: * This interval has already converged and does not need refinement. 144: * (Note that the gaps might change through refining the 145: * eigenvalues, however, they can only get bigger.) 146: * Remove it from the list. 147: IWORK( K-1 ) = -1 148: * Make sure that I1 always points to the first unconverged interval 149: IF((I.EQ.I1).AND.(I.LT.I2)) I1 = I + 1 150: IF((PREV.GE.I1).AND.(I.LE.I2)) IWORK( 2*PREV-1 ) = I + 1 151: ELSE 152: * unconverged interval found 153: PREV = I 154: * Make sure that [LEFT,RIGHT] contains the desired eigenvalue 155: * 156: * Do while( CNT(LEFT).GT.I-1 ) 157: * 158: FAC = ONE 159: 20 CONTINUE 160: CNT = 0 161: S = LEFT 162: DPLUS = D( 1 ) - S 163: IF( DPLUS.LT.ZERO ) CNT = CNT + 1 164: DO 30 J = 2, N 165: DPLUS = D( J ) - S - E2( J-1 )/DPLUS 166: IF( DPLUS.LT.ZERO ) CNT = CNT + 1 167: 30 CONTINUE 168: IF( CNT.GT.I-1 ) THEN 169: LEFT = LEFT - WERR( II )*FAC 170: FAC = TWO*FAC 171: GO TO 20 172: END IF 173: * 174: * Do while( CNT(RIGHT).LT.I ) 175: * 176: FAC = ONE 177: 50 CONTINUE 178: CNT = 0 179: S = RIGHT 180: DPLUS = D( 1 ) - S 181: IF( DPLUS.LT.ZERO ) CNT = CNT + 1 182: DO 60 J = 2, N 183: DPLUS = D( J ) - S - E2( J-1 )/DPLUS 184: IF( DPLUS.LT.ZERO ) CNT = CNT + 1 185: 60 CONTINUE 186: IF( CNT.LT.I ) THEN 187: RIGHT = RIGHT + WERR( II )*FAC 188: FAC = TWO*FAC 189: GO TO 50 190: END IF 191: NINT = NINT + 1 192: IWORK( K-1 ) = I + 1 193: IWORK( K ) = CNT 194: END IF 195: WORK( K-1 ) = LEFT 196: WORK( K ) = RIGHT 197: 75 CONTINUE 198: 199: 200: SAVI1 = I1 201: * 202: * Do while( NINT.GT.0 ), i.e. there are still unconverged intervals 203: * and while (ITER.LT.MAXITR) 204: * 205: ITER = 0 206: 80 CONTINUE 207: PREV = I1 - 1 208: I = I1 209: OLNINT = NINT 210: 211: DO 100 P = 1, OLNINT 212: K = 2*I 213: II = I - OFFSET 214: NEXT = IWORK( K-1 ) 215: LEFT = WORK( K-1 ) 216: RIGHT = WORK( K ) 217: MID = HALF*( LEFT + RIGHT ) 218: 219: * semiwidth of interval 220: WIDTH = RIGHT - MID 221: TMP = MAX( ABS( LEFT ), ABS( RIGHT ) ) 222: 223: IF( ( WIDTH.LT.RTOL*TMP ) .OR. 224: $ (ITER.EQ.MAXITR) )THEN 225: * reduce number of unconverged intervals 226: NINT = NINT - 1 227: * Mark interval as converged. 228: IWORK( K-1 ) = 0 229: IF( I1.EQ.I ) THEN 230: I1 = NEXT 231: ELSE 232: * Prev holds the last unconverged interval previously examined 233: IF(PREV.GE.I1) IWORK( 2*PREV-1 ) = NEXT 234: END IF 235: I = NEXT 236: GO TO 100 237: END IF 238: PREV = I 239: * 240: * Perform one bisection step 241: * 242: CNT = 0 243: S = MID 244: DPLUS = D( 1 ) - S 245: IF( DPLUS.LT.ZERO ) CNT = CNT + 1 246: DO 90 J = 2, N 247: DPLUS = D( J ) - S - E2( J-1 )/DPLUS 248: IF( DPLUS.LT.ZERO ) CNT = CNT + 1 249: 90 CONTINUE 250: IF( CNT.LE.I-1 ) THEN 251: WORK( K-1 ) = MID 252: ELSE 253: WORK( K ) = MID 254: END IF 255: I = NEXT 256: 257: 100 CONTINUE 258: ITER = ITER + 1 259: * do another loop if there are still unconverged intervals 260: * However, in the last iteration, all intervals are accepted 261: * since this is the best we can do. 262: IF( ( NINT.GT.0 ).AND.(ITER.LE.MAXITR) ) GO TO 80 263: * 264: * 265: * At this point, all the intervals have converged 266: DO 110 I = SAVI1, ILAST 267: K = 2*I 268: II = I - OFFSET 269: * All intervals marked by '0' have been refined. 270: IF( IWORK( K-1 ).EQ.0 ) THEN 271: W( II ) = HALF*( WORK( K-1 )+WORK( K ) ) 272: WERR( II ) = WORK( K ) - W( II ) 273: END IF 274: 110 CONTINUE 275: * 276: 277: RETURN 278: * 279: * End of DLARRJ 280: * 281: END