SRC/zlar1v.f

   1 *> \brief \b ZLAR1V 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 ZLAR1V + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zlar1v.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zlar1v.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zlar1v.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE ZLAR1V( 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 *       COMPLEX*16       Z( * )
  36 *       ..
  37 *
  38 *
  39 *> \par Purpose:
  40 *  =============
  41 *>
  42 *> \verbatim
  43 *>
  44 *> ZLAR1V 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 COMPLEX*16 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 complex16OTHERauxiliary
 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 ZLAR1V( 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       COMPLEX*16       Z( * )
 249 *     ..
 250 *
 251 *  =====================================================================
 252 *
 253 *     .. Parameters ..
 254       DOUBLE PRECISION   ZERO, ONE
 255       PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0 )
 256       COMPLEX*16         CONE
 257       PARAMETER          ( CONE = ( 1.0D0, 0.0D0 ) )
 258
 259 *     ..
 260 *     .. Local Scalars ..
 261       LOGICAL            SAWNAN1, SAWNAN2
 262       INTEGER            I, INDLPL, INDP, INDS, INDUMN, NEG1, NEG2, R1,
 263      $                   R2
 264       DOUBLE PRECISION   DMINUS, DPLUS, EPS, S, TMP
 265 *     ..
 266 *     .. External Functions ..
 267       LOGICAL DISNAN
 268       DOUBLE PRECISION   DLAMCH
 269       EXTERNAL           DISNAN, DLAMCH
 270 *     ..
 271 *     .. Intrinsic Functions ..
 272       INTRINSIC          ABS, DBLE
 273 *     ..
 274 *     .. Executable Statements ..
 275 *
 276       EPS = DLAMCH( 'Precision' )
 277
 278
 279       IF( R.EQ.0 ) THEN
 280          R1 = B1
 281          R2 = BN
 282       ELSE
 283          R1 = R
 284          R2 = R
 285       END IF
 286
 287 *     Storage for LPLUS
 288       INDLPL = 0
 289 *     Storage for UMINUS
 290       INDUMN = N
 291       INDS = 2*N + 1
 292       INDP = 3*N + 1
 293
 294       IF( B1.EQ.1 ) THEN
 295          WORK( INDS ) = ZERO
 296       ELSE
 297          WORK( INDS+B1-1 ) = LLD( B1-1 )
 298       END IF
 299
 300 *
 301 *     Compute the stationary transform (using the differential form)
 302 *     until the index R2.
 303 *
 304       SAWNAN1 = .FALSE.
 305       NEG1 = 0
 306       S = WORK( INDS+B1-1 ) - LAMBDA
 307       DO 50 I = B1, R1 - 1
 308          DPLUS = D( I ) + S
 309          WORK( INDLPL+I ) = LD( I ) / DPLUS
 310          IF(DPLUS.LT.ZERO) NEG1 = NEG1 + 1
 311          WORK( INDS+I ) = S*WORK( INDLPL+I )*L( I )
 312          S = WORK( INDS+I ) - LAMBDA
 313  50   CONTINUE
 314       SAWNAN1 = DISNAN( S )
 315       IF( SAWNAN1 ) GOTO 60
 316       DO 51 I = R1, R2 - 1
 317          DPLUS = D( I ) + S
 318          WORK( INDLPL+I ) = LD( I ) / DPLUS
 319          WORK( INDS+I ) = S*WORK( INDLPL+I )*L( I )
 320          S = WORK( INDS+I ) - LAMBDA
 321  51   CONTINUE
 322       SAWNAN1 = DISNAN( S )
 323 *
 324  60   CONTINUE
 325       IF( SAWNAN1 ) THEN
 326 *        Runs a slower version of the above loop if a NaN is detected
 327          NEG1 = 0
 328          S = WORK( INDS+B1-1 ) - LAMBDA
 329          DO 70 I = B1, R1 - 1
 330             DPLUS = D( I ) + S
 331             IF(ABS(DPLUS).LT.PIVMIN) DPLUS = -PIVMIN
 332             WORK( INDLPL+I ) = LD( I ) / DPLUS
 333             IF(DPLUS.LT.ZERO) NEG1 = NEG1 + 1
 334             WORK( INDS+I ) = S*WORK( INDLPL+I )*L( I )
 335             IF( WORK( INDLPL+I ).EQ.ZERO )
 336      $                      WORK( INDS+I ) = LLD( I )
 337             S = WORK( INDS+I ) - LAMBDA
 338  70      CONTINUE
 339          DO 71 I = R1, R2 - 1
 340             DPLUS = D( I ) + S
 341             IF(ABS(DPLUS).LT.PIVMIN) DPLUS = -PIVMIN
 342             WORK( INDLPL+I ) = LD( I ) / DPLUS
 343             WORK( INDS+I ) = S*WORK( INDLPL+I )*L( I )
 344             IF( WORK( INDLPL+I ).EQ.ZERO )
 345      $                      WORK( INDS+I ) = LLD( I )
 346             S = WORK( INDS+I ) - LAMBDA
 347  71      CONTINUE
 348       END IF
 349 *
 350 *     Compute the progressive transform (using the differential form)
 351 *     until the index R1
 352 *
 353       SAWNAN2 = .FALSE.
 354       NEG2 = 0
 355       WORK( INDP+BN-1 ) = D( BN ) - LAMBDA
 356       DO 80 I = BN - 1, R1, -1
 357          DMINUS = LLD( I ) + WORK( INDP+I )
 358          TMP = D( I ) / DMINUS
 359          IF(DMINUS.LT.ZERO) NEG2 = NEG2 + 1
 360          WORK( INDUMN+I ) = L( I )*TMP
 361          WORK( INDP+I-1 ) = WORK( INDP+I )*TMP - LAMBDA
 362  80   CONTINUE
 363       TMP = WORK( INDP+R1-1 )
 364       SAWNAN2 = DISNAN( TMP )
 365
 366       IF( SAWNAN2 ) THEN
 367 *        Runs a slower version of the above loop if a NaN is detected
 368          NEG2 = 0
 369          DO 100 I = BN-1, R1, -1
 370             DMINUS = LLD( I ) + WORK( INDP+I )
 371             IF(ABS(DMINUS).LT.PIVMIN) DMINUS = -PIVMIN
 372             TMP = D( I ) / DMINUS
 373             IF(DMINUS.LT.ZERO) NEG2 = NEG2 + 1
 374             WORK( INDUMN+I ) = L( I )*TMP
 375             WORK( INDP+I-1 ) = WORK( INDP+I )*TMP - LAMBDA
 376             IF( TMP.EQ.ZERO )
 377      $          WORK( INDP+I-1 ) = D( I ) - LAMBDA
 378  100     CONTINUE
 379       END IF
 380 *
 381 *     Find the index (from R1 to R2) of the largest (in magnitude)
 382 *     diagonal element of the inverse
 383 *
 384       MINGMA = WORK( INDS+R1-1 ) + WORK( INDP+R1-1 )
 385       IF( MINGMA.LT.ZERO ) NEG1 = NEG1 + 1
 386       IF( WANTNC ) THEN
 387          NEGCNT = NEG1 + NEG2
 388       ELSE
 389          NEGCNT = -1
 390       ENDIF
 391       IF( ABS(MINGMA).EQ.ZERO )
 392      $   MINGMA = EPS*WORK( INDS+R1-1 )
 393       R = R1
 394       DO 110 I = R1, R2 - 1
 395          TMP = WORK( INDS+I ) + WORK( INDP+I )
 396          IF( TMP.EQ.ZERO )
 397      $      TMP = EPS*WORK( INDS+I )
 398          IF( ABS( TMP ).LE.ABS( MINGMA ) ) THEN
 399             MINGMA = TMP
 400             R = I + 1
 401          END IF
 402  110  CONTINUE
 403 *
 404 *     Compute the FP vector: solve N^T v = e_r
 405 *
 406       ISUPPZ( 1 ) = B1
 407       ISUPPZ( 2 ) = BN
 408       Z( R ) = CONE
 409       ZTZ = ONE
 410 *
 411 *     Compute the FP vector upwards from R
 412 *
 413       IF( .NOT.SAWNAN1 .AND. .NOT.SAWNAN2 ) THEN
 414          DO 210 I = R-1, B1, -1
 415             Z( I ) = -( WORK( INDLPL+I )*Z( I+1 ) )
 416             IF( (ABS(Z(I))+ABS(Z(I+1)))* ABS(LD(I)).LT.GAPTOL )
 417      $           THEN
 418                Z( I ) = ZERO
 419                ISUPPZ( 1 ) = I + 1
 420                GOTO 220
 421             ENDIF
 422             ZTZ = ZTZ + DBLE( Z( I )*Z( I ) )
 423  210     CONTINUE
 424  220     CONTINUE
 425       ELSE
 426 *        Run slower loop if NaN occurred.
 427          DO 230 I = R - 1, B1, -1
 428             IF( Z( I+1 ).EQ.ZERO ) THEN
 429                Z( I ) = -( LD( I+1 ) / LD( I ) )*Z( I+2 )
 430             ELSE
 431                Z( I ) = -( WORK( INDLPL+I )*Z( I+1 ) )
 432             END IF
 433             IF( (ABS(Z(I))+ABS(Z(I+1)))* ABS(LD(I)).LT.GAPTOL )
 434      $           THEN
 435                Z( I ) = ZERO
 436                ISUPPZ( 1 ) = I + 1
 437                GO TO 240
 438             END IF
 439             ZTZ = ZTZ + DBLE( Z( I )*Z( I ) )
 440  230     CONTINUE
 441  240     CONTINUE
 442       ENDIF
 443
 444 *     Compute the FP vector downwards from R in blocks of size BLKSIZ
 445       IF( .NOT.SAWNAN1 .AND. .NOT.SAWNAN2 ) THEN
 446          DO 250 I = R, BN-1
 447             Z( I+1 ) = -( WORK( INDUMN+I )*Z( I ) )
 448             IF( (ABS(Z(I))+ABS(Z(I+1)))* ABS(LD(I)).LT.GAPTOL )
 449      $         THEN
 450                Z( I+1 ) = ZERO
 451                ISUPPZ( 2 ) = I
 452                GO TO 260
 453             END IF
 454             ZTZ = ZTZ + DBLE( Z( I+1 )*Z( I+1 ) )
 455  250     CONTINUE
 456  260     CONTINUE
 457       ELSE
 458 *        Run slower loop if NaN occurred.
 459          DO 270 I = R, BN - 1
 460             IF( Z( I ).EQ.ZERO ) THEN
 461                Z( I+1 ) = -( LD( I-1 ) / LD( I ) )*Z( I-1 )
 462             ELSE
 463                Z( I+1 ) = -( WORK( INDUMN+I )*Z( I ) )
 464             END IF
 465             IF( (ABS(Z(I))+ABS(Z(I+1)))* ABS(LD(I)).LT.GAPTOL )
 466      $           THEN
 467                Z( I+1 ) = ZERO
 468                ISUPPZ( 2 ) = I
 469                GO TO 280
 470             END IF
 471             ZTZ = ZTZ + DBLE( Z( I+1 )*Z( I+1 ) )
 472  270     CONTINUE
 473  280     CONTINUE
 474       END IF
 475 *
 476 *     Compute quantities for convergence test
 477 *
 478       TMP = ONE / ZTZ
 479       NRMINV = SQRT( TMP )
 480       RESID = ABS( MINGMA )*NRMINV
 481       RQCORR = MINGMA*TMP
 482 *
 483 *
 484       RETURN
 485 *
 486 *     End of ZLAR1V
 487 *
 488       END