SRC/dlarrf.f

   1 *> \brief \b DLARRF finds a new relatively robust representation such that at least one of the eigenvalues is relatively isolated.
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download DLARRF + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlarrf.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlarrf.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlarrf.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE DLARRF( N, D, L, LD, CLSTRT, CLEND,
  22 *                          W, WGAP, WERR,
  23 *                          SPDIAM, CLGAPL, CLGAPR, PIVMIN, SIGMA,
  24 *                          DPLUS, LPLUS, WORK, INFO )
  25 *
  26 *       .. Scalar Arguments ..
  27 *       INTEGER            CLSTRT, CLEND, INFO, N
  28 *       DOUBLE PRECISION   CLGAPL, CLGAPR, PIVMIN, SIGMA, SPDIAM
  29 *       ..
  30 *       .. Array Arguments ..
  31 *       DOUBLE PRECISION   D( * ), DPLUS( * ), L( * ), LD( * ),
  32 *      $          LPLUS( * ), W( * ), WGAP( * ), WERR( * ), WORK( * )
  33 *       ..
  34 *
  35 *
  36 *> \par Purpose:
  37 *  =============
  38 *>
  39 *> \verbatim
  40 *>
  41 *> Given the initial representation L D L^T and its cluster of close
  42 *> eigenvalues (in a relative measure), W( CLSTRT ), W( CLSTRT+1 ), ...
  43 *> W( CLEND ), DLARRF finds a new relatively robust representation
  44 *> L D L^T - SIGMA I = L(+) D(+) L(+)^T such that at least one of the
  45 *> eigenvalues of L(+) D(+) L(+)^T is relatively isolated.
  46 *> \endverbatim
  47 *
  48 *  Arguments:
  49 *  ==========
  50 *
  51 *> \param[in] N
  52 *> \verbatim
  53 *>          N is INTEGER
  54 *>          The order of the matrix (subblock, if the matrix split).
  55 *> \endverbatim
  56 *>
  57 *> \param[in] D
  58 *> \verbatim
  59 *>          D is DOUBLE PRECISION array, dimension (N)
  60 *>          The N diagonal elements of the diagonal matrix D.
  61 *> \endverbatim
  62 *>
  63 *> \param[in] L
  64 *> \verbatim
  65 *>          L is DOUBLE PRECISION array, dimension (N-1)
  66 *>          The (N-1) subdiagonal elements of the unit bidiagonal
  67 *>          matrix L.
  68 *> \endverbatim
  69 *>
  70 *> \param[in] LD
  71 *> \verbatim
  72 *>          LD is DOUBLE PRECISION array, dimension (N-1)
  73 *>          The (N-1) elements L(i)*D(i).
  74 *> \endverbatim
  75 *>
  76 *> \param[in] CLSTRT
  77 *> \verbatim
  78 *>          CLSTRT is INTEGER
  79 *>          The index of the first eigenvalue in the cluster.
  80 *> \endverbatim
  81 *>
  82 *> \param[in] CLEND
  83 *> \verbatim
  84 *>          CLEND is INTEGER
  85 *>          The index of the last eigenvalue in the cluster.
  86 *> \endverbatim
  87 *>
  88 *> \param[in] W
  89 *> \verbatim
  90 *>          W is DOUBLE PRECISION array, dimension
  91 *>          dimension is >=  (CLEND-CLSTRT+1)
  92 *>          The eigenvalue APPROXIMATIONS of L D L^T in ascending order.
  93 *>          W( CLSTRT ) through W( CLEND ) form the cluster of relatively
  94 *>          close eigenalues.
  95 *> \endverbatim
  96 *>
  97 *> \param[in,out] WGAP
  98 *> \verbatim
  99 *>          WGAP is DOUBLE PRECISION array, dimension
 100 *>          dimension is >=  (CLEND-CLSTRT+1)
 101 *>          The separation from the right neighbor eigenvalue in W.
 102 *> \endverbatim
 103 *>
 104 *> \param[in] WERR
 105 *> \verbatim
 106 *>          WERR is DOUBLE PRECISION array, dimension
 107 *>          dimension is  >=  (CLEND-CLSTRT+1)
 108 *>          WERR contain the semiwidth of the uncertainty
 109 *>          interval of the corresponding eigenvalue APPROXIMATION in W
 110 *> \endverbatim
 111 *>
 112 *> \param[in] SPDIAM
 113 *> \verbatim
 114 *>          SPDIAM is DOUBLE PRECISION
 115 *>          estimate of the spectral diameter obtained from the
 116 *>          Gerschgorin intervals
 117 *> \endverbatim
 118 *>
 119 *> \param[in] CLGAPL
 120 *> \verbatim
 121 *>          CLGAPL is DOUBLE PRECISION
 122 *> \endverbatim
 123 *>
 124 *> \param[in] CLGAPR
 125 *> \verbatim
 126 *>          CLGAPR is DOUBLE PRECISION
 127 *>          absolute gap on each end of the cluster.
 128 *>          Set by the calling routine to protect against shifts too close
 129 *>          to eigenvalues outside the cluster.
 130 *> \endverbatim
 131 *>
 132 *> \param[in] PIVMIN
 133 *> \verbatim
 134 *>          PIVMIN is DOUBLE PRECISION
 135 *>          The minimum pivot allowed in the Sturm sequence.
 136 *> \endverbatim
 137 *>
 138 *> \param[out] SIGMA
 139 *> \verbatim
 140 *>          SIGMA is DOUBLE PRECISION
 141 *>          The shift used to form L(+) D(+) L(+)^T.
 142 *> \endverbatim
 143 *>
 144 *> \param[out] DPLUS
 145 *> \verbatim
 146 *>          DPLUS is DOUBLE PRECISION array, dimension (N)
 147 *>          The N diagonal elements of the diagonal matrix D(+).
 148 *> \endverbatim
 149 *>
 150 *> \param[out] LPLUS
 151 *> \verbatim
 152 *>          LPLUS is DOUBLE PRECISION array, dimension (N-1)
 153 *>          The first (N-1) elements of LPLUS contain the subdiagonal
 154 *>          elements of the unit bidiagonal matrix L(+).
 155 *> \endverbatim
 156 *>
 157 *> \param[out] WORK
 158 *> \verbatim
 159 *>          WORK is DOUBLE PRECISION array, dimension (2*N)
 160 *>          Workspace.
 161 *> \endverbatim
 162 *>
 163 *> \param[out] INFO
 164 *> \verbatim
 165 *>          INFO is INTEGER
 166 *>          Signals processing OK (=0) or failure (=1)
 167 *> \endverbatim
 168 *
 169 *  Authors:
 170 *  ========
 171 *
 172 *> \author Univ. of Tennessee
 173 *> \author Univ. of California Berkeley
 174 *> \author Univ. of Colorado Denver
 175 *> \author NAG Ltd.
 176 *
 177 *> \date June 2016
 178 *
 179 *> \ingroup OTHERauxiliary
 180 *
 181 *> \par Contributors:
 182 *  ==================
 183 *>
 184 *> Beresford Parlett, University of California, Berkeley, USA \n
 185 *> Jim Demmel, University of California, Berkeley, USA \n
 186 *> Inderjit Dhillon, University of Texas, Austin, USA \n
 187 *> Osni Marques, LBNL/NERSC, USA \n
 188 *> Christof Voemel, University of California, Berkeley, USA
 189 *
 190 *  =====================================================================
 191       SUBROUTINE DLARRF( N, D, L, LD, CLSTRT, CLEND,
 192      $                   W, WGAP, WERR,
 193      $                   SPDIAM, CLGAPL, CLGAPR, PIVMIN, SIGMA,
 194      $                   DPLUS, LPLUS, WORK, INFO )
 195 *
 196 *  -- LAPACK auxiliary routine (version 3.6.1) --
 197 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 198 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 199 *     June 2016
 200 *
 201 *     .. Scalar Arguments ..
 202       INTEGER            CLSTRT, CLEND, INFO, N
 203       DOUBLE PRECISION   CLGAPL, CLGAPR, PIVMIN, SIGMA, SPDIAM
 204 *     ..
 205 *     .. Array Arguments ..
 206       DOUBLE PRECISION   D( * ), DPLUS( * ), L( * ), LD( * ),
 207      $          LPLUS( * ), W( * ), WGAP( * ), WERR( * ), WORK( * )
 208 *     ..
 209 *
 210 *  =====================================================================
 211 *
 212 *     .. Parameters ..
 213       DOUBLE PRECISION   FOUR, MAXGROWTH1, MAXGROWTH2, ONE, QUART, TWO
 214       PARAMETER          ( ONE = 1.0D0, TWO = 2.0D0, FOUR = 4.0D0,
 215      $                     QUART = 0.25D0,
 216      $                     MAXGROWTH1 = 8.D0,
 217      $                     MAXGROWTH2 = 8.D0 )
 218 *     ..
 219 *     .. Local Scalars ..
 220       LOGICAL   DORRR1, FORCER, NOFAIL, SAWNAN1, SAWNAN2, TRYRRR1
 221       INTEGER            I, INDX, KTRY, KTRYMAX, SLEFT, SRIGHT, SHIFT
 222       PARAMETER          ( KTRYMAX = 1, SLEFT = 1, SRIGHT = 2 )
 223       DOUBLE PRECISION   AVGAP, BESTSHIFT, CLWDTH, EPS, FACT, FAIL,
 224      $                   FAIL2, GROWTHBOUND, LDELTA, LDMAX, LSIGMA,
 225      $                   MAX1, MAX2, MINGAP, OLDP, PROD, RDELTA, RDMAX,
 226      $                   RRR1, RRR2, RSIGMA, S, SMLGROWTH, TMP, ZNM2
 227 *     ..
 228 *     .. External Functions ..
 229       LOGICAL DISNAN
 230       DOUBLE PRECISION   DLAMCH
 231       EXTERNAL           DISNAN, DLAMCH
 232 *     ..
 233 *     .. External Subroutines ..
 234       EXTERNAL           DCOPY
 235 *     ..
 236 *     .. Intrinsic Functions ..
 237       INTRINSIC          ABS
 238 *     ..
 239 *     .. Executable Statements ..
 240 *
 241       INFO = 0
 242       FACT = DBLE(2**KTRYMAX)
 243       EPS = DLAMCH( 'Precision' )
 244       SHIFT = 0
 245       FORCER = .FALSE.
 246
 247
 248 *     Note that we cannot guarantee that for any of the shifts tried,
 249 *     the factorization has a small or even moderate element growth.
 250 *     There could be Ritz values at both ends of the cluster and despite
 251 *     backing off, there are examples where all factorizations tried
 252 *     (in IEEE mode, allowing zero pivots & infinities) have INFINITE
 253 *     element growth.
 254 *     For this reason, we should use PIVMIN in this subroutine so that at
 255 *     least the L D L^T factorization exists. It can be checked afterwards
 256 *     whether the element growth caused bad residuals/orthogonality.
 257
 258 *     Decide whether the code should accept the best among all
 259 *     representations despite large element growth or signal INFO=1
 260 *     Setting NOFAIL to .FALSE. for quick fix for bug 113
 261       NOFAIL = .FALSE.
 262 *
 263
 264 *     Compute the average gap length of the cluster
 265       CLWDTH = ABS(W(CLEND)-W(CLSTRT)) + WERR(CLEND) + WERR(CLSTRT)
 266       AVGAP = CLWDTH / DBLE(CLEND-CLSTRT)
 267       MINGAP = MIN(CLGAPL, CLGAPR)
 268 *     Initial values for shifts to both ends of cluster
 269       LSIGMA = MIN(W( CLSTRT ),W( CLEND )) - WERR( CLSTRT )
 270       RSIGMA = MAX(W( CLSTRT ),W( CLEND )) + WERR( CLEND )
 271
 272 *     Use a small fudge to make sure that we really shift to the outside
 273       LSIGMA = LSIGMA - ABS(LSIGMA)* FOUR * EPS
 274       RSIGMA = RSIGMA + ABS(RSIGMA)* FOUR * EPS
 275
 276 *     Compute upper bounds for how much to back off the initial shifts
 277       LDMAX = QUART * MINGAP + TWO * PIVMIN
 278       RDMAX = QUART * MINGAP + TWO * PIVMIN
 279
 280       LDELTA = MAX(AVGAP,WGAP( CLSTRT ))/FACT
 281       RDELTA = MAX(AVGAP,WGAP( CLEND-1 ))/FACT
 282 *
 283 *     Initialize the record of the best representation found
 284 *
 285       S = DLAMCH( 'S' )
 286       SMLGROWTH = ONE / S
 287       FAIL = DBLE(N-1)*MINGAP/(SPDIAM*EPS)
 288       FAIL2 = DBLE(N-1)*MINGAP/(SPDIAM*SQRT(EPS))
 289       BESTSHIFT = LSIGMA
 290 *
 291 *     while (KTRY <= KTRYMAX)
 292       KTRY = 0
 293       GROWTHBOUND = MAXGROWTH1*SPDIAM
 294
 295  5    CONTINUE
 296       SAWNAN1 = .FALSE.
 297       SAWNAN2 = .FALSE.
 298 *     Ensure that we do not back off too much of the initial shifts
 299       LDELTA = MIN(LDMAX,LDELTA)
 300       RDELTA = MIN(RDMAX,RDELTA)
 301
 302 *     Compute the element growth when shifting to both ends of the cluster
 303 *     accept the shift if there is no element growth at one of the two ends
 304
 305 *     Left end
 306       S = -LSIGMA
 307       DPLUS( 1 ) = D( 1 ) + S
 308       IF(ABS(DPLUS(1)).LT.PIVMIN) THEN
 309          DPLUS(1) = -PIVMIN
 310 *        Need to set SAWNAN1 because refined RRR test should not be used
 311 *        in this case
 312          SAWNAN1 = .TRUE.
 313       ENDIF
 314       MAX1 = ABS( DPLUS( 1 ) )
 315       DO 6 I = 1, N - 1
 316          LPLUS( I ) = LD( I ) / DPLUS( I )
 317          S = S*LPLUS( I )*L( I ) - LSIGMA
 318          DPLUS( I+1 ) = D( I+1 ) + S
 319          IF(ABS(DPLUS(I+1)).LT.PIVMIN) THEN
 320             DPLUS(I+1) = -PIVMIN
 321 *           Need to set SAWNAN1 because refined RRR test should not be used
 322 *           in this case
 323             SAWNAN1 = .TRUE.
 324          ENDIF
 325          MAX1 = MAX( MAX1,ABS(DPLUS(I+1)) )
 326  6    CONTINUE
 327       SAWNAN1 = SAWNAN1 .OR.  DISNAN( MAX1 )
 328
 329       IF( FORCER .OR.
 330      $   (MAX1.LE.GROWTHBOUND .AND. .NOT.SAWNAN1 ) ) THEN
 331          SIGMA = LSIGMA
 332          SHIFT = SLEFT
 333          GOTO 100
 334       ENDIF
 335
 336 *     Right end
 337       S = -RSIGMA
 338       WORK( 1 ) = D( 1 ) + S
 339       IF(ABS(WORK(1)).LT.PIVMIN) THEN
 340          WORK(1) = -PIVMIN
 341 *        Need to set SAWNAN2 because refined RRR test should not be used
 342 *        in this case
 343          SAWNAN2 = .TRUE.
 344       ENDIF
 345       MAX2 = ABS( WORK( 1 ) )
 346       DO 7 I = 1, N - 1
 347          WORK( N+I ) = LD( I ) / WORK( I )
 348          S = S*WORK( N+I )*L( I ) - RSIGMA
 349          WORK( I+1 ) = D( I+1 ) + S
 350          IF(ABS(WORK(I+1)).LT.PIVMIN) THEN
 351             WORK(I+1) = -PIVMIN
 352 *           Need to set SAWNAN2 because refined RRR test should not be used
 353 *           in this case
 354             SAWNAN2 = .TRUE.
 355          ENDIF
 356          MAX2 = MAX( MAX2,ABS(WORK(I+1)) )
 357  7    CONTINUE
 358       SAWNAN2 = SAWNAN2 .OR.  DISNAN( MAX2 )
 359
 360       IF( FORCER .OR.
 361      $   (MAX2.LE.GROWTHBOUND .AND. .NOT.SAWNAN2 ) ) THEN
 362          SIGMA = RSIGMA
 363          SHIFT = SRIGHT
 364          GOTO 100
 365       ENDIF
 366 *     If we are at this point, both shifts led to too much element growth
 367
 368 *     Record the better of the two shifts (provided it didn't lead to NaN)
 369       IF(SAWNAN1.AND.SAWNAN2) THEN
 370 *        both MAX1 and MAX2 are NaN
 371          GOTO 50
 372       ELSE
 373          IF( .NOT.SAWNAN1 ) THEN
 374             INDX = 1
 375             IF(MAX1.LE.SMLGROWTH) THEN
 376                SMLGROWTH = MAX1
 377                BESTSHIFT = LSIGMA
 378             ENDIF
 379          ENDIF
 380          IF( .NOT.SAWNAN2 ) THEN
 381             IF(SAWNAN1 .OR. MAX2.LE.MAX1) INDX = 2
 382             IF(MAX2.LE.SMLGROWTH) THEN
 383                SMLGROWTH = MAX2
 384                BESTSHIFT = RSIGMA
 385             ENDIF
 386          ENDIF
 387       ENDIF
 388
 389 *     If we are here, both the left and the right shift led to
 390 *     element growth. If the element growth is moderate, then
 391 *     we may still accept the representation, if it passes a
 392 *     refined test for RRR. This test supposes that no NaN occurred.
 393 *     Moreover, we use the refined RRR test only for isolated clusters.
 394       IF((CLWDTH.LT.MINGAP/DBLE(128)) .AND.
 395      $   (MIN(MAX1,MAX2).LT.FAIL2)
 396      $  .AND.(.NOT.SAWNAN1).AND.(.NOT.SAWNAN2)) THEN
 397          DORRR1 = .TRUE.
 398       ELSE
 399          DORRR1 = .FALSE.
 400       ENDIF
 401       TRYRRR1 = .TRUE.
 402       IF( TRYRRR1 .AND. DORRR1 ) THEN
 403       IF(INDX.EQ.1) THEN
 404          TMP = ABS( DPLUS( N ) )
 405          ZNM2 = ONE
 406          PROD = ONE
 407          OLDP = ONE
 408          DO 15 I = N-1, 1, -1
 409             IF( PROD .LE. EPS ) THEN
 410                PROD =
 411      $         ((DPLUS(I+1)*WORK(N+I+1))/(DPLUS(I)*WORK(N+I)))*OLDP
 412             ELSE
 413                PROD = PROD*ABS(WORK(N+I))
 414             END IF
 415             OLDP = PROD
 416             ZNM2 = ZNM2 + PROD**2
 417             TMP = MAX( TMP, ABS( DPLUS( I ) * PROD ))
 418  15      CONTINUE
 419          RRR1 = TMP/( SPDIAM * SQRT( ZNM2 ) )
 420          IF (RRR1.LE.MAXGROWTH2) THEN
 421             SIGMA = LSIGMA
 422             SHIFT = SLEFT
 423             GOTO 100
 424          ENDIF
 425       ELSE IF(INDX.EQ.2) THEN
 426          TMP = ABS( WORK( N ) )
 427          ZNM2 = ONE
 428          PROD = ONE
 429          OLDP = ONE
 430          DO 16 I = N-1, 1, -1
 431             IF( PROD .LE. EPS ) THEN
 432                PROD = ((WORK(I+1)*LPLUS(I+1))/(WORK(I)*LPLUS(I)))*OLDP
 433             ELSE
 434                PROD = PROD*ABS(LPLUS(I))
 435             END IF
 436             OLDP = PROD
 437             ZNM2 = ZNM2 + PROD**2
 438             TMP = MAX( TMP, ABS( WORK( I ) * PROD ))
 439  16      CONTINUE
 440          RRR2 = TMP/( SPDIAM * SQRT( ZNM2 ) )
 441          IF (RRR2.LE.MAXGROWTH2) THEN
 442             SIGMA = RSIGMA
 443             SHIFT = SRIGHT
 444             GOTO 100
 445          ENDIF
 446       END IF
 447       ENDIF
 448
 449  50   CONTINUE
 450
 451       IF (KTRY.LT.KTRYMAX) THEN
 452 *        If we are here, both shifts failed also the RRR test.
 453 *        Back off to the outside
 454          LSIGMA = MAX( LSIGMA - LDELTA,
 455      $     LSIGMA - LDMAX)
 456          RSIGMA = MIN( RSIGMA + RDELTA,
 457      $     RSIGMA + RDMAX )
 458          LDELTA = TWO * LDELTA
 459          RDELTA = TWO * RDELTA
 460          KTRY = KTRY + 1
 461          GOTO 5
 462       ELSE
 463 *        None of the representations investigated satisfied our
 464 *        criteria. Take the best one we found.
 465          IF((SMLGROWTH.LT.FAIL).OR.NOFAIL) THEN
 466             LSIGMA = BESTSHIFT
 467             RSIGMA = BESTSHIFT
 468             FORCER = .TRUE.
 469             GOTO 5
 470          ELSE
 471             INFO = 1
 472             RETURN
 473          ENDIF
 474       END IF
 475
 476  100  CONTINUE
 477       IF (SHIFT.EQ.SLEFT) THEN
 478       ELSEIF (SHIFT.EQ.SRIGHT) THEN
 479 *        store new L and D back into DPLUS, LPLUS
 480          CALL DCOPY( N, WORK, 1, DPLUS, 1 )
 481          CALL DCOPY( N-1, WORK(N+1), 1, LPLUS, 1 )
 482       ENDIF
 483
 484       RETURN
 485 *
 486 *     End of DLARRF
 487 *
 488       END