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