LAPACK: SRC/slarrf.f Source File

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