LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
slarrb.f
Go to the documentation of this file.
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
 All Files Functions