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*> \ingroup doubleOTHERauxiliary 216* 217*> \par Contributors: 218* ================== 219*> 220*> Beresford Parlett, University of California, Berkeley, USA \n 221*> Jim Demmel, University of California, Berkeley, USA \n 222*> Inderjit Dhillon, University of Texas, Austin, USA \n 223*> Osni Marques, LBNL/NERSC, USA \n 224*> Christof Voemel, University of California, Berkeley, USA 225* 226* ===================================================================== 227 SUBROUTINE DLAR1V( N, B1, BN, LAMBDA, D, L, LD, LLD, 228 $ PIVMIN, GAPTOL, Z, WANTNC, NEGCNT, ZTZ, MINGMA, 229 $ R, ISUPPZ, NRMINV, RESID, RQCORR, WORK ) 230* 231* -- LAPACK auxiliary routine -- 232* -- LAPACK is a software package provided by Univ. of Tennessee, -- 233* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 234* 235* .. Scalar Arguments .. 236 LOGICAL WANTNC 237 INTEGER B1, BN, N, NEGCNT, R 238 DOUBLE PRECISION GAPTOL, LAMBDA, MINGMA, NRMINV, PIVMIN, RESID, 239 $ RQCORR, ZTZ 240* .. 241* .. Array Arguments .. 242 INTEGER ISUPPZ( * ) 243 DOUBLE PRECISION D( * ), L( * ), LD( * ), LLD( * ), 244 $ WORK( * ) 245 DOUBLE PRECISION Z( * ) 246* .. 247* 248* ===================================================================== 249* 250* .. Parameters .. 251 DOUBLE PRECISION ZERO, ONE 252 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 ) 253 254* .. 255* .. Local Scalars .. 256 LOGICAL SAWNAN1, SAWNAN2 257 INTEGER I, INDLPL, INDP, INDS, INDUMN, NEG1, NEG2, R1, 258 $ R2 259 DOUBLE PRECISION DMINUS, DPLUS, EPS, S, TMP 260* .. 261* .. External Functions .. 262 LOGICAL DISNAN 263 DOUBLE PRECISION DLAMCH 264 EXTERNAL DISNAN, DLAMCH 265* .. 266* .. Intrinsic Functions .. 267 INTRINSIC ABS 268* .. 269* .. Executable Statements .. 270* 271 EPS = DLAMCH( 'Precision' ) 272 273 274 IF( R.EQ.0 ) THEN 275 R1 = B1 276 R2 = BN 277 ELSE 278 R1 = R 279 R2 = R 280 END IF 281 282* Storage for LPLUS 283 INDLPL = 0 284* Storage for UMINUS 285 INDUMN = N 286 INDS = 2*N + 1 287 INDP = 3*N + 1 288 289 IF( B1.EQ.1 ) THEN 290 WORK( INDS ) = ZERO 291 ELSE 292 WORK( INDS+B1-1 ) = LLD( B1-1 ) 293 END IF 294 295* 296* Compute the stationary transform (using the differential form) 297* until the index R2. 298* 299 SAWNAN1 = .FALSE. 300 NEG1 = 0 301 S = WORK( INDS+B1-1 ) - LAMBDA 302 DO 50 I = B1, R1 - 1 303 DPLUS = D( I ) + S 304 WORK( INDLPL+I ) = LD( I ) / DPLUS 305 IF(DPLUS.LT.ZERO) NEG1 = NEG1 + 1 306 WORK( INDS+I ) = S*WORK( INDLPL+I )*L( I ) 307 S = WORK( INDS+I ) - LAMBDA 308 50 CONTINUE 309 SAWNAN1 = DISNAN( S ) 310 IF( SAWNAN1 ) GOTO 60 311 DO 51 I = R1, R2 - 1 312 DPLUS = D( I ) + S 313 WORK( INDLPL+I ) = LD( I ) / DPLUS 314 WORK( INDS+I ) = S*WORK( INDLPL+I )*L( I ) 315 S = WORK( INDS+I ) - LAMBDA 316 51 CONTINUE 317 SAWNAN1 = DISNAN( S ) 318* 319 60 CONTINUE 320 IF( SAWNAN1 ) THEN 321* Runs a slower version of the above loop if a NaN is detected 322 NEG1 = 0 323 S = WORK( INDS+B1-1 ) - LAMBDA 324 DO 70 I = B1, R1 - 1 325 DPLUS = D( I ) + S 326 IF(ABS(DPLUS).LT.PIVMIN) DPLUS = -PIVMIN 327 WORK( INDLPL+I ) = LD( I ) / DPLUS 328 IF(DPLUS.LT.ZERO) NEG1 = NEG1 + 1 329 WORK( INDS+I ) = S*WORK( INDLPL+I )*L( I ) 330 IF( WORK( INDLPL+I ).EQ.ZERO ) 331 $ WORK( INDS+I ) = LLD( I ) 332 S = WORK( INDS+I ) - LAMBDA 333 70 CONTINUE 334 DO 71 I = R1, R2 - 1 335 DPLUS = D( I ) + S 336 IF(ABS(DPLUS).LT.PIVMIN) DPLUS = -PIVMIN 337 WORK( INDLPL+I ) = LD( I ) / DPLUS 338 WORK( INDS+I ) = S*WORK( INDLPL+I )*L( I ) 339 IF( WORK( INDLPL+I ).EQ.ZERO ) 340 $ WORK( INDS+I ) = LLD( I ) 341 S = WORK( INDS+I ) - LAMBDA 342 71 CONTINUE 343 END IF 344* 345* Compute the progressive transform (using the differential form) 346* until the index R1 347* 348 SAWNAN2 = .FALSE. 349 NEG2 = 0 350 WORK( INDP+BN-1 ) = D( BN ) - LAMBDA 351 DO 80 I = BN - 1, R1, -1 352 DMINUS = LLD( I ) + WORK( INDP+I ) 353 TMP = D( I ) / DMINUS 354 IF(DMINUS.LT.ZERO) NEG2 = NEG2 + 1 355 WORK( INDUMN+I ) = L( I )*TMP 356 WORK( INDP+I-1 ) = WORK( INDP+I )*TMP - LAMBDA 357 80 CONTINUE 358 TMP = WORK( INDP+R1-1 ) 359 SAWNAN2 = DISNAN( TMP ) 360 361 IF( SAWNAN2 ) THEN 362* Runs a slower version of the above loop if a NaN is detected 363 NEG2 = 0 364 DO 100 I = BN-1, R1, -1 365 DMINUS = LLD( I ) + WORK( INDP+I ) 366 IF(ABS(DMINUS).LT.PIVMIN) DMINUS = -PIVMIN 367 TMP = D( I ) / DMINUS 368 IF(DMINUS.LT.ZERO) NEG2 = NEG2 + 1 369 WORK( INDUMN+I ) = L( I )*TMP 370 WORK( INDP+I-1 ) = WORK( INDP+I )*TMP - LAMBDA 371 IF( TMP.EQ.ZERO ) 372 $ WORK( INDP+I-1 ) = D( I ) - LAMBDA 373 100 CONTINUE 374 END IF 375* 376* Find the index (from R1 to R2) of the largest (in magnitude) 377* diagonal element of the inverse 378* 379 MINGMA = WORK( INDS+R1-1 ) + WORK( INDP+R1-1 ) 380 IF( MINGMA.LT.ZERO ) NEG1 = NEG1 + 1 381 IF( WANTNC ) THEN 382 NEGCNT = NEG1 + NEG2 383 ELSE 384 NEGCNT = -1 385 ENDIF 386 IF( ABS(MINGMA).EQ.ZERO ) 387 $ MINGMA = EPS*WORK( INDS+R1-1 ) 388 R = R1 389 DO 110 I = R1, R2 - 1 390 TMP = WORK( INDS+I ) + WORK( INDP+I ) 391 IF( TMP.EQ.ZERO ) 392 $ TMP = EPS*WORK( INDS+I ) 393 IF( ABS( TMP ).LE.ABS( MINGMA ) ) THEN 394 MINGMA = TMP 395 R = I + 1 396 END IF 397 110 CONTINUE 398* 399* Compute the FP vector: solve N^T v = e_r 400* 401 ISUPPZ( 1 ) = B1 402 ISUPPZ( 2 ) = BN 403 Z( R ) = ONE 404 ZTZ = ONE 405* 406* Compute the FP vector upwards from R 407* 408 IF( .NOT.SAWNAN1 .AND. .NOT.SAWNAN2 ) THEN 409 DO 210 I = R-1, B1, -1 410 Z( I ) = -( WORK( INDLPL+I )*Z( I+1 ) ) 411 IF( (ABS(Z(I))+ABS(Z(I+1)))* ABS(LD(I)).LT.GAPTOL ) 412 $ THEN 413 Z( I ) = ZERO 414 ISUPPZ( 1 ) = I + 1 415 GOTO 220 416 ENDIF 417 ZTZ = ZTZ + Z( I )*Z( I ) 418 210 CONTINUE 419 220 CONTINUE 420 ELSE 421* Run slower loop if NaN occurred. 422 DO 230 I = R - 1, B1, -1 423 IF( Z( I+1 ).EQ.ZERO ) THEN 424 Z( I ) = -( LD( I+1 ) / LD( I ) )*Z( I+2 ) 425 ELSE 426 Z( I ) = -( WORK( INDLPL+I )*Z( I+1 ) ) 427 END IF 428 IF( (ABS(Z(I))+ABS(Z(I+1)))* ABS(LD(I)).LT.GAPTOL ) 429 $ THEN 430 Z( I ) = ZERO 431 ISUPPZ( 1 ) = I + 1 432 GO TO 240 433 END IF 434 ZTZ = ZTZ + Z( I )*Z( I ) 435 230 CONTINUE 436 240 CONTINUE 437 ENDIF 438 439* Compute the FP vector downwards from R in blocks of size BLKSIZ 440 IF( .NOT.SAWNAN1 .AND. .NOT.SAWNAN2 ) THEN 441 DO 250 I = R, BN-1 442 Z( I+1 ) = -( WORK( INDUMN+I )*Z( I ) ) 443 IF( (ABS(Z(I))+ABS(Z(I+1)))* ABS(LD(I)).LT.GAPTOL ) 444 $ THEN 445 Z( I+1 ) = ZERO 446 ISUPPZ( 2 ) = I 447 GO TO 260 448 END IF 449 ZTZ = ZTZ + Z( I+1 )*Z( I+1 ) 450 250 CONTINUE 451 260 CONTINUE 452 ELSE 453* Run slower loop if NaN occurred. 454 DO 270 I = R, BN - 1 455 IF( Z( I ).EQ.ZERO ) THEN 456 Z( I+1 ) = -( LD( I-1 ) / LD( I ) )*Z( I-1 ) 457 ELSE 458 Z( I+1 ) = -( WORK( INDUMN+I )*Z( I ) ) 459 END IF 460 IF( (ABS(Z(I))+ABS(Z(I+1)))* ABS(LD(I)).LT.GAPTOL ) 461 $ THEN 462 Z( I+1 ) = ZERO 463 ISUPPZ( 2 ) = I 464 GO TO 280 465 END IF 466 ZTZ = ZTZ + Z( I+1 )*Z( I+1 ) 467 270 CONTINUE 468 280 CONTINUE 469 END IF 470* 471* Compute quantities for convergence test 472* 473 TMP = ONE / ZTZ 474 NRMINV = SQRT( TMP ) 475 RESID = ABS( MINGMA )*NRMINV 476 RQCORR = MINGMA*TMP 477* 478* 479 RETURN 480* 481* End of DLAR1V 482* 483 END 484