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 December 2016 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.7.0) -- 234* -- LAPACK is a software package provided by Univ. of Tennessee, -- 235* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 236* December 2016 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 487