SRC/dlar1v.f

   1 *> \brief \b DLAR1V computes the (scaled) r-th column of the inverse of the submatrix in rows b1 through bn of the tridiagonal matrix LDLT - λI.
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download DLAR1V + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlar1v.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlar1v.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlar1v.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE DLAR1V( N, B1, BN, LAMBDA, D, L, LD, LLD,
  22 *                  PIVMIN, GAPTOL, Z, WANTNC, NEGCNT, ZTZ, MINGMA,
  23 *                  R, ISUPPZ, NRMINV, RESID, RQCORR, WORK )
  24 *
  25 *       .. Scalar Arguments ..
  26 *       LOGICAL            WANTNC
  27 *       INTEGER   B1, BN, N, NEGCNT, R
  28 *       DOUBLE PRECISION   GAPTOL, LAMBDA, MINGMA, NRMINV, PIVMIN, RESID,
  29 *      $                   RQCORR, ZTZ
  30 *       ..
  31 *       .. Array Arguments ..
  32 *       INTEGER            ISUPPZ( * )
  33 *       DOUBLE PRECISION   D( * ), L( * ), LD( * ), LLD( * ),
  34 *      $                  WORK( * )
  35 *       DOUBLE PRECISION Z( * )
  36 *       ..
  37 *
  38 *
  39 *> \par Purpose:
  40 *  =============
  41 *>
  42 *> \verbatim
  43 *>
  44 *> DLAR1V computes the (scaled) r-th column of the inverse of
  45 *> the sumbmatrix in rows B1 through BN of the tridiagonal matrix
  46 *> L D L**T - sigma I. When sigma is close to an eigenvalue, the
  47 *> computed vector is an accurate eigenvector. Usually, r corresponds
  48 *> to the index where the eigenvector is largest in magnitude.
  49 *> The following steps accomplish this computation :
  50 *> (a) Stationary qd transform,  L D L**T - sigma I = L(+) D(+) L(+)**T,
  51 *> (b) Progressive qd transform, L D L**T - sigma I = U(-) D(-) U(-)**T,
  52 *> (c) Computation of the diagonal elements of the inverse of
  53 *>     L D L**T - sigma I by combining the above transforms, and choosing
  54 *>     r as the index where the diagonal of the inverse is (one of the)
  55 *>     largest in magnitude.
  56 *> (d) Computation of the (scaled) r-th column of the inverse using the
  57 *>     twisted factorization obtained by combining the top part of the
  58 *>     the stationary and the bottom part of the progressive transform.
  59 *> \endverbatim
  60 *
  61 *  Arguments:
  62 *  ==========
  63 *
  64 *> \param[in] N
  65 *> \verbatim
  66 *>          N is INTEGER
  67 *>           The order of the matrix L D L**T.
  68 *> \endverbatim
  69 *>
  70 *> \param[in] B1
  71 *> \verbatim
  72 *>          B1 is INTEGER
  73 *>           First index of the submatrix of L D L**T.
  74 *> \endverbatim
  75 *>
  76 *> \param[in] BN
  77 *> \verbatim
  78 *>          BN is INTEGER
  79 *>           Last index of the submatrix of L D L**T.
  80 *> \endverbatim
  81 *>
  82 *> \param[in] LAMBDA
  83 *> \verbatim
  84 *>          LAMBDA is DOUBLE PRECISION
  85 *>           The shift. In order to compute an accurate eigenvector,
  86 *>           LAMBDA should be a good approximation to an eigenvalue
  87 *>           of L D L**T.
  88 *> \endverbatim
  89 *>
  90 *> \param[in] L
  91 *> \verbatim
  92 *>          L is DOUBLE PRECISION array, dimension (N-1)
  93 *>           The (n-1) subdiagonal elements of the unit bidiagonal matrix
  94 *>           L, in elements 1 to N-1.
  95 *> \endverbatim
  96 *>
  97 *> \param[in] D
  98 *> \verbatim
  99 *>          D is DOUBLE PRECISION array, dimension (N)
 100 *>           The n diagonal elements of the diagonal matrix D.
 101 *> \endverbatim
 102 *>
 103 *> \param[in] LD
 104 *> \verbatim
 105 *>          LD is DOUBLE PRECISION array, dimension (N-1)
 106 *>           The n-1 elements L(i)*D(i).
 107 *> \endverbatim
 108 *>
 109 *> \param[in] LLD
 110 *> \verbatim
 111 *>          LLD is DOUBLE PRECISION array, dimension (N-1)
 112 *>           The n-1 elements L(i)*L(i)*D(i).
 113 *> \endverbatim
 114 *>
 115 *> \param[in] PIVMIN
 116 *> \verbatim
 117 *>          PIVMIN is DOUBLE PRECISION
 118 *>           The minimum pivot in the Sturm sequence.
 119 *> \endverbatim
 120 *>
 121 *> \param[in] GAPTOL
 122 *> \verbatim
 123 *>          GAPTOL is DOUBLE PRECISION
 124 *>           Tolerance that indicates when eigenvector entries are negligible
 125 *>           w.r.t. their contribution to the residual.
 126 *> \endverbatim
 127 *>
 128 *> \param[in,out] Z
 129 *> \verbatim
 130 *>          Z is DOUBLE PRECISION array, dimension (N)
 131 *>           On input, all entries of Z must be set to 0.
 132 *>           On output, Z contains the (scaled) r-th column of the
 133 *>           inverse. The scaling is such that Z(R) equals 1.
 134 *> \endverbatim
 135 *>
 136 *> \param[in] WANTNC
 137 *> \verbatim
 138 *>          WANTNC is LOGICAL
 139 *>           Specifies whether NEGCNT has to be computed.
 140 *> \endverbatim
 141 *>
 142 *> \param[out] NEGCNT
 143 *> \verbatim
 144 *>          NEGCNT is INTEGER
 145 *>           If WANTNC is .TRUE. then NEGCNT = the number of pivots < pivmin
 146 *>           in the  matrix factorization L D L**T, and NEGCNT = -1 otherwise.
 147 *> \endverbatim
 148 *>
 149 *> \param[out] ZTZ
 150 *> \verbatim
 151 *>          ZTZ is DOUBLE PRECISION
 152 *>           The square of the 2-norm of Z.
 153 *> \endverbatim
 154 *>
 155 *> \param[out] MINGMA
 156 *> \verbatim
 157 *>          MINGMA is DOUBLE PRECISION
 158 *>           The reciprocal of the largest (in magnitude) diagonal
 159 *>           element of the inverse of L D L**T - sigma I.
 160 *> \endverbatim
 161 *>
 162 *> \param[in,out] R
 163 *> \verbatim
 164 *>          R is INTEGER
 165 *>           The twist index for the twisted factorization used to
 166 *>           compute Z.
 167 *>           On input, 0 <= R <= N. If R is input as 0, R is set to
 168 *>           the index where (L D L**T - sigma I)^{-1} is largest
 169 *>           in magnitude. If 1 <= R <= N, R is unchanged.
 170 *>           On output, R contains the twist index used to compute Z.
 171 *>           Ideally, R designates the position of the maximum entry in the
 172 *>           eigenvector.
 173 *> \endverbatim
 174 *>
 175 *> \param[out] ISUPPZ
 176 *> \verbatim
 177 *>          ISUPPZ is INTEGER array, dimension (2)
 178 *>           The support of the vector in Z, i.e., the vector Z is
 179 *>           nonzero only in elements ISUPPZ(1) through ISUPPZ( 2 ).
 180 *> \endverbatim
 181 *>
 182 *> \param[out] NRMINV
 183 *> \verbatim
 184 *>          NRMINV is DOUBLE PRECISION
 185 *>           NRMINV = 1/SQRT( ZTZ )
 186 *> \endverbatim
 187 *>
 188 *> \param[out] RESID
 189 *> \verbatim
 190 *>          RESID is DOUBLE PRECISION
 191 *>           The residual of the FP vector.
 192 *>           RESID = ABS( MINGMA )/SQRT( ZTZ )
 193 *> \endverbatim
 194 *>
 195 *> \param[out] RQCORR
 196 *> \verbatim
 197 *>          RQCORR is DOUBLE PRECISION
 198 *>           The Rayleigh Quotient correction to LAMBDA.
 199 *>           RQCORR = MINGMA*TMP
 200 *> \endverbatim
 201 *>
 202 *> \param[out] WORK
 203 *> \verbatim
 204 *>          WORK is DOUBLE PRECISION array, dimension (4*N)
 205 *> \endverbatim
 206 *
 207 *  Authors:
 208 *  ========
 209 *
 210 *> \author Univ. of Tennessee
 211 *> \author Univ. of California Berkeley
 212 *> \author Univ. of Colorado Denver
 213 *> \author NAG Ltd.
 214 *
 215 *> \date September 2012
 216 *
 217 *> \ingroup doubleOTHERauxiliary
 218 *
 219 *> \par Contributors:
 220 *  ==================
 221 *>
 222 *> Beresford Parlett, University of California, Berkeley, USA \n
 223 *> Jim Demmel, University of California, Berkeley, USA \n
 224 *> Inderjit Dhillon, University of Texas, Austin, USA \n
 225 *> Osni Marques, LBNL/NERSC, USA \n
 226 *> Christof Voemel, University of California, Berkeley, USA
 227 *
 228 *  =====================================================================
 229       SUBROUTINE DLAR1V( N, B1, BN, LAMBDA, D, L, LD, LLD,
 230      $           PIVMIN, GAPTOL, Z, WANTNC, NEGCNT, ZTZ, MINGMA,
 231      $           R, ISUPPZ, NRMINV, RESID, RQCORR, WORK )
 232 *
 233 *  -- LAPACK auxiliary routine (version 3.4.2) --
 234 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 235 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 236 *     September 2012
 237 *
 238 *     .. Scalar Arguments ..
 239       LOGICAL            WANTNC
 240       INTEGER   B1, BN, N, NEGCNT, R
 241       DOUBLE PRECISION   GAPTOL, LAMBDA, MINGMA, NRMINV, PIVMIN, RESID,
 242      $                   RQCORR, ZTZ
 243 *     ..
 244 *     .. Array Arguments ..
 245       INTEGER            ISUPPZ( * )
 246       DOUBLE PRECISION   D( * ), L( * ), LD( * ), LLD( * ),
 247      $                  WORK( * )
 248       DOUBLE PRECISION Z( * )
 249 *     ..
 250 *
 251 *  =====================================================================
 252 *
 253 *     .. Parameters ..
 254       DOUBLE PRECISION   ZERO, ONE
 255       PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0 )
 256
 257 *     ..
 258 *     .. Local Scalars ..
 259       LOGICAL            SAWNAN1, SAWNAN2
 260       INTEGER            I, INDLPL, INDP, INDS, INDUMN, NEG1, NEG2, R1,
 261      $                   R2
 262       DOUBLE PRECISION   DMINUS, DPLUS, EPS, S, TMP
 263 *     ..
 264 *     .. External Functions ..
 265       LOGICAL DISNAN
 266       DOUBLE PRECISION   DLAMCH
 267       EXTERNAL           DISNAN, DLAMCH
 268 *     ..
 269 *     .. Intrinsic Functions ..
 270       INTRINSIC          ABS
 271 *     ..
 272 *     .. Executable Statements ..
 273 *
 274       EPS = DLAMCH( 'Precision' )
 275
 276
 277       IF( R.EQ.0 ) THEN
 278          R1 = B1
 279          R2 = BN
 280       ELSE
 281          R1 = R
 282          R2 = R
 283       END IF
 284
 285 *     Storage for LPLUS
 286       INDLPL = 0
 287 *     Storage for UMINUS
 288       INDUMN = N
 289       INDS = 2*N + 1
 290       INDP = 3*N + 1
 291
 292       IF( B1.EQ.1 ) THEN
 293          WORK( INDS ) = ZERO
 294       ELSE
 295          WORK( INDS+B1-1 ) = LLD( B1-1 )
 296       END IF
 297
 298 *
 299 *     Compute the stationary transform (using the differential form)
 300 *     until the index R2.
 301 *
 302       SAWNAN1 = .FALSE.
 303       NEG1 = 0
 304       S = WORK( INDS+B1-1 ) - LAMBDA
 305       DO 50 I = B1, R1 - 1
 306          DPLUS = D( I ) + S
 307          WORK( INDLPL+I ) = LD( I ) / DPLUS
 308          IF(DPLUS.LT.ZERO) NEG1 = NEG1 + 1
 309          WORK( INDS+I ) = S*WORK( INDLPL+I )*L( I )
 310          S = WORK( INDS+I ) - LAMBDA
 311  50   CONTINUE
 312       SAWNAN1 = DISNAN( S )
 313       IF( SAWNAN1 ) GOTO 60
 314       DO 51 I = R1, R2 - 1
 315          DPLUS = D( I ) + S
 316          WORK( INDLPL+I ) = LD( I ) / DPLUS
 317          WORK( INDS+I ) = S*WORK( INDLPL+I )*L( I )
 318          S = WORK( INDS+I ) - LAMBDA
 319  51   CONTINUE
 320       SAWNAN1 = DISNAN( S )
 321 *
 322  60   CONTINUE
 323       IF( SAWNAN1 ) THEN
 324 *        Runs a slower version of the above loop if a NaN is detected
 325          NEG1 = 0
 326          S = WORK( INDS+B1-1 ) - LAMBDA
 327          DO 70 I = B1, R1 - 1
 328             DPLUS = D( I ) + S
 329             IF(ABS(DPLUS).LT.PIVMIN) DPLUS = -PIVMIN
 330             WORK( INDLPL+I ) = LD( I ) / DPLUS
 331             IF(DPLUS.LT.ZERO) NEG1 = NEG1 + 1
 332             WORK( INDS+I ) = S*WORK( INDLPL+I )*L( I )
 333             IF( WORK( INDLPL+I ).EQ.ZERO )
 334      $                      WORK( INDS+I ) = LLD( I )
 335             S = WORK( INDS+I ) - LAMBDA
 336  70      CONTINUE
 337          DO 71 I = R1, R2 - 1
 338             DPLUS = D( I ) + S
 339             IF(ABS(DPLUS).LT.PIVMIN) DPLUS = -PIVMIN
 340             WORK( INDLPL+I ) = LD( I ) / DPLUS
 341             WORK( INDS+I ) = S*WORK( INDLPL+I )*L( I )
 342             IF( WORK( INDLPL+I ).EQ.ZERO )
 343      $                      WORK( INDS+I ) = LLD( I )
 344             S = WORK( INDS+I ) - LAMBDA
 345  71      CONTINUE
 346       END IF
 347 *
 348 *     Compute the progressive transform (using the differential form)
 349 *     until the index R1
 350 *
 351       SAWNAN2 = .FALSE.
 352       NEG2 = 0
 353       WORK( INDP+BN-1 ) = D( BN ) - LAMBDA
 354       DO 80 I = BN - 1, R1, -1
 355          DMINUS = LLD( I ) + WORK( INDP+I )
 356          TMP = D( I ) / DMINUS
 357          IF(DMINUS.LT.ZERO) NEG2 = NEG2 + 1
 358          WORK( INDUMN+I ) = L( I )*TMP
 359          WORK( INDP+I-1 ) = WORK( INDP+I )*TMP - LAMBDA
 360  80   CONTINUE
 361       TMP = WORK( INDP+R1-1 )
 362       SAWNAN2 = DISNAN( TMP )
 363
 364       IF( SAWNAN2 ) THEN
 365 *        Runs a slower version of the above loop if a NaN is detected
 366          NEG2 = 0
 367          DO 100 I = BN-1, R1, -1
 368             DMINUS = LLD( I ) + WORK( INDP+I )
 369             IF(ABS(DMINUS).LT.PIVMIN) DMINUS = -PIVMIN
 370             TMP = D( I ) / DMINUS
 371             IF(DMINUS.LT.ZERO) NEG2 = NEG2 + 1
 372             WORK( INDUMN+I ) = L( I )*TMP
 373             WORK( INDP+I-1 ) = WORK( INDP+I )*TMP - LAMBDA
 374             IF( TMP.EQ.ZERO )
 375      $          WORK( INDP+I-1 ) = D( I ) - LAMBDA
 376  100     CONTINUE
 377       END IF
 378 *
 379 *     Find the index (from R1 to R2) of the largest (in magnitude)
 380 *     diagonal element of the inverse
 381 *
 382       MINGMA = WORK( INDS+R1-1 ) + WORK( INDP+R1-1 )
 383       IF( MINGMA.LT.ZERO ) NEG1 = NEG1 + 1
 384       IF( WANTNC ) THEN
 385          NEGCNT = NEG1 + NEG2
 386       ELSE
 387          NEGCNT = -1
 388       ENDIF
 389       IF( ABS(MINGMA).EQ.ZERO )
 390      $   MINGMA = EPS*WORK( INDS+R1-1 )
 391       R = R1
 392       DO 110 I = R1, R2 - 1
 393          TMP = WORK( INDS+I ) + WORK( INDP+I )
 394          IF( TMP.EQ.ZERO )
 395      $      TMP = EPS*WORK( INDS+I )
 396          IF( ABS( TMP ).LE.ABS( MINGMA ) ) THEN
 397             MINGMA = TMP
 398             R = I + 1
 399          END IF
 400  110  CONTINUE
 401 *
 402 *     Compute the FP vector: solve N^T v = e_r
 403 *
 404       ISUPPZ( 1 ) = B1
 405       ISUPPZ( 2 ) = BN
 406       Z( R ) = ONE
 407       ZTZ = ONE
 408 *
 409 *     Compute the FP vector upwards from R
 410 *
 411       IF( .NOT.SAWNAN1 .AND. .NOT.SAWNAN2 ) THEN
 412          DO 210 I = R-1, B1, -1
 413             Z( I ) = -( WORK( INDLPL+I )*Z( I+1 ) )
 414             IF( (ABS(Z(I))+ABS(Z(I+1)))* ABS(LD(I)).LT.GAPTOL )
 415      $           THEN
 416                Z( I ) = ZERO
 417                ISUPPZ( 1 ) = I + 1
 418                GOTO 220
 419             ENDIF
 420             ZTZ = ZTZ + Z( I )*Z( I )
 421  210     CONTINUE
 422  220     CONTINUE
 423       ELSE
 424 *        Run slower loop if NaN occurred.
 425          DO 230 I = R - 1, B1, -1
 426             IF( Z( I+1 ).EQ.ZERO ) THEN
 427                Z( I ) = -( LD( I+1 ) / LD( I ) )*Z( I+2 )
 428             ELSE
 429                Z( I ) = -( WORK( INDLPL+I )*Z( I+1 ) )
 430             END IF
 431             IF( (ABS(Z(I))+ABS(Z(I+1)))* ABS(LD(I)).LT.GAPTOL )
 432      $           THEN
 433                Z( I ) = ZERO
 434                ISUPPZ( 1 ) = I + 1
 435                GO TO 240
 436             END IF
 437             ZTZ = ZTZ + Z( I )*Z( I )
 438  230     CONTINUE
 439  240     CONTINUE
 440       ENDIF
 441
 442 *     Compute the FP vector downwards from R in blocks of size BLKSIZ
 443       IF( .NOT.SAWNAN1 .AND. .NOT.SAWNAN2 ) THEN
 444          DO 250 I = R, BN-1
 445             Z( I+1 ) = -( WORK( INDUMN+I )*Z( I ) )
 446             IF( (ABS(Z(I))+ABS(Z(I+1)))* ABS(LD(I)).LT.GAPTOL )
 447      $         THEN
 448                Z( I+1 ) = ZERO
 449                ISUPPZ( 2 ) = I
 450                GO TO 260
 451             END IF
 452             ZTZ = ZTZ + Z( I+1 )*Z( I+1 )
 453  250     CONTINUE
 454  260     CONTINUE
 455       ELSE
 456 *        Run slower loop if NaN occurred.
 457          DO 270 I = R, BN - 1
 458             IF( Z( I ).EQ.ZERO ) THEN
 459                Z( I+1 ) = -( LD( I-1 ) / LD( I ) )*Z( I-1 )
 460             ELSE
 461                Z( I+1 ) = -( WORK( INDUMN+I )*Z( I ) )
 462             END IF
 463             IF( (ABS(Z(I))+ABS(Z(I+1)))* ABS(LD(I)).LT.GAPTOL )
 464      $           THEN
 465                Z( I+1 ) = ZERO
 466                ISUPPZ( 2 ) = I
 467                GO TO 280
 468             END IF
 469             ZTZ = ZTZ + Z( I+1 )*Z( I+1 )
 470  270     CONTINUE
 471  280     CONTINUE
 472       END IF
 473 *
 474 *     Compute quantities for convergence test
 475 *
 476       TMP = ONE / ZTZ
 477       NRMINV = SQRT( TMP )
 478       RESID = ABS( MINGMA )*NRMINV
 479       RQCORR = MINGMA*TMP
 480 *
 481 *
 482       RETURN
 483 *
 484 *     End of DLAR1V
 485 *
 486       END