![]() ![]() | ![]() |
Mise à jour de lapack vers la version 3.3.0.
1: SUBROUTINE DLAEBZ( IJOB, NITMAX, N, MMAX, MINP, NBMIN, ABSTOL, 2: $ RELTOL, PIVMIN, D, E, E2, NVAL, AB, C, MOUT, 3: $ NAB, WORK, IWORK, INFO ) 4: * 5: * -- LAPACK auxiliary routine (version 3.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: * November 2006 9: * 10: * .. Scalar Arguments .. 11: INTEGER IJOB, INFO, MINP, MMAX, MOUT, N, NBMIN, NITMAX 12: DOUBLE PRECISION ABSTOL, PIVMIN, RELTOL 13: * .. 14: * .. Array Arguments .. 15: INTEGER IWORK( * ), NAB( MMAX, * ), NVAL( * ) 16: DOUBLE PRECISION AB( MMAX, * ), C( * ), D( * ), E( * ), E2( * ), 17: $ WORK( * ) 18: * .. 19: * 20: * Purpose 21: * ======= 22: * 23: * DLAEBZ contains the iteration loops which compute and use the 24: * function N(w), which is the count of eigenvalues of a symmetric 25: * tridiagonal matrix T less than or equal to its argument w. It 26: * performs a choice of two types of loops: 27: * 28: * IJOB=1, followed by 29: * IJOB=2: It takes as input a list of intervals and returns a list of 30: * sufficiently small intervals whose union contains the same 31: * eigenvalues as the union of the original intervals. 32: * The input intervals are (AB(j,1),AB(j,2)], j=1,...,MINP. 33: * The output interval (AB(j,1),AB(j,2)] will contain 34: * eigenvalues NAB(j,1)+1,...,NAB(j,2), where 1 <= j <= MOUT. 35: * 36: * IJOB=3: It performs a binary search in each input interval 37: * (AB(j,1),AB(j,2)] for a point w(j) such that 38: * N(w(j))=NVAL(j), and uses C(j) as the starting point of 39: * the search. If such a w(j) is found, then on output 40: * AB(j,1)=AB(j,2)=w. If no such w(j) is found, then on output 41: * (AB(j,1),AB(j,2)] will be a small interval containing the 42: * point where N(w) jumps through NVAL(j), unless that point 43: * lies outside the initial interval. 44: * 45: * Note that the intervals are in all cases half-open intervals, 46: * i.e., of the form (a,b] , which includes b but not a . 47: * 48: * To avoid underflow, the matrix should be scaled so that its largest 49: * element is no greater than overflow**(1/2) * underflow**(1/4) 50: * in absolute value. To assure the most accurate computation 51: * of small eigenvalues, the matrix should be scaled to be 52: * not much smaller than that, either. 53: * 54: * See W. Kahan "Accurate Eigenvalues of a Symmetric Tridiagonal 55: * Matrix", Report CS41, Computer Science Dept., Stanford 56: * University, July 21, 1966 57: * 58: * Note: the arguments are, in general, *not* checked for unreasonable 59: * values. 60: * 61: * Arguments 62: * ========= 63: * 64: * IJOB (input) INTEGER 65: * Specifies what is to be done: 66: * = 1: Compute NAB for the initial intervals. 67: * = 2: Perform bisection iteration to find eigenvalues of T. 68: * = 3: Perform bisection iteration to invert N(w), i.e., 69: * to find a point which has a specified number of 70: * eigenvalues of T to its left. 71: * Other values will cause DLAEBZ to return with INFO=-1. 72: * 73: * NITMAX (input) INTEGER 74: * The maximum number of "levels" of bisection to be 75: * performed, i.e., an interval of width W will not be made 76: * smaller than 2^(-NITMAX) * W. If not all intervals 77: * have converged after NITMAX iterations, then INFO is set 78: * to the number of non-converged intervals. 79: * 80: * N (input) INTEGER 81: * The dimension n of the tridiagonal matrix T. It must be at 82: * least 1. 83: * 84: * MMAX (input) INTEGER 85: * The maximum number of intervals. If more than MMAX intervals 86: * are generated, then DLAEBZ will quit with INFO=MMAX+1. 87: * 88: * MINP (input) INTEGER 89: * The initial number of intervals. It may not be greater than 90: * MMAX. 91: * 92: * NBMIN (input) INTEGER 93: * The smallest number of intervals that should be processed 94: * using a vector loop. If zero, then only the scalar loop 95: * will be used. 96: * 97: * ABSTOL (input) DOUBLE PRECISION 98: * The minimum (absolute) width of an interval. When an 99: * interval is narrower than ABSTOL, or than RELTOL times the 100: * larger (in magnitude) endpoint, then it is considered to be 101: * sufficiently small, i.e., converged. This must be at least 102: * zero. 103: * 104: * RELTOL (input) DOUBLE PRECISION 105: * The minimum relative width of an interval. When an interval 106: * is narrower than ABSTOL, or than RELTOL times the larger (in 107: * magnitude) endpoint, then it is considered to be 108: * sufficiently small, i.e., converged. Note: this should 109: * always be at least radix*machine epsilon. 110: * 111: * PIVMIN (input) DOUBLE PRECISION 112: * The minimum absolute value of a "pivot" in the Sturm 113: * sequence loop. This *must* be at least max |e(j)**2| * 114: * safe_min and at least safe_min, where safe_min is at least 115: * the smallest number that can divide one without overflow. 116: * 117: * D (input) DOUBLE PRECISION array, dimension (N) 118: * The diagonal elements of the tridiagonal matrix T. 119: * 120: * E (input) DOUBLE PRECISION array, dimension (N) 121: * The offdiagonal elements of the tridiagonal matrix T in 122: * positions 1 through N-1. E(N) is arbitrary. 123: * 124: * E2 (input) DOUBLE PRECISION array, dimension (N) 125: * The squares of the offdiagonal elements of the tridiagonal 126: * matrix T. E2(N) is ignored. 127: * 128: * NVAL (input/output) INTEGER array, dimension (MINP) 129: * If IJOB=1 or 2, not referenced. 130: * If IJOB=3, the desired values of N(w). The elements of NVAL 131: * will be reordered to correspond with the intervals in AB. 132: * Thus, NVAL(j) on output will not, in general be the same as 133: * NVAL(j) on input, but it will correspond with the interval 134: * (AB(j,1),AB(j,2)] on output. 135: * 136: * AB (input/output) DOUBLE PRECISION array, dimension (MMAX,2) 137: * The endpoints of the intervals. AB(j,1) is a(j), the left 138: * endpoint of the j-th interval, and AB(j,2) is b(j), the 139: * right endpoint of the j-th interval. The input intervals 140: * will, in general, be modified, split, and reordered by the 141: * calculation. 142: * 143: * C (input/output) DOUBLE PRECISION array, dimension (MMAX) 144: * If IJOB=1, ignored. 145: * If IJOB=2, workspace. 146: * If IJOB=3, then on input C(j) should be initialized to the 147: * first search point in the binary search. 148: * 149: * MOUT (output) INTEGER 150: * If IJOB=1, the number of eigenvalues in the intervals. 151: * If IJOB=2 or 3, the number of intervals output. 152: * If IJOB=3, MOUT will equal MINP. 153: * 154: * NAB (input/output) INTEGER array, dimension (MMAX,2) 155: * If IJOB=1, then on output NAB(i,j) will be set to N(AB(i,j)). 156: * If IJOB=2, then on input, NAB(i,j) should be set. It must 157: * satisfy the condition: 158: * N(AB(i,1)) <= NAB(i,1) <= NAB(i,2) <= N(AB(i,2)), 159: * which means that in interval i only eigenvalues 160: * NAB(i,1)+1,...,NAB(i,2) will be considered. Usually, 161: * NAB(i,j)=N(AB(i,j)), from a previous call to DLAEBZ with 162: * IJOB=1. 163: * On output, NAB(i,j) will contain 164: * max(na(k),min(nb(k),N(AB(i,j)))), where k is the index of 165: * the input interval that the output interval 166: * (AB(j,1),AB(j,2)] came from, and na(k) and nb(k) are the 167: * the input values of NAB(k,1) and NAB(k,2). 168: * If IJOB=3, then on output, NAB(i,j) contains N(AB(i,j)), 169: * unless N(w) > NVAL(i) for all search points w , in which 170: * case NAB(i,1) will not be modified, i.e., the output 171: * value will be the same as the input value (modulo 172: * reorderings -- see NVAL and AB), or unless N(w) < NVAL(i) 173: * for all search points w , in which case NAB(i,2) will 174: * not be modified. Normally, NAB should be set to some 175: * distinctive value(s) before DLAEBZ is called. 176: * 177: * WORK (workspace) DOUBLE PRECISION array, dimension (MMAX) 178: * Workspace. 179: * 180: * IWORK (workspace) INTEGER array, dimension (MMAX) 181: * Workspace. 182: * 183: * INFO (output) INTEGER 184: * = 0: All intervals converged. 185: * = 1--MMAX: The last INFO intervals did not converge. 186: * = MMAX+1: More than MMAX intervals were generated. 187: * 188: * Further Details 189: * =============== 190: * 191: * This routine is intended to be called only by other LAPACK 192: * routines, thus the interface is less user-friendly. It is intended 193: * for two purposes: 194: * 195: * (a) finding eigenvalues. In this case, DLAEBZ should have one or 196: * more initial intervals set up in AB, and DLAEBZ should be called 197: * with IJOB=1. This sets up NAB, and also counts the eigenvalues. 198: * Intervals with no eigenvalues would usually be thrown out at 199: * this point. Also, if not all the eigenvalues in an interval i 200: * are desired, NAB(i,1) can be increased or NAB(i,2) decreased. 201: * For example, set NAB(i,1)=NAB(i,2)-1 to get the largest 202: * eigenvalue. DLAEBZ is then called with IJOB=2 and MMAX 203: * no smaller than the value of MOUT returned by the call with 204: * IJOB=1. After this (IJOB=2) call, eigenvalues NAB(i,1)+1 205: * through NAB(i,2) are approximately AB(i,1) (or AB(i,2)) to the 206: * tolerance specified by ABSTOL and RELTOL. 207: * 208: * (b) finding an interval (a',b'] containing eigenvalues w(f),...,w(l). 209: * In this case, start with a Gershgorin interval (a,b). Set up 210: * AB to contain 2 search intervals, both initially (a,b). One 211: * NVAL element should contain f-1 and the other should contain l 212: * , while C should contain a and b, resp. NAB(i,1) should be -1 213: * and NAB(i,2) should be N+1, to flag an error if the desired 214: * interval does not lie in (a,b). DLAEBZ is then called with 215: * IJOB=3. On exit, if w(f-1) < w(f), then one of the intervals -- 216: * j -- will have AB(j,1)=AB(j,2) and NAB(j,1)=NAB(j,2)=f-1, while 217: * if, to the specified tolerance, w(f-k)=...=w(f+r), k > 0 and r 218: * >= 0, then the interval will have N(AB(j,1))=NAB(j,1)=f-k and 219: * N(AB(j,2))=NAB(j,2)=f+r. The cases w(l) < w(l+1) and 220: * w(l-r)=...=w(l+k) are handled similarly. 221: * 222: * ===================================================================== 223: * 224: * .. Parameters .. 225: DOUBLE PRECISION ZERO, TWO, HALF 226: PARAMETER ( ZERO = 0.0D0, TWO = 2.0D0, 227: $ HALF = 1.0D0 / TWO ) 228: * .. 229: * .. Local Scalars .. 230: INTEGER ITMP1, ITMP2, J, JI, JIT, JP, KF, KFNEW, KL, 231: $ KLNEW 232: DOUBLE PRECISION TMP1, TMP2 233: * .. 234: * .. Intrinsic Functions .. 235: INTRINSIC ABS, MAX, MIN 236: * .. 237: * .. Executable Statements .. 238: * 239: * Check for Errors 240: * 241: INFO = 0 242: IF( IJOB.LT.1 .OR. IJOB.GT.3 ) THEN 243: INFO = -1 244: RETURN 245: END IF 246: * 247: * Initialize NAB 248: * 249: IF( IJOB.EQ.1 ) THEN 250: * 251: * Compute the number of eigenvalues in the initial intervals. 252: * 253: MOUT = 0 254: *DIR$ NOVECTOR 255: DO 30 JI = 1, MINP 256: DO 20 JP = 1, 2 257: TMP1 = D( 1 ) - AB( JI, JP ) 258: IF( ABS( TMP1 ).LT.PIVMIN ) 259: $ TMP1 = -PIVMIN 260: NAB( JI, JP ) = 0 261: IF( TMP1.LE.ZERO ) 262: $ NAB( JI, JP ) = 1 263: * 264: DO 10 J = 2, N 265: TMP1 = D( J ) - E2( J-1 ) / TMP1 - AB( JI, JP ) 266: IF( ABS( TMP1 ).LT.PIVMIN ) 267: $ TMP1 = -PIVMIN 268: IF( TMP1.LE.ZERO ) 269: $ NAB( JI, JP ) = NAB( JI, JP ) + 1 270: 10 CONTINUE 271: 20 CONTINUE 272: MOUT = MOUT + NAB( JI, 2 ) - NAB( JI, 1 ) 273: 30 CONTINUE 274: RETURN 275: END IF 276: * 277: * Initialize for loop 278: * 279: * KF and KL have the following meaning: 280: * Intervals 1,...,KF-1 have converged. 281: * Intervals KF,...,KL still need to be refined. 282: * 283: KF = 1 284: KL = MINP 285: * 286: * If IJOB=2, initialize C. 287: * If IJOB=3, use the user-supplied starting point. 288: * 289: IF( IJOB.EQ.2 ) THEN 290: DO 40 JI = 1, MINP 291: C( JI ) = HALF*( AB( JI, 1 )+AB( JI, 2 ) ) 292: 40 CONTINUE 293: END IF 294: * 295: * Iteration loop 296: * 297: DO 130 JIT = 1, NITMAX 298: * 299: * Loop over intervals 300: * 301: IF( KL-KF+1.GE.NBMIN .AND. NBMIN.GT.0 ) THEN 302: * 303: * Begin of Parallel Version of the loop 304: * 305: DO 60 JI = KF, KL 306: * 307: * Compute N(c), the number of eigenvalues less than c 308: * 309: WORK( JI ) = D( 1 ) - C( JI ) 310: IWORK( JI ) = 0 311: IF( WORK( JI ).LE.PIVMIN ) THEN 312: IWORK( JI ) = 1 313: WORK( JI ) = MIN( WORK( JI ), -PIVMIN ) 314: END IF 315: * 316: DO 50 J = 2, N 317: WORK( JI ) = D( J ) - E2( J-1 ) / WORK( JI ) - C( JI ) 318: IF( WORK( JI ).LE.PIVMIN ) THEN 319: IWORK( JI ) = IWORK( JI ) + 1 320: WORK( JI ) = MIN( WORK( JI ), -PIVMIN ) 321: END IF 322: 50 CONTINUE 323: 60 CONTINUE 324: * 325: IF( IJOB.LE.2 ) THEN 326: * 327: * IJOB=2: Choose all intervals containing eigenvalues. 328: * 329: KLNEW = KL 330: DO 70 JI = KF, KL 331: * 332: * Insure that N(w) is monotone 333: * 334: IWORK( JI ) = MIN( NAB( JI, 2 ), 335: $ MAX( NAB( JI, 1 ), IWORK( JI ) ) ) 336: * 337: * Update the Queue -- add intervals if both halves 338: * contain eigenvalues. 339: * 340: IF( IWORK( JI ).EQ.NAB( JI, 2 ) ) THEN 341: * 342: * No eigenvalue in the upper interval: 343: * just use the lower interval. 344: * 345: AB( JI, 2 ) = C( JI ) 346: * 347: ELSE IF( IWORK( JI ).EQ.NAB( JI, 1 ) ) THEN 348: * 349: * No eigenvalue in the lower interval: 350: * just use the upper interval. 351: * 352: AB( JI, 1 ) = C( JI ) 353: ELSE 354: KLNEW = KLNEW + 1 355: IF( KLNEW.LE.MMAX ) THEN 356: * 357: * Eigenvalue in both intervals -- add upper to 358: * queue. 359: * 360: AB( KLNEW, 2 ) = AB( JI, 2 ) 361: NAB( KLNEW, 2 ) = NAB( JI, 2 ) 362: AB( KLNEW, 1 ) = C( JI ) 363: NAB( KLNEW, 1 ) = IWORK( JI ) 364: AB( JI, 2 ) = C( JI ) 365: NAB( JI, 2 ) = IWORK( JI ) 366: ELSE 367: INFO = MMAX + 1 368: END IF 369: END IF 370: 70 CONTINUE 371: IF( INFO.NE.0 ) 372: $ RETURN 373: KL = KLNEW 374: ELSE 375: * 376: * IJOB=3: Binary search. Keep only the interval containing 377: * w s.t. N(w) = NVAL 378: * 379: DO 80 JI = KF, KL 380: IF( IWORK( JI ).LE.NVAL( JI ) ) THEN 381: AB( JI, 1 ) = C( JI ) 382: NAB( JI, 1 ) = IWORK( JI ) 383: END IF 384: IF( IWORK( JI ).GE.NVAL( JI ) ) THEN 385: AB( JI, 2 ) = C( JI ) 386: NAB( JI, 2 ) = IWORK( JI ) 387: END IF 388: 80 CONTINUE 389: END IF 390: * 391: ELSE 392: * 393: * End of Parallel Version of the loop 394: * 395: * Begin of Serial Version of the loop 396: * 397: KLNEW = KL 398: DO 100 JI = KF, KL 399: * 400: * Compute N(w), the number of eigenvalues less than w 401: * 402: TMP1 = C( JI ) 403: TMP2 = D( 1 ) - TMP1 404: ITMP1 = 0 405: IF( TMP2.LE.PIVMIN ) THEN 406: ITMP1 = 1 407: TMP2 = MIN( TMP2, -PIVMIN ) 408: END IF 409: * 410: * A series of compiler directives to defeat vectorization 411: * for the next loop 412: * 413: *$PL$ CMCHAR=' ' 414: CDIR$ NEXTSCALAR 415: C$DIR SCALAR 416: CDIR$ NEXT SCALAR 417: CVD$L NOVECTOR 418: CDEC$ NOVECTOR 419: CVD$ NOVECTOR 420: *VDIR NOVECTOR 421: *VOCL LOOP,SCALAR 422: CIBM PREFER SCALAR 423: *$PL$ CMCHAR='*' 424: * 425: DO 90 J = 2, N 426: TMP2 = D( J ) - E2( J-1 ) / TMP2 - TMP1 427: IF( TMP2.LE.PIVMIN ) THEN 428: ITMP1 = ITMP1 + 1 429: TMP2 = MIN( TMP2, -PIVMIN ) 430: END IF 431: 90 CONTINUE 432: * 433: IF( IJOB.LE.2 ) THEN 434: * 435: * IJOB=2: Choose all intervals containing eigenvalues. 436: * 437: * Insure that N(w) is monotone 438: * 439: ITMP1 = MIN( NAB( JI, 2 ), 440: $ MAX( NAB( JI, 1 ), ITMP1 ) ) 441: * 442: * Update the Queue -- add intervals if both halves 443: * contain eigenvalues. 444: * 445: IF( ITMP1.EQ.NAB( JI, 2 ) ) THEN 446: * 447: * No eigenvalue in the upper interval: 448: * just use the lower interval. 449: * 450: AB( JI, 2 ) = TMP1 451: * 452: ELSE IF( ITMP1.EQ.NAB( JI, 1 ) ) THEN 453: * 454: * No eigenvalue in the lower interval: 455: * just use the upper interval. 456: * 457: AB( JI, 1 ) = TMP1 458: ELSE IF( KLNEW.LT.MMAX ) THEN 459: * 460: * Eigenvalue in both intervals -- add upper to queue. 461: * 462: KLNEW = KLNEW + 1 463: AB( KLNEW, 2 ) = AB( JI, 2 ) 464: NAB( KLNEW, 2 ) = NAB( JI, 2 ) 465: AB( KLNEW, 1 ) = TMP1 466: NAB( KLNEW, 1 ) = ITMP1 467: AB( JI, 2 ) = TMP1 468: NAB( JI, 2 ) = ITMP1 469: ELSE 470: INFO = MMAX + 1 471: RETURN 472: END IF 473: ELSE 474: * 475: * IJOB=3: Binary search. Keep only the interval 476: * containing w s.t. N(w) = NVAL 477: * 478: IF( ITMP1.LE.NVAL( JI ) ) THEN 479: AB( JI, 1 ) = TMP1 480: NAB( JI, 1 ) = ITMP1 481: END IF 482: IF( ITMP1.GE.NVAL( JI ) ) THEN 483: AB( JI, 2 ) = TMP1 484: NAB( JI, 2 ) = ITMP1 485: END IF 486: END IF 487: 100 CONTINUE 488: KL = KLNEW 489: * 490: * End of Serial Version of the loop 491: * 492: END IF 493: * 494: * Check for convergence 495: * 496: KFNEW = KF 497: DO 110 JI = KF, KL 498: TMP1 = ABS( AB( JI, 2 )-AB( JI, 1 ) ) 499: TMP2 = MAX( ABS( AB( JI, 2 ) ), ABS( AB( JI, 1 ) ) ) 500: IF( TMP1.LT.MAX( ABSTOL, PIVMIN, RELTOL*TMP2 ) .OR. 501: $ NAB( JI, 1 ).GE.NAB( JI, 2 ) ) THEN 502: * 503: * Converged -- Swap with position KFNEW, 504: * then increment KFNEW 505: * 506: IF( JI.GT.KFNEW ) THEN 507: TMP1 = AB( JI, 1 ) 508: TMP2 = AB( JI, 2 ) 509: ITMP1 = NAB( JI, 1 ) 510: ITMP2 = NAB( JI, 2 ) 511: AB( JI, 1 ) = AB( KFNEW, 1 ) 512: AB( JI, 2 ) = AB( KFNEW, 2 ) 513: NAB( JI, 1 ) = NAB( KFNEW, 1 ) 514: NAB( JI, 2 ) = NAB( KFNEW, 2 ) 515: AB( KFNEW, 1 ) = TMP1 516: AB( KFNEW, 2 ) = TMP2 517: NAB( KFNEW, 1 ) = ITMP1 518: NAB( KFNEW, 2 ) = ITMP2 519: IF( IJOB.EQ.3 ) THEN 520: ITMP1 = NVAL( JI ) 521: NVAL( JI ) = NVAL( KFNEW ) 522: NVAL( KFNEW ) = ITMP1 523: END IF 524: END IF 525: KFNEW = KFNEW + 1 526: END IF 527: 110 CONTINUE 528: KF = KFNEW 529: * 530: * Choose Midpoints 531: * 532: DO 120 JI = KF, KL 533: C( JI ) = HALF*( AB( JI, 1 )+AB( JI, 2 ) ) 534: 120 CONTINUE 535: * 536: * If no more intervals to refine, quit. 537: * 538: IF( KF.GT.KL ) 539: $ GO TO 140 540: 130 CONTINUE 541: * 542: * Converged 543: * 544: 140 CONTINUE 545: INFO = MAX( KL+1-KF, 0 ) 546: MOUT = KL 547: * 548: RETURN 549: * 550: * End of DLAEBZ 551: * 552: END