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