 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 SUBROUTINE SLAEBZ( 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.3.1) -- 00006 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00007 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00008 * -- April 2011 -- 00009 * 00010 * .. Scalar Arguments .. 00011 INTEGER IJOB, INFO, MINP, MMAX, MOUT, N, NBMIN, NITMAX 00012 REAL ABSTOL, PIVMIN, RELTOL 00013 * .. 00014 * .. Array Arguments .. 00015 INTEGER IWORK( * ), NAB( MMAX, * ), NVAL( * ) 00016 REAL AB( MMAX, * ), C( * ), D( * ), E( * ), E2( * ), 00017 $ WORK( * ) 00018 * .. 00019 * 00020 * Purpose 00021 * ======= 00022 * 00023 * SLAEBZ 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 SLAEBZ 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 SLAEBZ 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) REAL 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) REAL 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) REAL 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) REAL array, dimension (N) 00118 * The diagonal elements of the tridiagonal matrix T. 00119 * 00120 * E (input) REAL 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) REAL 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) REAL 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) REAL 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 SLAEBZ 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 SLAEBZ is called. 00176 * 00177 * WORK (workspace) REAL 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, SLAEBZ should have one or 00196 * more initial intervals set up in AB, and SLAEBZ 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. SLAEBZ 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). SLAEBZ 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 REAL ZERO, TWO, HALF 00226 PARAMETER ( ZERO = 0.0E0, TWO = 2.0E0, 00227 $ HALF = 1.0E0 / TWO ) 00228 * .. 00229 * .. Local Scalars .. 00230 INTEGER ITMP1, ITMP2, J, JI, JIT, JP, KF, KFNEW, KL, 00231 $ KLNEW 00232 REAL 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 DO 30 JI = 1, MINP 00255 DO 20 JP = 1, 2 00256 TMP1 = D( 1 ) - AB( JI, JP ) 00257 IF( ABS( TMP1 ).LT.PIVMIN ) 00258 $ TMP1 = -PIVMIN 00259 NAB( JI, JP ) = 0 00260 IF( TMP1.LE.ZERO ) 00261 $ NAB( JI, JP ) = 1 00262 * 00263 DO 10 J = 2, N 00264 TMP1 = D( J ) - E2( J-1 ) / TMP1 - AB( JI, JP ) 00265 IF( ABS( TMP1 ).LT.PIVMIN ) 00266 $ TMP1 = -PIVMIN 00267 IF( TMP1.LE.ZERO ) 00268 $ NAB( JI, JP ) = NAB( JI, JP ) + 1 00269 10 CONTINUE 00270 20 CONTINUE 00271 MOUT = MOUT + NAB( JI, 2 ) - NAB( JI, 1 ) 00272 30 CONTINUE 00273 RETURN 00274 END IF 00275 * 00276 * Initialize for loop 00277 * 00278 * KF and KL have the following meaning: 00279 * Intervals 1,...,KF-1 have converged. 00280 * Intervals KF,...,KL still need to be refined. 00281 * 00282 KF = 1 00283 KL = MINP 00284 * 00285 * If IJOB=2, initialize C. 00286 * If IJOB=3, use the user-supplied starting point. 00287 * 00288 IF( IJOB.EQ.2 ) THEN 00289 DO 40 JI = 1, MINP 00290 C( JI ) = HALF*( AB( JI, 1 )+AB( JI, 2 ) ) 00291 40 CONTINUE 00292 END IF 00293 * 00294 * Iteration loop 00295 * 00296 DO 130 JIT = 1, NITMAX 00297 * 00298 * Loop over intervals 00299 * 00300 IF( KL-KF+1.GE.NBMIN .AND. NBMIN.GT.0 ) THEN 00301 * 00302 * Begin of Parallel Version of the loop 00303 * 00304 DO 60 JI = KF, KL 00305 * 00306 * Compute N(c), the number of eigenvalues less than c 00307 * 00308 WORK( JI ) = D( 1 ) - C( JI ) 00309 IWORK( JI ) = 0 00310 IF( WORK( JI ).LE.PIVMIN ) THEN 00311 IWORK( JI ) = 1 00312 WORK( JI ) = MIN( WORK( JI ), -PIVMIN ) 00313 END IF 00314 * 00315 DO 50 J = 2, N 00316 WORK( JI ) = D( J ) - E2( J-1 ) / WORK( JI ) - C( JI ) 00317 IF( WORK( JI ).LE.PIVMIN ) THEN 00318 IWORK( JI ) = IWORK( JI ) + 1 00319 WORK( JI ) = MIN( WORK( JI ), -PIVMIN ) 00320 END IF 00321 50 CONTINUE 00322 60 CONTINUE 00323 * 00324 IF( IJOB.LE.2 ) THEN 00325 * 00326 * IJOB=2: Choose all intervals containing eigenvalues. 00327 * 00328 KLNEW = KL 00329 DO 70 JI = KF, KL 00330 * 00331 * Insure that N(w) is monotone 00332 * 00333 IWORK( JI ) = MIN( NAB( JI, 2 ), 00334 $ MAX( NAB( JI, 1 ), IWORK( JI ) ) ) 00335 * 00336 * Update the Queue -- add intervals if both halves 00337 * contain eigenvalues. 00338 * 00339 IF( IWORK( JI ).EQ.NAB( JI, 2 ) ) THEN 00340 * 00341 * No eigenvalue in the upper interval: 00342 * just use the lower interval. 00343 * 00344 AB( JI, 2 ) = C( JI ) 00345 * 00346 ELSE IF( IWORK( JI ).EQ.NAB( JI, 1 ) ) THEN 00347 * 00348 * No eigenvalue in the lower interval: 00349 * just use the upper interval. 00350 * 00351 AB( JI, 1 ) = C( JI ) 00352 ELSE 00353 KLNEW = KLNEW + 1 00354 IF( KLNEW.LE.MMAX ) THEN 00355 * 00356 * Eigenvalue in both intervals -- add upper to 00357 * queue. 00358 * 00359 AB( KLNEW, 2 ) = AB( JI, 2 ) 00360 NAB( KLNEW, 2 ) = NAB( JI, 2 ) 00361 AB( KLNEW, 1 ) = C( JI ) 00362 NAB( KLNEW, 1 ) = IWORK( JI ) 00363 AB( JI, 2 ) = C( JI ) 00364 NAB( JI, 2 ) = IWORK( JI ) 00365 ELSE 00366 INFO = MMAX + 1 00367 END IF 00368 END IF 00369 70 CONTINUE 00370 IF( INFO.NE.0 ) 00371 $ RETURN 00372 KL = KLNEW 00373 ELSE 00374 * 00375 * IJOB=3: Binary search. Keep only the interval containing 00376 * w s.t. N(w) = NVAL 00377 * 00378 DO 80 JI = KF, KL 00379 IF( IWORK( JI ).LE.NVAL( JI ) ) THEN 00380 AB( JI, 1 ) = C( JI ) 00381 NAB( JI, 1 ) = IWORK( JI ) 00382 END IF 00383 IF( IWORK( JI ).GE.NVAL( JI ) ) THEN 00384 AB( JI, 2 ) = C( JI ) 00385 NAB( JI, 2 ) = IWORK( JI ) 00386 END IF 00387 80 CONTINUE 00388 END IF 00389 * 00390 ELSE 00391 * 00392 * End of Parallel Version of the loop 00393 * 00394 * Begin of Serial Version of the loop 00395 * 00396 KLNEW = KL 00397 DO 100 JI = KF, KL 00398 * 00399 * Compute N(w), the number of eigenvalues less than w 00400 * 00401 TMP1 = C( JI ) 00402 TMP2 = D( 1 ) - TMP1 00403 ITMP1 = 0 00404 IF( TMP2.LE.PIVMIN ) THEN 00405 ITMP1 = 1 00406 TMP2 = MIN( TMP2, -PIVMIN ) 00407 END IF 00408 * 00409 DO 90 J = 2, N 00410 TMP2 = D( J ) - E2( J-1 ) / TMP2 - TMP1 00411 IF( TMP2.LE.PIVMIN ) THEN 00412 ITMP1 = ITMP1 + 1 00413 TMP2 = MIN( TMP2, -PIVMIN ) 00414 END IF 00415 90 CONTINUE 00416 * 00417 IF( IJOB.LE.2 ) THEN 00418 * 00419 * IJOB=2: Choose all intervals containing eigenvalues. 00420 * 00421 * Insure that N(w) is monotone 00422 * 00423 ITMP1 = MIN( NAB( JI, 2 ), 00424 $ MAX( NAB( JI, 1 ), ITMP1 ) ) 00425 * 00426 * Update the Queue -- add intervals if both halves 00427 * contain eigenvalues. 00428 * 00429 IF( ITMP1.EQ.NAB( JI, 2 ) ) THEN 00430 * 00431 * No eigenvalue in the upper interval: 00432 * just use the lower interval. 00433 * 00434 AB( JI, 2 ) = TMP1 00435 * 00436 ELSE IF( ITMP1.EQ.NAB( JI, 1 ) ) THEN 00437 * 00438 * No eigenvalue in the lower interval: 00439 * just use the upper interval. 00440 * 00441 AB( JI, 1 ) = TMP1 00442 ELSE IF( KLNEW.LT.MMAX ) THEN 00443 * 00444 * Eigenvalue in both intervals -- add upper to queue. 00445 * 00446 KLNEW = KLNEW + 1 00447 AB( KLNEW, 2 ) = AB( JI, 2 ) 00448 NAB( KLNEW, 2 ) = NAB( JI, 2 ) 00449 AB( KLNEW, 1 ) = TMP1 00450 NAB( KLNEW, 1 ) = ITMP1 00451 AB( JI, 2 ) = TMP1 00452 NAB( JI, 2 ) = ITMP1 00453 ELSE 00454 INFO = MMAX + 1 00455 RETURN 00456 END IF 00457 ELSE 00458 * 00459 * IJOB=3: Binary search. Keep only the interval 00460 * containing w s.t. N(w) = NVAL 00461 * 00462 IF( ITMP1.LE.NVAL( JI ) ) THEN 00463 AB( JI, 1 ) = TMP1 00464 NAB( JI, 1 ) = ITMP1 00465 END IF 00466 IF( ITMP1.GE.NVAL( JI ) ) THEN 00467 AB( JI, 2 ) = TMP1 00468 NAB( JI, 2 ) = ITMP1 00469 END IF 00470 END IF 00471 100 CONTINUE 00472 KL = KLNEW 00473 * 00474 END IF 00475 * 00476 * Check for convergence 00477 * 00478 KFNEW = KF 00479 DO 110 JI = KF, KL 00480 TMP1 = ABS( AB( JI, 2 )-AB( JI, 1 ) ) 00481 TMP2 = MAX( ABS( AB( JI, 2 ) ), ABS( AB( JI, 1 ) ) ) 00482 IF( TMP1.LT.MAX( ABSTOL, PIVMIN, RELTOL*TMP2 ) .OR. 00483 $ NAB( JI, 1 ).GE.NAB( JI, 2 ) ) THEN 00484 * 00485 * Converged -- Swap with position KFNEW, 00486 * then increment KFNEW 00487 * 00488 IF( JI.GT.KFNEW ) THEN 00489 TMP1 = AB( JI, 1 ) 00490 TMP2 = AB( JI, 2 ) 00491 ITMP1 = NAB( JI, 1 ) 00492 ITMP2 = NAB( JI, 2 ) 00493 AB( JI, 1 ) = AB( KFNEW, 1 ) 00494 AB( JI, 2 ) = AB( KFNEW, 2 ) 00495 NAB( JI, 1 ) = NAB( KFNEW, 1 ) 00496 NAB( JI, 2 ) = NAB( KFNEW, 2 ) 00497 AB( KFNEW, 1 ) = TMP1 00498 AB( KFNEW, 2 ) = TMP2 00499 NAB( KFNEW, 1 ) = ITMP1 00500 NAB( KFNEW, 2 ) = ITMP2 00501 IF( IJOB.EQ.3 ) THEN 00502 ITMP1 = NVAL( JI ) 00503 NVAL( JI ) = NVAL( KFNEW ) 00504 NVAL( KFNEW ) = ITMP1 00505 END IF 00506 END IF 00507 KFNEW = KFNEW + 1 00508 END IF 00509 110 CONTINUE 00510 KF = KFNEW 00511 * 00512 * Choose Midpoints 00513 * 00514 DO 120 JI = KF, KL 00515 C( JI ) = HALF*( AB( JI, 1 )+AB( JI, 2 ) ) 00516 120 CONTINUE 00517 * 00518 * If no more intervals to refine, quit. 00519 * 00520 IF( KF.GT.KL ) 00521 $ GO TO 140 00522 130 CONTINUE 00523 * 00524 * Converged 00525 * 00526 140 CONTINUE 00527 INFO = MAX( KL+1-KF, 0 ) 00528 MOUT = KL 00529 * 00530 RETURN 00531 * 00532 * End of SLAEBZ 00533 * 00534 END
1.7.3 
