1 SUBROUTINE DLAR1V( N, B1, BN, SIGMA, D, L, LD, LLD, GERSCH, Z, 2 $ ZTZ, MINGMA, R, ISUPPZ, WORK ) 3* 4* -- LAPACK auxiliary routine (version 3.0) -- 5* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., 6* Courant Institute, Argonne National Lab, and Rice University 7* June 30, 1999 8* 9* .. Scalar Arguments .. 10 INTEGER B1, BN, N, R 11 DOUBLE PRECISION MINGMA, SIGMA, ZTZ 12* .. 13* .. Array Arguments .. 14 INTEGER ISUPPZ( * ) 15 DOUBLE PRECISION D( * ), GERSCH( * ), L( * ), LD( * ), LLD( * ), 16 $ WORK( * ), Z( * ) 17* .. 18* 19* Purpose 20* ======= 21* 22* DLAR1V computes the (scaled) r-th column of the inverse of 23* the sumbmatrix in rows B1 through BN of the tridiagonal matrix 24* L D L^T - sigma I. The following steps accomplish this computation : 25* (a) Stationary qd transform, L D L^T - sigma I = L(+) D(+) L(+)^T, 26* (b) Progressive qd transform, L D L^T - sigma I = U(-) D(-) U(-)^T, 27* (c) Computation of the diagonal elements of the inverse of 28* L D L^T - sigma I by combining the above transforms, and choosing 29* r as the index where the diagonal of the inverse is (one of the) 30* largest in magnitude. 31* (d) Computation of the (scaled) r-th column of the inverse using the 32* twisted factorization obtained by combining the top part of the 33* the stationary and the bottom part of the progressive transform. 34* 35* Arguments 36* ========= 37* 38* N (input) INTEGER 39* The order of the matrix L D L^T. 40* 41* B1 (input) INTEGER 42* First index of the submatrix of L D L^T. 43* 44* BN (input) INTEGER 45* Last index of the submatrix of L D L^T. 46* 47* SIGMA (input) DOUBLE PRECISION 48* The shift. Initially, when R = 0, SIGMA should be a good 49* approximation to an eigenvalue of L D L^T. 50* 51* L (input) DOUBLE PRECISION array, dimension (N-1) 52* The (n-1) subdiagonal elements of the unit bidiagonal matrix 53* L, in elements 1 to N-1. 54* 55* D (input) DOUBLE PRECISION array, dimension (N) 56* The n diagonal elements of the diagonal matrix D. 57* 58* LD (input) DOUBLE PRECISION array, dimension (N-1) 59* The n-1 elements L(i)*D(i). 60* 61* LLD (input) DOUBLE PRECISION array, dimension (N-1) 62* The n-1 elements L(i)*L(i)*D(i). 63* 64* GERSCH (input) DOUBLE PRECISION array, dimension (2*N) 65* The n Gerschgorin intervals. These are used to restrict 66* the initial search for R, when R is input as 0. 67* 68* Z (output) DOUBLE PRECISION array, dimension (N) 69* The (scaled) r-th column of the inverse. Z(R) is returned 70* to be 1. 71* 72* ZTZ (output) DOUBLE PRECISION 73* The square of the norm of Z. 74* 75* MINGMA (output) DOUBLE PRECISION 76* The reciprocal of the largest (in magnitude) diagonal 77* element of the inverse of L D L^T - sigma I. 78* 79* R (input/output) INTEGER 80* Initially, R should be input to be 0 and is then output as 81* the index where the diagonal element of the inverse is 82* largest in magnitude. In later iterations, this same value 83* of R should be input. 84* 85* ISUPPZ (output) INTEGER array, dimension (2) 86* The support of the vector in Z, i.e., the vector Z is 87* nonzero only in elements ISUPPZ(1) through ISUPPZ( 2 ). 88* 89* WORK (workspace) DOUBLE PRECISION array, dimension (4*N) 90* 91* Further Details 92* =============== 93* 94* Based on contributions by 95* Inderjit Dhillon, IBM Almaden, USA 96* Osni Marques, LBNL/NERSC, USA 97* 98* ===================================================================== 99* 100* .. Parameters .. 101 INTEGER BLKSIZ 102 PARAMETER ( BLKSIZ = 32 ) 103 DOUBLE PRECISION ZERO, ONE 104 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 ) 105* .. 106* .. Local Scalars .. 107 LOGICAL SAWNAN 108 INTEGER FROM, I, INDP, INDS, INDUMN, J, R1, R2, TO 109 DOUBLE PRECISION DMINUS, DPLUS, EPS, S, TMP 110* .. 111* .. External Functions .. 112 DOUBLE PRECISION DLAMCH 113 EXTERNAL DLAMCH 114* .. 115* .. Intrinsic Functions .. 116 INTRINSIC ABS, MAX, MIN 117* .. 118* .. Executable Statements .. 119* 120 EPS = DLAMCH( 'Precision' ) 121 IF( R.EQ.0 ) THEN 122* 123* Eliminate the top and bottom indices from the possible values 124* of R where the desired eigenvector is largest in magnitude. 125* 126 R1 = B1 127 DO 10 I = B1, BN 128 IF( SIGMA.GE.GERSCH( 2*I-1 ) .OR. SIGMA.LE.GERSCH( 2*I ) ) 129 $ THEN 130 R1 = I 131 GO TO 20 132 END IF 133 10 CONTINUE 134 20 CONTINUE 135 R2 = BN 136 DO 30 I = BN, B1, -1 137 IF( SIGMA.GE.GERSCH( 2*I-1 ) .OR. SIGMA.LE.GERSCH( 2*I ) ) 138 $ THEN 139 R2 = I 140 GO TO 40 141 END IF 142 30 CONTINUE 143 40 CONTINUE 144 ELSE 145 R1 = R 146 R2 = R 147 END IF 148* 149 INDUMN = N 150 INDS = 2*N + 1 151 INDP = 3*N + 1 152 SAWNAN = .FALSE. 153* 154* Compute the stationary transform (using the differential form) 155* untill the index R2 156* 157 IF( B1.EQ.1 ) THEN 158 WORK( INDS ) = ZERO 159 ELSE 160 WORK( INDS ) = LLD( B1-1 ) 161 END IF 162 S = WORK( INDS ) - SIGMA 163 DO 50 I = B1, R2 - 1 164 DPLUS = D( I ) + S 165 WORK( I ) = LD( I ) / DPLUS 166 WORK( INDS+I ) = S*WORK( I )*L( I ) 167 S = WORK( INDS+I ) - SIGMA 168 50 CONTINUE 169* 170 IF( .NOT.( S.GT.ZERO .OR. S.LT.ONE ) ) THEN 171* 172* Run a slower version of the above loop if a NaN is detected 173* 174 SAWNAN = .TRUE. 175 J = B1 + 1 176 60 CONTINUE 177 IF( WORK( INDS+J ).GT.ZERO .OR. WORK( INDS+J ).LT.ONE ) THEN 178 J = J + 1 179 GO TO 60 180 END IF 181 WORK( INDS+J ) = LLD( J ) 182 S = WORK( INDS+J ) - SIGMA 183 DO 70 I = J + 1, R2 - 1 184 DPLUS = D( I ) + S 185 WORK( I ) = LD( I ) / DPLUS 186 IF( WORK( I ).EQ.ZERO ) THEN 187 WORK( INDS+I ) = LLD( I ) 188 ELSE 189 WORK( INDS+I ) = S*WORK( I )*L( I ) 190 END IF 191 S = WORK( INDS+I ) - SIGMA 192 70 CONTINUE 193 END IF 194 WORK( INDP+BN-1 ) = D( BN ) - SIGMA 195 DO 80 I = BN - 1, R1, -1 196 DMINUS = LLD( I ) + WORK( INDP+I ) 197 TMP = D( I ) / DMINUS 198 WORK( INDUMN+I ) = L( I )*TMP 199 WORK( INDP+I-1 ) = WORK( INDP+I )*TMP - SIGMA 200 80 CONTINUE 201 TMP = WORK( INDP+R1-1 ) 202 IF( .NOT.( TMP.GT.ZERO .OR. TMP.LT.ONE ) ) THEN 203* 204* Run a slower version of the above loop if a NaN is detected 205* 206 SAWNAN = .TRUE. 207 J = BN - 3 208 90 CONTINUE 209 IF( WORK( INDP+J ).GT.ZERO .OR. WORK( INDP+J ).LT.ONE ) THEN 210 J = J - 1 211 GO TO 90 212 END IF 213 WORK( INDP+J ) = D( J+1 ) - SIGMA 214 DO 100 I = J, R1, -1 215 DMINUS = LLD( I ) + WORK( INDP+I ) 216 TMP = D( I ) / DMINUS 217 WORK( INDUMN+I ) = L( I )*TMP 218 IF( TMP.EQ.ZERO ) THEN 219 WORK( INDP+I-1 ) = D( I ) - SIGMA 220 ELSE 221 WORK( INDP+I-1 ) = WORK( INDP+I )*TMP - SIGMA 222 END IF 223 100 CONTINUE 224 END IF 225* 226* Find the index (from R1 to R2) of the largest (in magnitude) 227* diagonal element of the inverse 228* 229 MINGMA = WORK( INDS+R1-1 ) + WORK( INDP+R1-1 ) 230 IF( MINGMA.EQ.ZERO ) 231 $ MINGMA = EPS*WORK( INDS+R1-1 ) 232 R = R1 233 DO 110 I = R1, R2 - 1 234 TMP = WORK( INDS+I ) + WORK( INDP+I ) 235 IF( TMP.EQ.ZERO ) 236 $ TMP = EPS*WORK( INDS+I ) 237 IF( ABS( TMP ).LT.ABS( MINGMA ) ) THEN 238 MINGMA = TMP 239 R = I + 1 240 END IF 241 110 CONTINUE 242* 243* Compute the (scaled) r-th column of the inverse 244* 245 ISUPPZ( 1 ) = B1 246 ISUPPZ( 2 ) = BN 247 Z( R ) = ONE 248 ZTZ = ONE 249 IF( .NOT.SAWNAN ) THEN 250 FROM = R - 1 251 TO = MAX( R-BLKSIZ, B1 ) 252 120 CONTINUE 253 IF( FROM.GE.B1 ) THEN 254 DO 130 I = FROM, TO, -1 255 Z( I ) = -( WORK( I )*Z( I+1 ) ) 256 ZTZ = ZTZ + Z( I )*Z( I ) 257 130 CONTINUE 258 IF( ABS( Z( TO ) ).LE.EPS .AND. ABS( Z( TO+1 ) ).LE.EPS ) 259 $ THEN 260 ISUPPZ( 1 ) = TO + 2 261 ELSE 262 FROM = TO - 1 263 TO = MAX( TO-BLKSIZ, B1 ) 264 GO TO 120 265 END IF 266 END IF 267 FROM = R + 1 268 TO = MIN( R+BLKSIZ, BN ) 269 140 CONTINUE 270 IF( FROM.LE.BN ) THEN 271 DO 150 I = FROM, TO 272 Z( I ) = -( WORK( INDUMN+I-1 )*Z( I-1 ) ) 273 ZTZ = ZTZ + Z( I )*Z( I ) 274 150 CONTINUE 275 IF( ABS( Z( TO ) ).LE.EPS .AND. ABS( Z( TO-1 ) ).LE.EPS ) 276 $ THEN 277 ISUPPZ( 2 ) = TO - 2 278 ELSE 279 FROM = TO + 1 280 TO = MIN( TO+BLKSIZ, BN ) 281 GO TO 140 282 END IF 283 END IF 284 ELSE 285 DO 160 I = R - 1, B1, -1 286 IF( Z( I+1 ).EQ.ZERO ) THEN 287 Z( I ) = -( LD( I+1 ) / LD( I ) )*Z( I+2 ) 288 ELSE IF( ABS( Z( I+1 ) ).LE.EPS .AND. ABS( Z( I+2 ) ).LE. 289 $ EPS ) THEN 290 ISUPPZ( 1 ) = I + 3 291 GO TO 170 292 ELSE 293 Z( I ) = -( WORK( I )*Z( I+1 ) ) 294 END IF 295 ZTZ = ZTZ + Z( I )*Z( I ) 296 160 CONTINUE 297 170 CONTINUE 298 DO 180 I = R, BN - 1 299 IF( Z( I ).EQ.ZERO ) THEN 300 Z( I+1 ) = -( LD( I-1 ) / LD( I ) )*Z( I-1 ) 301 ELSE IF( ABS( Z( I ) ).LE.EPS .AND. ABS( Z( I-1 ) ).LE.EPS ) 302 $ THEN 303 ISUPPZ( 2 ) = I - 2 304 GO TO 190 305 ELSE 306 Z( I+1 ) = -( WORK( INDUMN+I )*Z( I ) ) 307 END IF 308 ZTZ = ZTZ + Z( I+1 )*Z( I+1 ) 309 180 CONTINUE 310 190 CONTINUE 311 END IF 312 DO 200 I = B1, ISUPPZ( 1 ) - 3 313 Z( I ) = ZERO 314 200 CONTINUE 315 DO 210 I = ISUPPZ( 2 ) + 3, BN 316 Z( I ) = ZERO 317 210 CONTINUE 318* 319 RETURN 320* 321* End of DLAR1V 322* 323 END 324