00001 SUBROUTINE DLAEBZ( IJOB, NITMAX, N, MMAX, MINP, NBMIN, ABSTOL, 00002 $ RELTOL, PIVMIN, D, E, E2, NVAL, AB, C, MOUT, 00003 $ NAB, WORK, IWORK, INFO ) 00004 * 00005 * -- LAPACK auxiliary routine (version 3.2) -- 00006 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00007 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00008 * November 2006 00009 * 00010 * .. Scalar Arguments .. 00011 INTEGER IJOB, INFO, MINP, MMAX, MOUT, N, NBMIN, NITMAX 00012 DOUBLE PRECISION ABSTOL, PIVMIN, RELTOL 00013 * .. 00014 * .. Array Arguments .. 00015 INTEGER IWORK( * ), NAB( MMAX, * ), NVAL( * ) 00016 DOUBLE PRECISION AB( MMAX, * ), C( * ), D( * ), E( * ), E2( * ), 00017 $ WORK( * ) 00018 * .. 00019 * 00020 * Purpose 00021 * ======= 00022 * 00023 * DLAEBZ contains the iteration loops which compute and use the 00024 * function N(w), which is the count of eigenvalues of a symmetric 00025 * tridiagonal matrix T less than or equal to its argument w. It 00026 * performs a choice of two types of loops: 00027 * 00028 * IJOB=1, followed by 00029 * IJOB=2: It takes as input a list of intervals and returns a list of 00030 * sufficiently small intervals whose union contains the same 00031 * eigenvalues as the union of the original intervals. 00032 * The input intervals are (AB(j,1),AB(j,2)], j=1,...,MINP. 00033 * The output interval (AB(j,1),AB(j,2)] will contain 00034 * eigenvalues NAB(j,1)+1,...,NAB(j,2), where 1 <= j <= MOUT. 00035 * 00036 * IJOB=3: It performs a binary search in each input interval 00037 * (AB(j,1),AB(j,2)] for a point w(j) such that 00038 * N(w(j))=NVAL(j), and uses C(j) as the starting point of 00039 * the search. If such a w(j) is found, then on output 00040 * AB(j,1)=AB(j,2)=w. If no such w(j) is found, then on output 00041 * (AB(j,1),AB(j,2)] will be a small interval containing the 00042 * point where N(w) jumps through NVAL(j), unless that point 00043 * lies outside the initial interval. 00044 * 00045 * Note that the intervals are in all cases half-open intervals, 00046 * i.e., of the form (a,b] , which includes b but not a . 00047 * 00048 * To avoid underflow, the matrix should be scaled so that its largest 00049 * element is no greater than overflow**(1/2) * underflow**(1/4) 00050 * in absolute value. To assure the most accurate computation 00051 * of small eigenvalues, the matrix should be scaled to be 00052 * not much smaller than that, either. 00053 * 00054 * See W. Kahan "Accurate Eigenvalues of a Symmetric Tridiagonal 00055 * Matrix", Report CS41, Computer Science Dept., Stanford 00056 * University, July 21, 1966 00057 * 00058 * Note: the arguments are, in general, *not* checked for unreasonable 00059 * values. 00060 * 00061 * Arguments 00062 * ========= 00063 * 00064 * IJOB (input) INTEGER 00065 * Specifies what is to be done: 00066 * = 1: Compute NAB for the initial intervals. 00067 * = 2: Perform bisection iteration to find eigenvalues of T. 00068 * = 3: Perform bisection iteration to invert N(w), i.e., 00069 * to find a point which has a specified number of 00070 * eigenvalues of T to its left. 00071 * Other values will cause DLAEBZ to return with INFO=-1. 00072 * 00073 * NITMAX (input) INTEGER 00074 * The maximum number of "levels" of bisection to be 00075 * performed, i.e., an interval of width W will not be made 00076 * smaller than 2^(-NITMAX) * W. If not all intervals 00077 * have converged after NITMAX iterations, then INFO is set 00078 * to the number of non-converged intervals. 00079 * 00080 * N (input) INTEGER 00081 * The dimension n of the tridiagonal matrix T. It must be at 00082 * least 1. 00083 * 00084 * MMAX (input) INTEGER 00085 * The maximum number of intervals. If more than MMAX intervals 00086 * are generated, then DLAEBZ will quit with INFO=MMAX+1. 00087 * 00088 * MINP (input) INTEGER 00089 * The initial number of intervals. It may not be greater than 00090 * MMAX. 00091 * 00092 * NBMIN (input) INTEGER 00093 * The smallest number of intervals that should be processed 00094 * using a vector loop. If zero, then only the scalar loop 00095 * will be used. 00096 * 00097 * ABSTOL (input) DOUBLE PRECISION 00098 * The minimum (absolute) width of an interval. When an 00099 * interval is narrower than ABSTOL, or than RELTOL times the 00100 * larger (in magnitude) endpoint, then it is considered to be 00101 * sufficiently small, i.e., converged. This must be at least 00102 * zero. 00103 * 00104 * RELTOL (input) DOUBLE PRECISION 00105 * The minimum relative width of an interval. When an interval 00106 * is narrower than ABSTOL, or than RELTOL times the larger (in 00107 * magnitude) endpoint, then it is considered to be 00108 * sufficiently small, i.e., converged. Note: this should 00109 * always be at least radix*machine epsilon. 00110 * 00111 * PIVMIN (input) DOUBLE PRECISION 00112 * The minimum absolute value of a "pivot" in the Sturm 00113 * sequence loop. This *must* be at least max |e(j)**2| * 00114 * safe_min and at least safe_min, where safe_min is at least 00115 * the smallest number that can divide one without overflow. 00116 * 00117 * D (input) DOUBLE PRECISION array, dimension (N) 00118 * The diagonal elements of the tridiagonal matrix T. 00119 * 00120 * E (input) DOUBLE PRECISION array, dimension (N) 00121 * The offdiagonal elements of the tridiagonal matrix T in 00122 * positions 1 through N-1. E(N) is arbitrary. 00123 * 00124 * E2 (input) DOUBLE PRECISION array, dimension (N) 00125 * The squares of the offdiagonal elements of the tridiagonal 00126 * matrix T. E2(N) is ignored. 00127 * 00128 * NVAL (input/output) INTEGER array, dimension (MINP) 00129 * If IJOB=1 or 2, not referenced. 00130 * If IJOB=3, the desired values of N(w). The elements of NVAL 00131 * will be reordered to correspond with the intervals in AB. 00132 * Thus, NVAL(j) on output will not, in general be the same as 00133 * NVAL(j) on input, but it will correspond with the interval 00134 * (AB(j,1),AB(j,2)] on output. 00135 * 00136 * AB (input/output) DOUBLE PRECISION array, dimension (MMAX,2) 00137 * The endpoints of the intervals. AB(j,1) is a(j), the left 00138 * endpoint of the j-th interval, and AB(j,2) is b(j), the 00139 * right endpoint of the j-th interval. The input intervals 00140 * will, in general, be modified, split, and reordered by the 00141 * calculation. 00142 * 00143 * C (input/output) DOUBLE PRECISION array, dimension (MMAX) 00144 * If IJOB=1, ignored. 00145 * If IJOB=2, workspace. 00146 * If IJOB=3, then on input C(j) should be initialized to the 00147 * first search point in the binary search. 00148 * 00149 * MOUT (output) INTEGER 00150 * If IJOB=1, the number of eigenvalues in the intervals. 00151 * If IJOB=2 or 3, the number of intervals output. 00152 * If IJOB=3, MOUT will equal MINP. 00153 * 00154 * NAB (input/output) INTEGER array, dimension (MMAX,2) 00155 * If IJOB=1, then on output NAB(i,j) will be set to N(AB(i,j)). 00156 * If IJOB=2, then on input, NAB(i,j) should be set. It must 00157 * satisfy the condition: 00158 * N(AB(i,1)) <= NAB(i,1) <= NAB(i,2) <= N(AB(i,2)), 00159 * which means that in interval i only eigenvalues 00160 * NAB(i,1)+1,...,NAB(i,2) will be considered. Usually, 00161 * NAB(i,j)=N(AB(i,j)), from a previous call to DLAEBZ with 00162 * IJOB=1. 00163 * On output, NAB(i,j) will contain 00164 * max(na(k),min(nb(k),N(AB(i,j)))), where k is the index of 00165 * the input interval that the output interval 00166 * (AB(j,1),AB(j,2)] came from, and na(k) and nb(k) are the 00167 * the input values of NAB(k,1) and NAB(k,2). 00168 * If IJOB=3, then on output, NAB(i,j) contains N(AB(i,j)), 00169 * unless N(w) > NVAL(i) for all search points w , in which 00170 * case NAB(i,1) will not be modified, i.e., the output 00171 * value will be the same as the input value (modulo 00172 * reorderings -- see NVAL and AB), or unless N(w) < NVAL(i) 00173 * for all search points w , in which case NAB(i,2) will 00174 * not be modified. Normally, NAB should be set to some 00175 * distinctive value(s) before DLAEBZ is called. 00176 * 00177 * WORK (workspace) DOUBLE PRECISION array, dimension (MMAX) 00178 * Workspace. 00179 * 00180 * IWORK (workspace) INTEGER array, dimension (MMAX) 00181 * Workspace. 00182 * 00183 * INFO (output) INTEGER 00184 * = 0: All intervals converged. 00185 * = 1--MMAX: The last INFO intervals did not converge. 00186 * = MMAX+1: More than MMAX intervals were generated. 00187 * 00188 * Further Details 00189 * =============== 00190 * 00191 * This routine is intended to be called only by other LAPACK 00192 * routines, thus the interface is less user-friendly. It is intended 00193 * for two purposes: 00194 * 00195 * (a) finding eigenvalues. In this case, DLAEBZ should have one or 00196 * more initial intervals set up in AB, and DLAEBZ should be called 00197 * with IJOB=1. This sets up NAB, and also counts the eigenvalues. 00198 * Intervals with no eigenvalues would usually be thrown out at 00199 * this point. Also, if not all the eigenvalues in an interval i 00200 * are desired, NAB(i,1) can be increased or NAB(i,2) decreased. 00201 * For example, set NAB(i,1)=NAB(i,2)-1 to get the largest 00202 * eigenvalue. DLAEBZ is then called with IJOB=2 and MMAX 00203 * no smaller than the value of MOUT returned by the call with 00204 * IJOB=1. After this (IJOB=2) call, eigenvalues NAB(i,1)+1 00205 * through NAB(i,2) are approximately AB(i,1) (or AB(i,2)) to the 00206 * tolerance specified by ABSTOL and RELTOL. 00207 * 00208 * (b) finding an interval (a',b'] containing eigenvalues w(f),...,w(l). 00209 * In this case, start with a Gershgorin interval (a,b). Set up 00210 * AB to contain 2 search intervals, both initially (a,b). One 00211 * NVAL element should contain f-1 and the other should contain l 00212 * , while C should contain a and b, resp. NAB(i,1) should be -1 00213 * and NAB(i,2) should be N+1, to flag an error if the desired 00214 * interval does not lie in (a,b). DLAEBZ is then called with 00215 * IJOB=3. On exit, if w(f-1) < w(f), then one of the intervals -- 00216 * j -- will have AB(j,1)=AB(j,2) and NAB(j,1)=NAB(j,2)=f-1, while 00217 * if, to the specified tolerance, w(f-k)=...=w(f+r), k > 0 and r 00218 * >= 0, then the interval will have N(AB(j,1))=NAB(j,1)=f-k and 00219 * N(AB(j,2))=NAB(j,2)=f+r. The cases w(l) < w(l+1) and 00220 * w(l-r)=...=w(l+k) are handled similarly. 00221 * 00222 * ===================================================================== 00223 * 00224 * .. Parameters .. 00225 DOUBLE PRECISION ZERO, TWO, HALF 00226 PARAMETER ( ZERO = 0.0D0, TWO = 2.0D0, 00227 $ HALF = 1.0D0 / TWO ) 00228 * .. 00229 * .. Local Scalars .. 00230 INTEGER ITMP1, ITMP2, J, JI, JIT, JP, KF, KFNEW, KL, 00231 $ KLNEW 00232 DOUBLE PRECISION TMP1, TMP2 00233 * .. 00234 * .. Intrinsic Functions .. 00235 INTRINSIC ABS, MAX, MIN 00236 * .. 00237 * .. Executable Statements .. 00238 * 00239 * Check for Errors 00240 * 00241 INFO = 0 00242 IF( IJOB.LT.1 .OR. IJOB.GT.3 ) THEN 00243 INFO = -1 00244 RETURN 00245 END IF 00246 * 00247 * Initialize NAB 00248 * 00249 IF( IJOB.EQ.1 ) THEN 00250 * 00251 * Compute the number of eigenvalues in the initial intervals. 00252 * 00253 MOUT = 0 00254 *DIR$ NOVECTOR 00255 DO 30 JI = 1, MINP 00256 DO 20 JP = 1, 2 00257 TMP1 = D( 1 ) - AB( JI, JP ) 00258 IF( ABS( TMP1 ).LT.PIVMIN ) 00259 $ TMP1 = -PIVMIN 00260 NAB( JI, JP ) = 0 00261 IF( TMP1.LE.ZERO ) 00262 $ NAB( JI, JP ) = 1 00263 * 00264 DO 10 J = 2, N 00265 TMP1 = D( J ) - E2( J-1 ) / TMP1 - AB( JI, JP ) 00266 IF( ABS( TMP1 ).LT.PIVMIN ) 00267 $ TMP1 = -PIVMIN 00268 IF( TMP1.LE.ZERO ) 00269 $ NAB( JI, JP ) = NAB( JI, JP ) + 1 00270 10 CONTINUE 00271 20 CONTINUE 00272 MOUT = MOUT + NAB( JI, 2 ) - NAB( JI, 1 ) 00273 30 CONTINUE 00274 RETURN 00275 END IF 00276 * 00277 * Initialize for loop 00278 * 00279 * KF and KL have the following meaning: 00280 * Intervals 1,...,KF-1 have converged. 00281 * Intervals KF,...,KL still need to be refined. 00282 * 00283 KF = 1 00284 KL = MINP 00285 * 00286 * If IJOB=2, initialize C. 00287 * If IJOB=3, use the user-supplied starting point. 00288 * 00289 IF( IJOB.EQ.2 ) THEN 00290 DO 40 JI = 1, MINP 00291 C( JI ) = HALF*( AB( JI, 1 )+AB( JI, 2 ) ) 00292 40 CONTINUE 00293 END IF 00294 * 00295 * Iteration loop 00296 * 00297 DO 130 JIT = 1, NITMAX 00298 * 00299 * Loop over intervals 00300 * 00301 IF( KL-KF+1.GE.NBMIN .AND. NBMIN.GT.0 ) THEN 00302 * 00303 * Begin of Parallel Version of the loop 00304 * 00305 DO 60 JI = KF, KL 00306 * 00307 * Compute N(c), the number of eigenvalues less than c 00308 * 00309 WORK( JI ) = D( 1 ) - C( JI ) 00310 IWORK( JI ) = 0 00311 IF( WORK( JI ).LE.PIVMIN ) THEN 00312 IWORK( JI ) = 1 00313 WORK( JI ) = MIN( WORK( JI ), -PIVMIN ) 00314 END IF 00315 * 00316 DO 50 J = 2, N 00317 WORK( JI ) = D( J ) - E2( J-1 ) / WORK( JI ) - C( JI ) 00318 IF( WORK( JI ).LE.PIVMIN ) THEN 00319 IWORK( JI ) = IWORK( JI ) + 1 00320 WORK( JI ) = MIN( WORK( JI ), -PIVMIN ) 00321 END IF 00322 50 CONTINUE 00323 60 CONTINUE 00324 * 00325 IF( IJOB.LE.2 ) THEN 00326 * 00327 * IJOB=2: Choose all intervals containing eigenvalues. 00328 * 00329 KLNEW = KL 00330 DO 70 JI = KF, KL 00331 * 00332 * Insure that N(w) is monotone 00333 * 00334 IWORK( JI ) = MIN( NAB( JI, 2 ), 00335 $ MAX( NAB( JI, 1 ), IWORK( JI ) ) ) 00336 * 00337 * Update the Queue -- add intervals if both halves 00338 * contain eigenvalues. 00339 * 00340 IF( IWORK( JI ).EQ.NAB( JI, 2 ) ) THEN 00341 * 00342 * No eigenvalue in the upper interval: 00343 * just use the lower interval. 00344 * 00345 AB( JI, 2 ) = C( JI ) 00346 * 00347 ELSE IF( IWORK( JI ).EQ.NAB( JI, 1 ) ) THEN 00348 * 00349 * No eigenvalue in the lower interval: 00350 * just use the upper interval. 00351 * 00352 AB( JI, 1 ) = C( JI ) 00353 ELSE 00354 KLNEW = KLNEW + 1 00355 IF( KLNEW.LE.MMAX ) THEN 00356 * 00357 * Eigenvalue in both intervals -- add upper to 00358 * queue. 00359 * 00360 AB( KLNEW, 2 ) = AB( JI, 2 ) 00361 NAB( KLNEW, 2 ) = NAB( JI, 2 ) 00362 AB( KLNEW, 1 ) = C( JI ) 00363 NAB( KLNEW, 1 ) = IWORK( JI ) 00364 AB( JI, 2 ) = C( JI ) 00365 NAB( JI, 2 ) = IWORK( JI ) 00366 ELSE 00367 INFO = MMAX + 1 00368 END IF 00369 END IF 00370 70 CONTINUE 00371 IF( INFO.NE.0 ) 00372 $ RETURN 00373 KL = KLNEW 00374 ELSE 00375 * 00376 * IJOB=3: Binary search. Keep only the interval containing 00377 * w s.t. N(w) = NVAL 00378 * 00379 DO 80 JI = KF, KL 00380 IF( IWORK( JI ).LE.NVAL( JI ) ) THEN 00381 AB( JI, 1 ) = C( JI ) 00382 NAB( JI, 1 ) = IWORK( JI ) 00383 END IF 00384 IF( IWORK( JI ).GE.NVAL( JI ) ) THEN 00385 AB( JI, 2 ) = C( JI ) 00386 NAB( JI, 2 ) = IWORK( JI ) 00387 END IF 00388 80 CONTINUE 00389 END IF 00390 * 00391 ELSE 00392 * 00393 * End of Parallel Version of the loop 00394 * 00395 * Begin of Serial Version of the loop 00396 * 00397 KLNEW = KL 00398 DO 100 JI = KF, KL 00399 * 00400 * Compute N(w), the number of eigenvalues less than w 00401 * 00402 TMP1 = C( JI ) 00403 TMP2 = D( 1 ) - TMP1 00404 ITMP1 = 0 00405 IF( TMP2.LE.PIVMIN ) THEN 00406 ITMP1 = 1 00407 TMP2 = MIN( TMP2, -PIVMIN ) 00408 END IF 00409 * 00410 * A series of compiler directives to defeat vectorization 00411 * for the next loop 00412 * 00413 *$PL$ CMCHAR=' ' 00414 CDIR$ NEXTSCALAR 00415 C$DIR SCALAR 00416 CDIR$ NEXT SCALAR 00417 CVD$L NOVECTOR 00418 CDEC$ NOVECTOR 00419 CVD$ NOVECTOR 00420 *VDIR NOVECTOR 00421 *VOCL LOOP,SCALAR 00422 CIBM PREFER SCALAR 00423 *$PL$ CMCHAR='*' 00424 * 00425 DO 90 J = 2, N 00426 TMP2 = D( J ) - E2( J-1 ) / TMP2 - TMP1 00427 IF( TMP2.LE.PIVMIN ) THEN 00428 ITMP1 = ITMP1 + 1 00429 TMP2 = MIN( TMP2, -PIVMIN ) 00430 END IF 00431 90 CONTINUE 00432 * 00433 IF( IJOB.LE.2 ) THEN 00434 * 00435 * IJOB=2: Choose all intervals containing eigenvalues. 00436 * 00437 * Insure that N(w) is monotone 00438 * 00439 ITMP1 = MIN( NAB( JI, 2 ), 00440 $ MAX( NAB( JI, 1 ), ITMP1 ) ) 00441 * 00442 * Update the Queue -- add intervals if both halves 00443 * contain eigenvalues. 00444 * 00445 IF( ITMP1.EQ.NAB( JI, 2 ) ) THEN 00446 * 00447 * No eigenvalue in the upper interval: 00448 * just use the lower interval. 00449 * 00450 AB( JI, 2 ) = TMP1 00451 * 00452 ELSE IF( ITMP1.EQ.NAB( JI, 1 ) ) THEN 00453 * 00454 * No eigenvalue in the lower interval: 00455 * just use the upper interval. 00456 * 00457 AB( JI, 1 ) = TMP1 00458 ELSE IF( KLNEW.LT.MMAX ) THEN 00459 * 00460 * Eigenvalue in both intervals -- add upper to queue. 00461 * 00462 KLNEW = KLNEW + 1 00463 AB( KLNEW, 2 ) = AB( JI, 2 ) 00464 NAB( KLNEW, 2 ) = NAB( JI, 2 ) 00465 AB( KLNEW, 1 ) = TMP1 00466 NAB( KLNEW, 1 ) = ITMP1 00467 AB( JI, 2 ) = TMP1 00468 NAB( JI, 2 ) = ITMP1 00469 ELSE 00470 INFO = MMAX + 1 00471 RETURN 00472 END IF 00473 ELSE 00474 * 00475 * IJOB=3: Binary search. Keep only the interval 00476 * containing w s.t. N(w) = NVAL 00477 * 00478 IF( ITMP1.LE.NVAL( JI ) ) THEN 00479 AB( JI, 1 ) = TMP1 00480 NAB( JI, 1 ) = ITMP1 00481 END IF 00482 IF( ITMP1.GE.NVAL( JI ) ) THEN 00483 AB( JI, 2 ) = TMP1 00484 NAB( JI, 2 ) = ITMP1 00485 END IF 00486 END IF 00487 100 CONTINUE 00488 KL = KLNEW 00489 * 00490 * End of Serial Version of the loop 00491 * 00492 END IF 00493 * 00494 * Check for convergence 00495 * 00496 KFNEW = KF 00497 DO 110 JI = KF, KL 00498 TMP1 = ABS( AB( JI, 2 )-AB( JI, 1 ) ) 00499 TMP2 = MAX( ABS( AB( JI, 2 ) ), ABS( AB( JI, 1 ) ) ) 00500 IF( TMP1.LT.MAX( ABSTOL, PIVMIN, RELTOL*TMP2 ) .OR. 00501 $ NAB( JI, 1 ).GE.NAB( JI, 2 ) ) THEN 00502 * 00503 * Converged -- Swap with position KFNEW, 00504 * then increment KFNEW 00505 * 00506 IF( JI.GT.KFNEW ) THEN 00507 TMP1 = AB( JI, 1 ) 00508 TMP2 = AB( JI, 2 ) 00509 ITMP1 = NAB( JI, 1 ) 00510 ITMP2 = NAB( JI, 2 ) 00511 AB( JI, 1 ) = AB( KFNEW, 1 ) 00512 AB( JI, 2 ) = AB( KFNEW, 2 ) 00513 NAB( JI, 1 ) = NAB( KFNEW, 1 ) 00514 NAB( JI, 2 ) = NAB( KFNEW, 2 ) 00515 AB( KFNEW, 1 ) = TMP1 00516 AB( KFNEW, 2 ) = TMP2 00517 NAB( KFNEW, 1 ) = ITMP1 00518 NAB( KFNEW, 2 ) = ITMP2 00519 IF( IJOB.EQ.3 ) THEN 00520 ITMP1 = NVAL( JI ) 00521 NVAL( JI ) = NVAL( KFNEW ) 00522 NVAL( KFNEW ) = ITMP1 00523 END IF 00524 END IF 00525 KFNEW = KFNEW + 1 00526 END IF 00527 110 CONTINUE 00528 KF = KFNEW 00529 * 00530 * Choose Midpoints 00531 * 00532 DO 120 JI = KF, KL 00533 C( JI ) = HALF*( AB( JI, 1 )+AB( JI, 2 ) ) 00534 120 CONTINUE 00535 * 00536 * If no more intervals to refine, quit. 00537 * 00538 IF( KF.GT.KL ) 00539 $ GO TO 140 00540 130 CONTINUE 00541 * 00542 * Converged 00543 * 00544 140 CONTINUE 00545 INFO = MAX( KL+1-KF, 0 ) 00546 MOUT = KL 00547 * 00548 RETURN 00549 * 00550 * End of DLAEBZ 00551 * 00552 END
1.7.2