1 SUBROUTINE DLARRV( N, D, L, ISPLIT, M, W, IBLOCK, GERSCH, TOL, Z, 2 $ LDZ, ISUPPZ, WORK, IWORK, INFO ) 3* 4* -- LAPACK auxiliary routine (instru to count ops, 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 INFO, LDZ, M, N 11 DOUBLE PRECISION TOL 12* .. 13* .. Array Arguments .. 14 INTEGER IBLOCK( * ), ISPLIT( * ), ISUPPZ( * ), 15 $ IWORK( * ) 16 DOUBLE PRECISION D( * ), GERSCH( * ), L( * ), W( * ), WORK( * ), 17 $ Z( LDZ, * ) 18* .. 19* Common block to return operation count 20* .. Common blocks .. 21 COMMON / LATIME / OPS, ITCNT 22* .. 23* .. Scalars in Common .. 24 DOUBLE PRECISION ITCNT, OPS 25* .. 26* 27* Purpose 28* ======= 29* 30* DLARRV computes the eigenvectors of the tridiagonal matrix 31* T = L D L^T given L, D and the eigenvalues of L D L^T. 32* The input eigenvalues should have high relative accuracy with 33* respect to the entries of L and D. The desired accuracy of the 34* output can be specified by the input parameter TOL. 35* 36* Arguments 37* ========= 38* 39* N (input) INTEGER 40* The order of the matrix. N >= 0. 41* 42* D (input/output) DOUBLE PRECISION array, dimension (N) 43* On entry, the n diagonal elements of the diagonal matrix D. 44* On exit, D may be overwritten. 45* 46* L (input/output) DOUBLE PRECISION array, dimension (N-1) 47* On entry, the (n-1) subdiagonal elements of the unit 48* bidiagonal matrix L in elements 1 to N-1 of L. L(N) need 49* not be set. On exit, L is overwritten. 50* 51* ISPLIT (input) INTEGER array, dimension (N) 52* The splitting points, at which T breaks up into submatrices. 53* The first submatrix consists of rows/columns 1 to 54* ISPLIT( 1 ), the second of rows/columns ISPLIT( 1 )+1 55* through ISPLIT( 2 ), etc. 56* 57* TOL (input) DOUBLE PRECISION 58* The absolute error tolerance for the 59* eigenvalues/eigenvectors. 60* Errors in the input eigenvalues must be bounded by TOL. 61* The eigenvectors output have residual norms 62* bounded by TOL, and the dot products between different 63* eigenvectors are bounded by TOL. TOL must be at least 64* N*EPS*|T|, where EPS is the machine precision and |T| is 65* the 1-norm of the tridiagonal matrix. 66* 67* M (input) INTEGER 68* The total number of eigenvalues found. 0 <= M <= N. 69* If RANGE = 'A', M = N, and if RANGE = 'I', M = IU-IL+1. 70* 71* W (input) DOUBLE PRECISION array, dimension (N) 72* The first M elements of W contain the eigenvalues for 73* which eigenvectors are to be computed. The eigenvalues 74* should be grouped by split-off block and ordered from 75* smallest to largest within the block ( The output array 76* W from DLARRE is expected here ). 77* Errors in W must be bounded by TOL (see above). 78* 79* IBLOCK (input) INTEGER array, dimension (N) 80* The submatrix indices associated with the corresponding 81* eigenvalues in W; IBLOCK(i)=1 if eigenvalue W(i) belongs to 82* the first submatrix from the top, =2 if W(i) belongs to 83* the second submatrix, etc. 84* 85* Z (output) DOUBLE PRECISION array, dimension (LDZ, max(1,M) ) 86* If JOBZ = 'V', then if INFO = 0, the first M columns of Z 87* contain the orthonormal eigenvectors of the matrix T 88* corresponding to the selected eigenvalues, with the i-th 89* column of Z holding the eigenvector associated with W(i). 90* If JOBZ = 'N', then Z is not referenced. 91* Note: the user must ensure that at least max(1,M) columns are 92* supplied in the array Z; if RANGE = 'V', the exact value of M 93* is not known in advance and an upper bound must be used. 94* 95* LDZ (input) INTEGER 96* The leading dimension of the array Z. LDZ >= 1, and if 97* JOBZ = 'V', LDZ >= max(1,N). 98* 99* ISUPPZ (output) INTEGER ARRAY, dimension ( 2*max(1,M) ) 100* The support of the eigenvectors in Z, i.e., the indices 101* indicating the nonzero elements in Z. The i-th eigenvector 102* is nonzero only in elements ISUPPZ( 2*i-1 ) through 103* ISUPPZ( 2*i ). 104* 105* WORK (workspace) DOUBLE PRECISION array, dimension (13*N) 106* 107* IWORK (workspace) INTEGER array, dimension (6*N) 108* 109* INFO (output) INTEGER 110* = 0: successful exit 111* < 0: if INFO = -i, the i-th argument had an illegal value 112* > 0: if INFO = 1, internal error in DLARRB 113* if INFO = 2, internal error in DSTEIN 114* 115* Further Details 116* =============== 117* 118* Based on contributions by 119* Inderjit Dhillon, IBM Almaden, USA 120* Osni Marques, LBNL/NERSC, USA 121* Ken Stanley, Computer Science Division, University of 122* California at Berkeley, USA 123* 124* ===================================================================== 125* 126* .. Parameters .. 127 INTEGER MGSSIZ 128 PARAMETER ( MGSSIZ = 20 ) 129 DOUBLE PRECISION ZERO, ONE, FOUR 130 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, FOUR = 4.0D0 ) 131* .. 132* .. Local Scalars .. 133 LOGICAL MGSCLS 134 INTEGER I, IBEGIN, IEND, IINDC1, IINDC2, IINDR, IINDWK, 135 $ IINFO, IM, IN, INDERR, INDGAP, INDLD, INDLLD, 136 $ INDWRK, ITER, ITMP1, ITMP2, J, JBLK, K, KTOT, 137 $ LSBDPT, MAXITR, NCLUS, NDEPTH, NDONE, NEWCLS, 138 $ NEWFRS, NEWFTT, NEWLST, NEWSIZ, NSPLIT, OLDCLS, 139 $ OLDFST, OLDIEN, OLDLST, OLDNCL, P, Q 140 DOUBLE PRECISION EPS, GAP, LAMBDA, MGSTOL, MINGMA, MINRGP, 141 $ NRMINV, RELGAP, RELTOL, RESID, RQCORR, SIGMA, 142 $ TMP1, ZTZ 143* .. 144* .. External Functions .. 145 DOUBLE PRECISION DDOT, DLAMCH, DNRM2 146 EXTERNAL DDOT, DLAMCH, DNRM2 147* .. 148* .. External Subroutines .. 149 EXTERNAL DAXPY, DCOPY, DLAR1V, DLARRB, DLARRF, DLASET, 150 $ DSCAL, DSTEIN 151* .. 152* .. Intrinsic Functions .. 153 INTRINSIC ABS, DBLE, MAX, MIN, SQRT 154* .. 155* .. Local Arrays .. 156 INTEGER TEMP( 1 ) 157* .. 158* .. Executable Statements .. 159* 160* Test the input parameters. 161* 162 INDERR = N + 1 163 INDLD = 2*N 164 INDLLD = 3*N 165 INDGAP = 4*N 166 INDWRK = 5*N + 1 167* 168 IINDR = N 169 IINDC1 = 2*N 170 IINDC2 = 3*N 171 IINDWK = 4*N + 1 172* 173 EPS = DLAMCH( 'Precision' ) 174* 175 DO 10 I = 1, 2*N 176 IWORK( I ) = 0 177 10 CONTINUE 178 OPS = OPS + DBLE( M+1 ) 179 DO 20 I = 1, M 180 WORK( INDERR+I-1 ) = EPS*ABS( W( I ) ) 181 20 CONTINUE 182 CALL DLASET( 'Full', N, N, ZERO, ZERO, Z, LDZ ) 183 MGSTOL = 5.0D0*EPS 184* 185 NSPLIT = IBLOCK( M ) 186 IBEGIN = 1 187 DO 170 JBLK = 1, NSPLIT 188 IEND = ISPLIT( JBLK ) 189* 190* Find the eigenvectors of the submatrix indexed IBEGIN 191* through IEND. 192* 193 IF( IBEGIN.EQ.IEND ) THEN 194 Z( IBEGIN, IBEGIN ) = ONE 195 ISUPPZ( 2*IBEGIN-1 ) = IBEGIN 196 ISUPPZ( 2*IBEGIN ) = IBEGIN 197 IBEGIN = IEND + 1 198 GO TO 170 199 END IF 200 OLDIEN = IBEGIN - 1 201 IN = IEND - OLDIEN 202 OPS = OPS + DBLE( 1 ) 203 RELTOL = MIN( 1.0D-2, ONE / DBLE( IN ) ) 204 IM = IN 205 CALL DCOPY( IM, W( IBEGIN ), 1, WORK, 1 ) 206 OPS = OPS + DBLE( IN-1 ) 207 DO 30 I = 1, IN - 1 208 WORK( INDGAP+I ) = WORK( I+1 ) - WORK( I ) 209 30 CONTINUE 210 WORK( INDGAP+IN ) = MAX( ABS( WORK( IN ) ), EPS ) 211 NDONE = 0 212* 213 NDEPTH = 0 214 LSBDPT = 1 215 NCLUS = 1 216 IWORK( IINDC1+1 ) = 1 217 IWORK( IINDC1+2 ) = IN 218* 219* While( NDONE.LT.IM ) do 220* 221 40 CONTINUE 222 IF( NDONE.LT.IM ) THEN 223 OLDNCL = NCLUS 224 NCLUS = 0 225 LSBDPT = 1 - LSBDPT 226 DO 150 I = 1, OLDNCL 227 IF( LSBDPT.EQ.0 ) THEN 228 OLDCLS = IINDC1 229 NEWCLS = IINDC2 230 ELSE 231 OLDCLS = IINDC2 232 NEWCLS = IINDC1 233 END IF 234* 235* If NDEPTH > 1, retrieve the relatively robust 236* representation (RRR) and perform limited bisection 237* (if necessary) to get approximate eigenvalues. 238* 239 J = OLDCLS + 2*I 240 OLDFST = IWORK( J-1 ) 241 OLDLST = IWORK( J ) 242 IF( NDEPTH.GT.0 ) THEN 243 J = OLDIEN + OLDFST 244 CALL DCOPY( IN, Z( IBEGIN, J ), 1, D( IBEGIN ), 1 ) 245 CALL DCOPY( IN, Z( IBEGIN, J+1 ), 1, L( IBEGIN ), 1 ) 246 SIGMA = L( IEND ) 247 END IF 248 K = IBEGIN 249 OPS = OPS + DBLE( 2*( IN-1 ) ) 250 DO 50 J = 1, IN - 1 251 WORK( INDLD+J ) = D( K )*L( K ) 252 WORK( INDLLD+J ) = WORK( INDLD+J )*L( K ) 253 K = K + 1 254 50 CONTINUE 255 IF( NDEPTH.GT.0 ) THEN 256 CALL DLARRB( IN, D( IBEGIN ), L( IBEGIN ), 257 $ WORK( INDLD+1 ), WORK( INDLLD+1 ), 258 $ OLDFST, OLDLST, SIGMA, RELTOL, WORK, 259 $ WORK( INDGAP+1 ), WORK( INDERR ), 260 $ WORK( INDWRK ), IWORK( IINDWK ), IINFO ) 261 IF( IINFO.NE.0 ) THEN 262 INFO = 1 263 RETURN 264 END IF 265 END IF 266* 267* Classify eigenvalues of the current representation (RRR) 268* as (i) isolated, (ii) loosely clustered or (iii) tightly 269* clustered 270* 271 NEWFRS = OLDFST 272 DO 140 J = OLDFST, OLDLST 273 OPS = OPS + DBLE( 1 ) 274 IF( J.EQ.OLDLST .OR. WORK( INDGAP+J ).GE.RELTOL* 275 $ ABS( WORK( J ) ) ) THEN 276 NEWLST = J 277 ELSE 278* 279* continue (to the next loop) 280* 281 OPS = OPS + DBLE( 1 ) 282 RELGAP = WORK( INDGAP+J ) / ABS( WORK( J ) ) 283 IF( J.EQ.NEWFRS ) THEN 284 MINRGP = RELGAP 285 ELSE 286 MINRGP = MIN( MINRGP, RELGAP ) 287 END IF 288 GO TO 140 289 END IF 290 NEWSIZ = NEWLST - NEWFRS + 1 291 MAXITR = 10 292 NEWFTT = OLDIEN + NEWFRS 293 IF( NEWSIZ.GT.1 ) THEN 294 MGSCLS = NEWSIZ.LE.MGSSIZ .AND. MINRGP.GE.MGSTOL 295 IF( .NOT.MGSCLS ) THEN 296 CALL DLARRF( IN, D( IBEGIN ), L( IBEGIN ), 297 $ WORK( INDLD+1 ), WORK( INDLLD+1 ), 298 $ NEWFRS, NEWLST, WORK, 299 $ Z( IBEGIN, NEWFTT ), 300 $ Z( IBEGIN, NEWFTT+1 ), 301 $ WORK( INDWRK ), IWORK( IINDWK ), 302 $ INFO ) 303 IF( INFO.EQ.0 ) THEN 304 NCLUS = NCLUS + 1 305 K = NEWCLS + 2*NCLUS 306 IWORK( K-1 ) = NEWFRS 307 IWORK( K ) = NEWLST 308 ELSE 309 INFO = 0 310 IF( MINRGP.GE.MGSTOL ) THEN 311 MGSCLS = .TRUE. 312 ELSE 313* 314* Call DSTEIN to process this tight cluster. 315* This happens only if MINRGP <= MGSTOL 316* and DLARRF returns INFO = 1. The latter 317* means that a new RRR to "break" the 318* cluster could not be found. 319* 320 WORK( INDWRK ) = D( IBEGIN ) 321 OPS = OPS + DBLE( IN-1 ) 322 DO 60 K = 1, IN - 1 323 WORK( INDWRK+K ) = D( IBEGIN+K ) + 324 $ WORK( INDLLD+K ) 325 60 CONTINUE 326 DO 70 K = 1, NEWSIZ 327 IWORK( IINDWK+K-1 ) = 1 328 70 CONTINUE 329 DO 80 K = NEWFRS, NEWLST 330 ISUPPZ( 2*( IBEGIN+K )-3 ) = 1 331 ISUPPZ( 2*( IBEGIN+K )-2 ) = IN 332 80 CONTINUE 333 TEMP( 1 ) = IN 334 CALL DSTEIN( IN, WORK( INDWRK ), 335 $ WORK( INDLD+1 ), NEWSIZ, 336 $ WORK( NEWFRS ), 337 $ IWORK( IINDWK ), TEMP( 1 ), 338 $ Z( IBEGIN, NEWFTT ), LDZ, 339 $ WORK( INDWRK+IN ), 340 $ IWORK( IINDWK+IN ), 341 $ IWORK( IINDWK+2*IN ), IINFO ) 342 IF( IINFO.NE.0 ) THEN 343 INFO = 2 344 RETURN 345 END IF 346 NDONE = NDONE + NEWSIZ 347 END IF 348 END IF 349 END IF 350 ELSE 351 MGSCLS = .FALSE. 352 END IF 353 IF( NEWSIZ.EQ.1 .OR. MGSCLS ) THEN 354 KTOT = NEWFTT 355 DO 100 K = NEWFRS, NEWLST 356 ITER = 0 357 90 CONTINUE 358 LAMBDA = WORK( K ) 359 CALL DLAR1V( IN, 1, IN, LAMBDA, D( IBEGIN ), 360 $ L( IBEGIN ), WORK( INDLD+1 ), 361 $ WORK( INDLLD+1 ), 362 $ GERSCH( 2*OLDIEN+1 ), 363 $ Z( IBEGIN, KTOT ), ZTZ, MINGMA, 364 $ IWORK( IINDR+KTOT ), 365 $ ISUPPZ( 2*KTOT-1 ), 366 $ WORK( INDWRK ) ) 367 OPS = OPS + DBLE( 4 ) 368 TMP1 = ONE / ZTZ 369 NRMINV = SQRT( TMP1 ) 370 RESID = ABS( MINGMA )*NRMINV 371 RQCORR = MINGMA*TMP1 372 IF( K.EQ.IN ) THEN 373 GAP = WORK( INDGAP+K-1 ) 374 ELSE IF( K.EQ.1 ) THEN 375 GAP = WORK( INDGAP+K ) 376 ELSE 377 GAP = MIN( WORK( INDGAP+K-1 ), 378 $ WORK( INDGAP+K ) ) 379 END IF 380 ITER = ITER + 1 381 OPS = OPS + DBLE( 3 ) 382 IF( RESID.GT.TOL*GAP .AND. ABS( RQCORR ).GT. 383 $ FOUR*EPS*ABS( LAMBDA ) ) THEN 384 OPS = OPS + DBLE( 1 ) 385 WORK( K ) = LAMBDA + RQCORR 386 IF( ITER.LT.MAXITR ) THEN 387 GO TO 90 388 END IF 389 END IF 390 IWORK( KTOT ) = 1 391 IF( NEWSIZ.EQ.1 ) 392 $ NDONE = NDONE + 1 393 OPS = OPS + DBLE( IN ) 394 CALL DSCAL( IN, NRMINV, Z( IBEGIN, KTOT ), 1 ) 395 KTOT = KTOT + 1 396 100 CONTINUE 397 IF( NEWSIZ.GT.1 ) THEN 398 ITMP1 = ISUPPZ( 2*NEWFTT-1 ) 399 ITMP2 = ISUPPZ( 2*NEWFTT ) 400 KTOT = OLDIEN + NEWLST 401 DO 120 P = NEWFTT + 1, KTOT 402 DO 110 Q = NEWFTT, P - 1 403 OPS = OPS + DBLE( 4*IN ) 404 TMP1 = -DDOT( IN, Z( IBEGIN, P ), 1, 405 $ Z( IBEGIN, Q ), 1 ) 406 CALL DAXPY( IN, TMP1, Z( IBEGIN, Q ), 1, 407 $ Z( IBEGIN, P ), 1 ) 408 110 CONTINUE 409 OPS = OPS + DBLE( 3*IN+1 ) 410 TMP1 = ONE / DNRM2( IN, Z( IBEGIN, P ), 1 ) 411 CALL DSCAL( IN, TMP1, Z( IBEGIN, P ), 1 ) 412 ITMP1 = MIN( ITMP1, ISUPPZ( 2*P-1 ) ) 413 ITMP2 = MAX( ITMP2, ISUPPZ( 2*P ) ) 414 120 CONTINUE 415 DO 130 P = NEWFTT, KTOT 416 ISUPPZ( 2*P-1 ) = ITMP1 417 ISUPPZ( 2*P ) = ITMP2 418 130 CONTINUE 419 NDONE = NDONE + NEWSIZ 420 END IF 421 END IF 422 NEWFRS = J + 1 423 140 CONTINUE 424 150 CONTINUE 425 NDEPTH = NDEPTH + 1 426 GO TO 40 427 END IF 428 J = 2*IBEGIN 429 DO 160 I = IBEGIN, IEND 430 ISUPPZ( J-1 ) = ISUPPZ( J-1 ) + OLDIEN 431 ISUPPZ( J ) = ISUPPZ( J ) + OLDIEN 432 J = J + 2 433 160 CONTINUE 434 IBEGIN = IEND + 1 435 170 CONTINUE 436* 437 RETURN 438* 439* End of DLARRV 440* 441 END 442