LAPACK: SRC/dlarrf.f Source File

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