Go to the documentation of this file.
00001 *> \brief \b SLARRF
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *> \htmlonly
00009 *> Download SLARRF + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slarrf.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slarrf.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slarrf.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       SUBROUTINE SLARRF( N, D, L, LD, CLSTRT, CLEND,
00022 *                          W, WGAP, WERR,
00023 *                          SPDIAM, CLGAPL, CLGAPR, PIVMIN, SIGMA,
00024 *                          DPLUS, LPLUS, WORK, INFO )
00025 * 
00026 *       .. Scalar Arguments ..
00027 *       INTEGER            CLSTRT, CLEND, INFO, N
00028 *       REAL               CLGAPL, CLGAPR, PIVMIN, SIGMA, SPDIAM
00029 *       ..
00030 *       .. Array Arguments ..
00031 *       REAL               D( * ), DPLUS( * ), L( * ), LD( * ),
00032 *      $          LPLUS( * ), W( * ), WGAP( * ), WERR( * ), WORK( * )
00033 *       ..
00034 *  
00035 *
00036 *> \par Purpose:
00037 *  =============
00038 *>
00039 *> \verbatim
00040 *>
00041 *> Given the initial representation L D L^T and its cluster of close
00042 *> eigenvalues (in a relative measure), W( CLSTRT ), W( CLSTRT+1 ), ...
00043 *> W( CLEND ), SLARRF finds a new relatively robust representation
00044 *> L D L^T - SIGMA I = L(+) D(+) L(+)^T such that at least one of the
00045 *> eigenvalues of L(+) D(+) L(+)^T is relatively isolated.
00046 *> \endverbatim
00047 *
00048 *  Arguments:
00049 *  ==========
00050 *
00051 *> \param[in] N
00052 *> \verbatim
00053 *>          N is INTEGER
00054 *>          The order of the matrix (subblock, if the matrix splitted).
00055 *> \endverbatim
00056 *>
00057 *> \param[in] D
00058 *> \verbatim
00059 *>          D is REAL array, dimension (N)
00060 *>          The N diagonal elements of the diagonal matrix D.
00061 *> \endverbatim
00062 *>
00063 *> \param[in] L
00064 *> \verbatim
00065 *>          L is REAL array, dimension (N-1)
00066 *>          The (N-1) subdiagonal elements of the unit bidiagonal
00067 *>          matrix L.
00068 *> \endverbatim
00069 *>
00070 *> \param[in] LD
00071 *> \verbatim
00072 *>          LD is REAL array, dimension (N-1)
00073 *>          The (N-1) elements L(i)*D(i).
00074 *> \endverbatim
00075 *>
00076 *> \param[in] CLSTRT
00077 *> \verbatim
00078 *>          CLSTRT is INTEGER
00079 *>          The index of the first eigenvalue in the cluster.
00080 *> \endverbatim
00081 *>
00082 *> \param[in] CLEND
00083 *> \verbatim
00084 *>          CLEND is INTEGER
00085 *>          The index of the last eigenvalue in the cluster.
00086 *> \endverbatim
00087 *>
00088 *> \param[in] W
00089 *> \verbatim
00090 *>          W is REAL array, dimension
00091 *>          dimension is >=  (CLEND-CLSTRT+1)
00092 *>          The eigenvalue APPROXIMATIONS of L D L^T in ascending order.
00093 *>          W( CLSTRT ) through W( CLEND ) form the cluster of relatively
00094 *>          close eigenalues.
00095 *> \endverbatim
00096 *>
00097 *> \param[in,out] WGAP
00098 *> \verbatim
00099 *>          WGAP is REAL array, dimension
00100 *>          dimension is >=  (CLEND-CLSTRT+1)
00101 *>          The separation from the right neighbor eigenvalue in W.
00102 *> \endverbatim
00103 *>
00104 *> \param[in] WERR
00105 *> \verbatim
00106 *>          WERR is REAL array, dimension
00107 *>          dimension is >=  (CLEND-CLSTRT+1)
00108 *>          WERR contain the semiwidth of the uncertainty
00109 *>          interval of the corresponding eigenvalue APPROXIMATION in W
00110 *> \endverbatim
00111 *>
00112 *> \param[in] SPDIAM
00113 *> \verbatim
00114 *>          SPDIAM is REAL
00115 *>          estimate of the spectral diameter obtained from the
00116 *>          Gerschgorin intervals
00117 *> \endverbatim
00118 *>
00119 *> \param[in] CLGAPL
00120 *> \verbatim
00121 *>          CLGAPL is REAL
00122 *> \endverbatim
00123 *>
00124 *> \param[in] CLGAPR
00125 *> \verbatim
00126 *>          CLGAPR is REAL
00127 *>          absolute gap on each end of the cluster.
00128 *>          Set by the calling routine to protect against shifts too close
00129 *>          to eigenvalues outside the cluster.
00130 *> \endverbatim
00131 *>
00132 *> \param[in] PIVMIN
00133 *> \verbatim
00134 *>          PIVMIN is REAL
00135 *>          The minimum pivot allowed in the Sturm sequence.
00136 *> \endverbatim
00137 *>
00138 *> \param[out] SIGMA
00139 *> \verbatim
00140 *>          SIGMA is REAL
00141 *>          The shift used to form L(+) D(+) L(+)^T.
00142 *> \endverbatim
00143 *>
00144 *> \param[out] DPLUS
00145 *> \verbatim
00146 *>          DPLUS is REAL array, dimension (N)
00147 *>          The N diagonal elements of the diagonal matrix D(+).
00148 *> \endverbatim
00149 *>
00150 *> \param[out] LPLUS
00151 *> \verbatim
00152 *>          LPLUS is REAL array, dimension (N-1)
00153 *>          The first (N-1) elements of LPLUS contain the subdiagonal
00154 *>          elements of the unit bidiagonal matrix L(+).
00155 *> \endverbatim
00156 *>
00157 *> \param[out] WORK
00158 *> \verbatim
00159 *>          WORK is REAL array, dimension (2*N)
00160 *>          Workspace.
00161 *> \endverbatim
00162 *>
00163 *> \param[out] INFO
00164 *> \verbatim
00165 *>          INFO is INTEGER
00166 *>          Signals processing OK (=0) or failure (=1)
00167 *> \endverbatim
00168 *
00169 *  Authors:
00170 *  ========
00171 *
00172 *> \author Univ. of Tennessee 
00173 *> \author Univ. of California Berkeley 
00174 *> \author Univ. of Colorado Denver 
00175 *> \author NAG Ltd. 
00176 *
00177 *> \date November 2011
00178 *
00179 *> \ingroup auxOTHERauxiliary
00180 *
00181 *> \par Contributors:
00182 *  ==================
00183 *>
00184 *> Beresford Parlett, University of California, Berkeley, USA \n
00185 *> Jim Demmel, University of California, Berkeley, USA \n
00186 *> Inderjit Dhillon, University of Texas, Austin, USA \n
00187 *> Osni Marques, LBNL/NERSC, USA \n
00188 *> Christof Voemel, University of California, Berkeley, USA
00189 *
00190 *  =====================================================================
00191       SUBROUTINE SLARRF( N, D, L, LD, CLSTRT, CLEND,
00192      $                   W, WGAP, WERR,
00193      $                   SPDIAM, CLGAPL, CLGAPR, PIVMIN, SIGMA,
00194      $                   DPLUS, LPLUS, WORK, INFO )
00195 *
00196 *  -- LAPACK auxiliary routine (version 3.4.0) --
00197 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00198 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00199 *     November 2011
00200 *
00201 *     .. Scalar Arguments ..
00202       INTEGER            CLSTRT, CLEND, INFO, N
00203       REAL               CLGAPL, CLGAPR, PIVMIN, SIGMA, SPDIAM
00204 *     ..
00205 *     .. Array Arguments ..
00206       REAL               D( * ), DPLUS( * ), L( * ), LD( * ),
00207      $          LPLUS( * ), W( * ), WGAP( * ), WERR( * ), WORK( * )
00208 *     ..
00209 *
00210 *  =====================================================================
00211 *
00212 *     .. Parameters ..
00213       REAL               FOUR, MAXGROWTH1, MAXGROWTH2, ONE, QUART, TWO,
00214      $                   ZERO
00215       PARAMETER          ( ZERO = 0.0E0, ONE = 1.0E0, TWO = 2.0E0,
00216      $                     FOUR = 4.0E0, QUART = 0.25E0,
00217      $                     MAXGROWTH1 = 8.E0,
00218      $                     MAXGROWTH2 = 8.E0 )
00219 *     ..
00220 *     .. Local Scalars ..
00221       LOGICAL   DORRR1, FORCER, NOFAIL, SAWNAN1, SAWNAN2, TRYRRR1
00222       INTEGER            I, INDX, KTRY, KTRYMAX, SLEFT, SRIGHT, SHIFT
00223       PARAMETER          ( KTRYMAX = 1, SLEFT = 1, SRIGHT = 2 )
00224       REAL               AVGAP, BESTSHIFT, CLWDTH, EPS, FACT, FAIL,
00225      $                   FAIL2, GROWTHBOUND, LDELTA, LDMAX, LSIGMA,
00226      $                   MAX1, MAX2, MINGAP, OLDP, PROD, RDELTA, RDMAX,
00227      $                   RRR1, RRR2, RSIGMA, S, SMLGROWTH, TMP, ZNM2
00228 *     ..
00229 *     .. External Functions ..
00230       LOGICAL SISNAN
00231       REAL               SLAMCH
00232       EXTERNAL           SISNAN, SLAMCH
00233 *     ..
00234 *     .. External Subroutines ..
00235       EXTERNAL           SCOPY
00236 *     ..
00237 *     .. Intrinsic Functions ..
00238       INTRINSIC          ABS
00239 *     ..
00240 *     .. Executable Statements ..
00241 *
00242       INFO = 0
00243       FACT = REAL(2**KTRYMAX)
00244       EPS = SLAMCH( 'Precision' )
00245       SHIFT = 0
00246       FORCER = .FALSE.
00247 
00248 
00249 *     Note that we cannot guarantee that for any of the shifts tried,
00250 *     the factorization has a small or even moderate element growth.
00251 *     There could be Ritz values at both ends of the cluster and despite
00252 *     backing off, there are examples where all factorizations tried
00253 *     (in IEEE mode, allowing zero pivots & infinities) have INFINITE
00254 *     element growth.
00255 *     For this reason, we should use PIVMIN in this subroutine so that at
00256 *     least the L D L^T factorization exists. It can be checked afterwards
00257 *     whether the element growth caused bad residuals/orthogonality.
00258 
00259 *     Decide whether the code should accept the best among all
00260 *     representations despite large element growth or signal INFO=1
00261       NOFAIL = .TRUE.
00262 *
00263 
00264 *     Compute the average gap length of the cluster
00265       CLWDTH = ABS(W(CLEND)-W(CLSTRT)) + WERR(CLEND) + WERR(CLSTRT)
00266       AVGAP = CLWDTH / REAL(CLEND-CLSTRT)
00267       MINGAP = MIN(CLGAPL, CLGAPR)
00268 *     Initial values for shifts to both ends of cluster
00269       LSIGMA = MIN(W( CLSTRT ),W( CLEND )) - WERR( CLSTRT )
00270       RSIGMA = MAX(W( CLSTRT ),W( CLEND )) + WERR( CLEND )
00271 
00272 *     Use a small fudge to make sure that we really shift to the outside
00273       LSIGMA = LSIGMA - ABS(LSIGMA)* TWO * EPS
00274       RSIGMA = RSIGMA + ABS(RSIGMA)* TWO * EPS
00275 
00276 *     Compute upper bounds for how much to back off the initial shifts
00277       LDMAX = QUART * MINGAP + TWO * PIVMIN
00278       RDMAX = QUART * MINGAP + TWO * PIVMIN
00279 
00280       LDELTA = MAX(AVGAP,WGAP( CLSTRT ))/FACT
00281       RDELTA = MAX(AVGAP,WGAP( CLEND-1 ))/FACT
00282 *
00283 *     Initialize the record of the best representation found
00284 *
00285       S = SLAMCH( 'S' )
00286       SMLGROWTH = ONE / S
00287       FAIL = REAL(N-1)*MINGAP/(SPDIAM*EPS)
00288       FAIL2 = REAL(N-1)*MINGAP/(SPDIAM*SQRT(EPS))
00289       BESTSHIFT = LSIGMA
00290 *
00291 *     while (KTRY <= KTRYMAX)
00292       KTRY = 0
00293       GROWTHBOUND = MAXGROWTH1*SPDIAM
00294 
00295  5    CONTINUE
00296       SAWNAN1 = .FALSE.
00297       SAWNAN2 = .FALSE.
00298 *     Ensure that we do not back off too much of the initial shifts
00299       LDELTA = MIN(LDMAX,LDELTA)
00300       RDELTA = MIN(RDMAX,RDELTA)
00301 
00302 *     Compute the element growth when shifting to both ends of the cluster
00303 *     accept the shift if there is no element growth at one of the two ends
00304 
00305 *     Left end
00306       S = -LSIGMA
00307       DPLUS( 1 ) = D( 1 ) + S
00308       IF(ABS(DPLUS(1)).LT.PIVMIN) THEN
00309          DPLUS(1) = -PIVMIN
00310 *        Need to set SAWNAN1 because refined RRR test should not be used
00311 *        in this case
00312          SAWNAN1 = .TRUE.
00313       ENDIF
00314       MAX1 = ABS( DPLUS( 1 ) )
00315       DO 6 I = 1, N - 1
00316          LPLUS( I ) = LD( I ) / DPLUS( I )
00317          S = S*LPLUS( I )*L( I ) - LSIGMA
00318          DPLUS( I+1 ) = D( I+1 ) + S
00319          IF(ABS(DPLUS(I+1)).LT.PIVMIN) THEN
00320             DPLUS(I+1) = -PIVMIN
00321 *           Need to set SAWNAN1 because refined RRR test should not be used
00322 *           in this case
00323             SAWNAN1 = .TRUE.
00324          ENDIF
00325          MAX1 = MAX( MAX1,ABS(DPLUS(I+1)) )
00326  6    CONTINUE
00327       SAWNAN1 = SAWNAN1 .OR.  SISNAN( MAX1 )
00328 
00329       IF( FORCER .OR.
00330      $   (MAX1.LE.GROWTHBOUND .AND. .NOT.SAWNAN1 ) ) THEN
00331          SIGMA = LSIGMA
00332          SHIFT = SLEFT
00333          GOTO 100
00334       ENDIF
00335 
00336 *     Right end
00337       S = -RSIGMA
00338       WORK( 1 ) = D( 1 ) + S
00339       IF(ABS(WORK(1)).LT.PIVMIN) THEN
00340          WORK(1) = -PIVMIN
00341 *        Need to set SAWNAN2 because refined RRR test should not be used
00342 *        in this case
00343          SAWNAN2 = .TRUE.
00344       ENDIF
00345       MAX2 = ABS( WORK( 1 ) )
00346       DO 7 I = 1, N - 1
00347          WORK( N+I ) = LD( I ) / WORK( I )
00348          S = S*WORK( N+I )*L( I ) - RSIGMA
00349          WORK( I+1 ) = D( I+1 ) + S
00350          IF(ABS(WORK(I+1)).LT.PIVMIN) THEN
00351             WORK(I+1) = -PIVMIN
00352 *           Need to set SAWNAN2 because refined RRR test should not be used
00353 *           in this case
00354             SAWNAN2 = .TRUE.
00355          ENDIF
00356          MAX2 = MAX( MAX2,ABS(WORK(I+1)) )
00357  7    CONTINUE
00358       SAWNAN2 = SAWNAN2 .OR.  SISNAN( MAX2 )
00359 
00360       IF( FORCER .OR.
00361      $   (MAX2.LE.GROWTHBOUND .AND. .NOT.SAWNAN2 ) ) THEN
00362          SIGMA = RSIGMA
00363          SHIFT = SRIGHT
00364          GOTO 100
00365       ENDIF
00366 *     If we are at this point, both shifts led to too much element growth
00367 
00368 *     Record the better of the two shifts (provided it didn't lead to NaN)
00369       IF(SAWNAN1.AND.SAWNAN2) THEN
00370 *        both MAX1 and MAX2 are NaN
00371          GOTO 50
00372       ELSE
00373          IF( .NOT.SAWNAN1 ) THEN
00374             INDX = 1
00375             IF(MAX1.LE.SMLGROWTH) THEN
00376                SMLGROWTH = MAX1
00377                BESTSHIFT = LSIGMA
00378             ENDIF
00379          ENDIF
00380          IF( .NOT.SAWNAN2 ) THEN
00381             IF(SAWNAN1 .OR. MAX2.LE.MAX1) INDX = 2
00382             IF(MAX2.LE.SMLGROWTH) THEN
00383                SMLGROWTH = MAX2
00384                BESTSHIFT = RSIGMA
00385             ENDIF
00386          ENDIF
00387       ENDIF
00388 
00389 *     If we are here, both the left and the right shift led to
00390 *     element growth. If the element growth is moderate, then
00391 *     we may still accept the representation, if it passes a
00392 *     refined test for RRR. This test supposes that no NaN occurred.
00393 *     Moreover, we use the refined RRR test only for isolated clusters.
00394       IF((CLWDTH.LT.MINGAP/REAL(128)) .AND.
00395      $   (MIN(MAX1,MAX2).LT.FAIL2)
00396      $  .AND.(.NOT.SAWNAN1).AND.(.NOT.SAWNAN2)) THEN
00397          DORRR1 = .TRUE.
00398       ELSE
00399          DORRR1 = .FALSE.
00400       ENDIF
00401       TRYRRR1 = .TRUE.
00402       IF( TRYRRR1 .AND. DORRR1 ) THEN
00403       IF(INDX.EQ.1) THEN
00404          TMP = ABS( DPLUS( N ) )
00405          ZNM2 = ONE
00406          PROD = ONE
00407          OLDP = ONE
00408          DO 15 I = N-1, 1, -1
00409             IF( PROD .LE. EPS ) THEN
00410                PROD =
00411      $         ((DPLUS(I+1)*WORK(N+I+1))/(DPLUS(I)*WORK(N+I)))*OLDP
00412             ELSE
00413                PROD = PROD*ABS(WORK(N+I))
00414             END IF
00415             OLDP = PROD
00416             ZNM2 = ZNM2 + PROD**2
00417             TMP = MAX( TMP, ABS( DPLUS( I ) * PROD ))
00418  15      CONTINUE
00419          RRR1 = TMP/( SPDIAM * SQRT( ZNM2 ) )
00420          IF (RRR1.LE.MAXGROWTH2) THEN
00421             SIGMA = LSIGMA
00422             SHIFT = SLEFT
00423             GOTO 100
00424          ENDIF
00425       ELSE IF(INDX.EQ.2) THEN
00426          TMP = ABS( WORK( N ) )
00427          ZNM2 = ONE
00428          PROD = ONE
00429          OLDP = ONE
00430          DO 16 I = N-1, 1, -1
00431             IF( PROD .LE. EPS ) THEN
00432                PROD = ((WORK(I+1)*LPLUS(I+1))/(WORK(I)*LPLUS(I)))*OLDP
00433             ELSE
00434                PROD = PROD*ABS(LPLUS(I))
00435             END IF
00436             OLDP = PROD
00437             ZNM2 = ZNM2 + PROD**2
00438             TMP = MAX( TMP, ABS( WORK( I ) * PROD ))
00439  16      CONTINUE
00440          RRR2 = TMP/( SPDIAM * SQRT( ZNM2 ) )
00441          IF (RRR2.LE.MAXGROWTH2) THEN
00442             SIGMA = RSIGMA
00443             SHIFT = SRIGHT
00444             GOTO 100
00445          ENDIF
00446       END IF
00447       ENDIF
00448 
00449  50   CONTINUE
00450 
00451       IF (KTRY.LT.KTRYMAX) THEN
00452 *        If we are here, both shifts failed also the RRR test.
00453 *        Back off to the outside
00454          LSIGMA = MAX( LSIGMA - LDELTA,
00455      $     LSIGMA - LDMAX)
00456          RSIGMA = MIN( RSIGMA + RDELTA,
00457      $     RSIGMA + RDMAX )
00458          LDELTA = TWO * LDELTA
00459          RDELTA = TWO * RDELTA
00460          KTRY = KTRY + 1
00461          GOTO 5
00462       ELSE
00463 *        None of the representations investigated satisfied our
00464 *        criteria. Take the best one we found.
00465          IF((SMLGROWTH.LT.FAIL).OR.NOFAIL) THEN
00466             LSIGMA = BESTSHIFT
00467             RSIGMA = BESTSHIFT
00468             FORCER = .TRUE.
00469             GOTO 5
00470          ELSE
00471             INFO = 1
00472             RETURN
00473          ENDIF
00474       END IF
00475 
00476  100  CONTINUE
00477       IF (SHIFT.EQ.SLEFT) THEN
00478       ELSEIF (SHIFT.EQ.SRIGHT) THEN
00479 *        store new L and D back into DPLUS, LPLUS
00480          CALL SCOPY( N, WORK, 1, DPLUS, 1 )
00481          CALL SCOPY( N-1, WORK(N+1), 1, LPLUS, 1 )
00482       ENDIF
00483 
00484       RETURN
00485 *
00486 *     End of SLARRF
00487 *
00488       END