![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b SLARRJ 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download SLARRJ + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slarrj.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slarrj.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slarrj.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE SLARRJ( N, D, E2, IFIRST, ILAST, 00022 * RTOL, OFFSET, W, WERR, WORK, IWORK, 00023 * PIVMIN, SPDIAM, INFO ) 00024 * 00025 * .. Scalar Arguments .. 00026 * INTEGER IFIRST, ILAST, INFO, N, OFFSET 00027 * REAL PIVMIN, RTOL, SPDIAM 00028 * .. 00029 * .. Array Arguments .. 00030 * INTEGER IWORK( * ) 00031 * REAL D( * ), E2( * ), W( * ), 00032 * $ WERR( * ), WORK( * ) 00033 * .. 00034 * 00035 * 00036 *> \par Purpose: 00037 * ============= 00038 *> 00039 *> \verbatim 00040 *> 00041 *> Given the initial eigenvalue approximations of T, SLARRJ 00042 *> does bisection to refine the eigenvalues of T, 00043 *> W( IFIRST-OFFSET ) through W( ILAST-OFFSET ), to more accuracy. Initial 00044 *> guesses for these eigenvalues are input in W, the corresponding estimate 00045 *> of the error in these guesses in WERR. During bisection, intervals 00046 *> [left, right] are maintained by storing their mid-points and 00047 *> semi-widths in the arrays W and WERR respectively. 00048 *> \endverbatim 00049 * 00050 * Arguments: 00051 * ========== 00052 * 00053 *> \param[in] N 00054 *> \verbatim 00055 *> N is INTEGER 00056 *> The order of the matrix. 00057 *> \endverbatim 00058 *> 00059 *> \param[in] D 00060 *> \verbatim 00061 *> D is REAL array, dimension (N) 00062 *> The N diagonal elements of T. 00063 *> \endverbatim 00064 *> 00065 *> \param[in] E2 00066 *> \verbatim 00067 *> E2 is REAL array, dimension (N-1) 00068 *> The Squares of the (N-1) subdiagonal elements of T. 00069 *> \endverbatim 00070 *> 00071 *> \param[in] IFIRST 00072 *> \verbatim 00073 *> IFIRST is INTEGER 00074 *> The index of the first eigenvalue to be computed. 00075 *> \endverbatim 00076 *> 00077 *> \param[in] ILAST 00078 *> \verbatim 00079 *> ILAST is INTEGER 00080 *> The index of the last eigenvalue to be computed. 00081 *> \endverbatim 00082 *> 00083 *> \param[in] RTOL 00084 *> \verbatim 00085 *> RTOL is REAL 00086 *> Tolerance for the convergence of the bisection intervals. 00087 *> An interval [LEFT,RIGHT] has converged if 00088 *> RIGHT-LEFT.LT.RTOL*MAX(|LEFT|,|RIGHT|). 00089 *> \endverbatim 00090 *> 00091 *> \param[in] OFFSET 00092 *> \verbatim 00093 *> OFFSET is INTEGER 00094 *> Offset for the arrays W and WERR, i.e., the IFIRST-OFFSET 00095 *> through ILAST-OFFSET elements of these arrays are to be used. 00096 *> \endverbatim 00097 *> 00098 *> \param[in,out] W 00099 *> \verbatim 00100 *> W is REAL array, dimension (N) 00101 *> On input, W( IFIRST-OFFSET ) through W( ILAST-OFFSET ) are 00102 *> estimates of the eigenvalues of L D L^T indexed IFIRST through 00103 *> ILAST. 00104 *> On output, these estimates are refined. 00105 *> \endverbatim 00106 *> 00107 *> \param[in,out] WERR 00108 *> \verbatim 00109 *> WERR is REAL array, dimension (N) 00110 *> On input, WERR( IFIRST-OFFSET ) through WERR( ILAST-OFFSET ) are 00111 *> the errors in the estimates of the corresponding elements in W. 00112 *> On output, these errors are refined. 00113 *> \endverbatim 00114 *> 00115 *> \param[out] WORK 00116 *> \verbatim 00117 *> WORK is REAL array, dimension (2*N) 00118 *> Workspace. 00119 *> \endverbatim 00120 *> 00121 *> \param[out] IWORK 00122 *> \verbatim 00123 *> IWORK is INTEGER array, dimension (2*N) 00124 *> Workspace. 00125 *> \endverbatim 00126 *> 00127 *> \param[in] PIVMIN 00128 *> \verbatim 00129 *> PIVMIN is REAL 00130 *> The minimum pivot in the Sturm sequence for T. 00131 *> \endverbatim 00132 *> 00133 *> \param[in] SPDIAM 00134 *> \verbatim 00135 *> SPDIAM is REAL 00136 *> The spectral diameter of T. 00137 *> \endverbatim 00138 *> 00139 *> \param[out] INFO 00140 *> \verbatim 00141 *> INFO is INTEGER 00142 *> Error flag. 00143 *> \endverbatim 00144 * 00145 * Authors: 00146 * ======== 00147 * 00148 *> \author Univ. of Tennessee 00149 *> \author Univ. of California Berkeley 00150 *> \author Univ. of Colorado Denver 00151 *> \author NAG Ltd. 00152 * 00153 *> \date November 2011 00154 * 00155 *> \ingroup auxOTHERauxiliary 00156 * 00157 *> \par Contributors: 00158 * ================== 00159 *> 00160 *> Beresford Parlett, University of California, Berkeley, USA \n 00161 *> Jim Demmel, University of California, Berkeley, USA \n 00162 *> Inderjit Dhillon, University of Texas, Austin, USA \n 00163 *> Osni Marques, LBNL/NERSC, USA \n 00164 *> Christof Voemel, University of California, Berkeley, USA 00165 * 00166 * ===================================================================== 00167 SUBROUTINE SLARRJ( N, D, E2, IFIRST, ILAST, 00168 $ RTOL, OFFSET, W, WERR, WORK, IWORK, 00169 $ PIVMIN, SPDIAM, INFO ) 00170 * 00171 * -- LAPACK auxiliary routine (version 3.4.0) -- 00172 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00173 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00174 * November 2011 00175 * 00176 * .. Scalar Arguments .. 00177 INTEGER IFIRST, ILAST, INFO, N, OFFSET 00178 REAL PIVMIN, RTOL, SPDIAM 00179 * .. 00180 * .. Array Arguments .. 00181 INTEGER IWORK( * ) 00182 REAL D( * ), E2( * ), W( * ), 00183 $ WERR( * ), WORK( * ) 00184 * .. 00185 * 00186 * ===================================================================== 00187 * 00188 * .. Parameters .. 00189 REAL ZERO, ONE, TWO, HALF 00190 PARAMETER ( ZERO = 0.0E0, ONE = 1.0E0, TWO = 2.0E0, 00191 $ HALF = 0.5E0 ) 00192 INTEGER MAXITR 00193 * .. 00194 * .. Local Scalars .. 00195 INTEGER CNT, I, I1, I2, II, ITER, J, K, NEXT, NINT, 00196 $ OLNINT, P, PREV, SAVI1 00197 REAL DPLUS, FAC, LEFT, MID, RIGHT, S, TMP, WIDTH 00198 * 00199 * .. 00200 * .. Intrinsic Functions .. 00201 INTRINSIC ABS, MAX 00202 * .. 00203 * .. Executable Statements .. 00204 * 00205 INFO = 0 00206 * 00207 MAXITR = INT( ( LOG( SPDIAM+PIVMIN )-LOG( PIVMIN ) ) / 00208 $ LOG( TWO ) ) + 2 00209 * 00210 * Initialize unconverged intervals in [ WORK(2*I-1), WORK(2*I) ]. 00211 * The Sturm Count, Count( WORK(2*I-1) ) is arranged to be I-1, while 00212 * Count( WORK(2*I) ) is stored in IWORK( 2*I ). The integer IWORK( 2*I-1 ) 00213 * for an unconverged interval is set to the index of the next unconverged 00214 * interval, and is -1 or 0 for a converged interval. Thus a linked 00215 * list of unconverged intervals is set up. 00216 * 00217 00218 I1 = IFIRST 00219 I2 = ILAST 00220 * The number of unconverged intervals 00221 NINT = 0 00222 * The last unconverged interval found 00223 PREV = 0 00224 DO 75 I = I1, I2 00225 K = 2*I 00226 II = I - OFFSET 00227 LEFT = W( II ) - WERR( II ) 00228 MID = W(II) 00229 RIGHT = W( II ) + WERR( II ) 00230 WIDTH = RIGHT - MID 00231 TMP = MAX( ABS( LEFT ), ABS( RIGHT ) ) 00232 00233 * The following test prevents the test of converged intervals 00234 IF( WIDTH.LT.RTOL*TMP ) THEN 00235 * This interval has already converged and does not need refinement. 00236 * (Note that the gaps might change through refining the 00237 * eigenvalues, however, they can only get bigger.) 00238 * Remove it from the list. 00239 IWORK( K-1 ) = -1 00240 * Make sure that I1 always points to the first unconverged interval 00241 IF((I.EQ.I1).AND.(I.LT.I2)) I1 = I + 1 00242 IF((PREV.GE.I1).AND.(I.LE.I2)) IWORK( 2*PREV-1 ) = I + 1 00243 ELSE 00244 * unconverged interval found 00245 PREV = I 00246 * Make sure that [LEFT,RIGHT] contains the desired eigenvalue 00247 * 00248 * Do while( CNT(LEFT).GT.I-1 ) 00249 * 00250 FAC = ONE 00251 20 CONTINUE 00252 CNT = 0 00253 S = LEFT 00254 DPLUS = D( 1 ) - S 00255 IF( DPLUS.LT.ZERO ) CNT = CNT + 1 00256 DO 30 J = 2, N 00257 DPLUS = D( J ) - S - E2( J-1 )/DPLUS 00258 IF( DPLUS.LT.ZERO ) CNT = CNT + 1 00259 30 CONTINUE 00260 IF( CNT.GT.I-1 ) THEN 00261 LEFT = LEFT - WERR( II )*FAC 00262 FAC = TWO*FAC 00263 GO TO 20 00264 END IF 00265 * 00266 * Do while( CNT(RIGHT).LT.I ) 00267 * 00268 FAC = ONE 00269 50 CONTINUE 00270 CNT = 0 00271 S = RIGHT 00272 DPLUS = D( 1 ) - S 00273 IF( DPLUS.LT.ZERO ) CNT = CNT + 1 00274 DO 60 J = 2, N 00275 DPLUS = D( J ) - S - E2( J-1 )/DPLUS 00276 IF( DPLUS.LT.ZERO ) CNT = CNT + 1 00277 60 CONTINUE 00278 IF( CNT.LT.I ) THEN 00279 RIGHT = RIGHT + WERR( II )*FAC 00280 FAC = TWO*FAC 00281 GO TO 50 00282 END IF 00283 NINT = NINT + 1 00284 IWORK( K-1 ) = I + 1 00285 IWORK( K ) = CNT 00286 END IF 00287 WORK( K-1 ) = LEFT 00288 WORK( K ) = RIGHT 00289 75 CONTINUE 00290 00291 00292 SAVI1 = I1 00293 * 00294 * Do while( NINT.GT.0 ), i.e. there are still unconverged intervals 00295 * and while (ITER.LT.MAXITR) 00296 * 00297 ITER = 0 00298 80 CONTINUE 00299 PREV = I1 - 1 00300 I = I1 00301 OLNINT = NINT 00302 00303 DO 100 P = 1, OLNINT 00304 K = 2*I 00305 II = I - OFFSET 00306 NEXT = IWORK( K-1 ) 00307 LEFT = WORK( K-1 ) 00308 RIGHT = WORK( K ) 00309 MID = HALF*( LEFT + RIGHT ) 00310 00311 * semiwidth of interval 00312 WIDTH = RIGHT - MID 00313 TMP = MAX( ABS( LEFT ), ABS( RIGHT ) ) 00314 00315 IF( ( WIDTH.LT.RTOL*TMP ) .OR. 00316 $ (ITER.EQ.MAXITR) )THEN 00317 * reduce number of unconverged intervals 00318 NINT = NINT - 1 00319 * Mark interval as converged. 00320 IWORK( K-1 ) = 0 00321 IF( I1.EQ.I ) THEN 00322 I1 = NEXT 00323 ELSE 00324 * Prev holds the last unconverged interval previously examined 00325 IF(PREV.GE.I1) IWORK( 2*PREV-1 ) = NEXT 00326 END IF 00327 I = NEXT 00328 GO TO 100 00329 END IF 00330 PREV = I 00331 * 00332 * Perform one bisection step 00333 * 00334 CNT = 0 00335 S = MID 00336 DPLUS = D( 1 ) - S 00337 IF( DPLUS.LT.ZERO ) CNT = CNT + 1 00338 DO 90 J = 2, N 00339 DPLUS = D( J ) - S - E2( J-1 )/DPLUS 00340 IF( DPLUS.LT.ZERO ) CNT = CNT + 1 00341 90 CONTINUE 00342 IF( CNT.LE.I-1 ) THEN 00343 WORK( K-1 ) = MID 00344 ELSE 00345 WORK( K ) = MID 00346 END IF 00347 I = NEXT 00348 00349 100 CONTINUE 00350 ITER = ITER + 1 00351 * do another loop if there are still unconverged intervals 00352 * However, in the last iteration, all intervals are accepted 00353 * since this is the best we can do. 00354 IF( ( NINT.GT.0 ).AND.(ITER.LE.MAXITR) ) GO TO 80 00355 * 00356 * 00357 * At this point, all the intervals have converged 00358 DO 110 I = SAVI1, ILAST 00359 K = 2*I 00360 II = I - OFFSET 00361 * All intervals marked by '0' have been refined. 00362 IF( IWORK( K-1 ).EQ.0 ) THEN 00363 W( II ) = HALF*( WORK( K-1 )+WORK( K ) ) 00364 WERR( II ) = WORK( K ) - W( II ) 00365 END IF 00366 110 CONTINUE 00367 * 00368 00369 RETURN 00370 * 00371 * End of SLARRJ 00372 * 00373 END