1 SUBROUTINE CLARRV( 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 REAL TOL 12* .. 13* .. Array Arguments .. 14 INTEGER IBLOCK( * ), ISPLIT( * ), ISUPPZ( * ), 15 $ IWORK( * ) 16 REAL D( * ), GERSCH( * ), L( * ), W( * ), WORK( * ) 17 COMPLEX Z( LDZ, * ) 18* .. 19* Common block to return operation count .. 20* .. Common blocks .. 21 COMMON / LATIME / OPS, ITCNT 22* .. 23* .. Scalars in Common .. 24 REAL ITCNT, OPS 25* .. 26* 27* Purpose 28* ======= 29* 30* CLARRV 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) REAL 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) REAL 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) REAL 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) REAL 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 SLARRE 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) COMPLEX 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) REAL 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 SLARRB 113* if INFO = 2, internal error in CSTEIN 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 REAL ZERO, ONE, FOUR 130 PARAMETER ( ZERO = 0.0E0, ONE = 1.0E0, FOUR = 4.0E0 ) 131 COMPLEX CZERO 132 PARAMETER ( CZERO = ( 0.0E0, 0.0E0 ) ) 133* .. 134* .. Local Scalars .. 135 LOGICAL MGSCLS 136 INTEGER I, IBEGIN, IEND, IINDC1, IINDC2, IINDR, IINDWK, 137 $ IINFO, IM, IN, INDERR, INDGAP, INDIN1, INDIN2, 138 $ INDLD, INDLLD, INDWRK, ITER, ITMP1, ITMP2, J, 139 $ JBLK, K, KTOT, LSBDPT, MAXITR, NCLUS, NDEPTH, 140 $ NDONE, NEWCLS, NEWFRS, NEWFTT, NEWLST, NEWSIZ, 141 $ NSPLIT, OLDCLS, OLDFST, OLDIEN, OLDLST, OLDNCL, 142 $ P, Q, TEMP( 1 ) 143 REAL EPS, GAP, LAMBDA, MGSTOL, MINGMA, MINRGP, 144 $ NRMINV, RELGAP, RELTOL, RESID, RQCORR, SIGMA, 145 $ TMP1, ZTZ 146* .. 147* .. External Functions .. 148 REAL SCNRM2, SLAMCH 149 COMPLEX CDOTU 150 EXTERNAL CDOTU, SCNRM2, SLAMCH 151* .. 152* .. External Subroutines .. 153 EXTERNAL CAXPY, CLAR1V, CLASET, CSTEIN, SCOPY, SLARRB 154* .. 155* .. Intrinsic Functions .. 156 INTRINSIC ABS, CMPLX, MAX, MIN, REAL, SQRT 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 INDIN1 = 5*N + 1 167 INDIN2 = 6*N + 1 168 INDWRK = 7*N + 1 169* 170 IINDR = N 171 IINDC1 = 2*N 172 IINDC2 = 3*N 173 IINDWK = 4*N + 1 174* 175 EPS = SLAMCH( 'Precision' ) 176* 177 DO 10 I = 1, 2*N 178 IWORK( I ) = 0 179 10 CONTINUE 180 OPS = OPS + REAL( M+1 ) 181 DO 20 I = 1, M 182 WORK( INDERR+I-1 ) = EPS*ABS( W( I ) ) 183 20 CONTINUE 184 CALL CLASET( 'Full', N, N, CZERO, CZERO, Z, LDZ ) 185 MGSTOL = 5.0E0*EPS 186* 187 NSPLIT = IBLOCK( M ) 188 IBEGIN = 1 189 DO 170 JBLK = 1, NSPLIT 190 IEND = ISPLIT( JBLK ) 191* 192* Find the eigenvectors of the submatrix indexed IBEGIN 193* through IEND. 194* 195 IF( IBEGIN.EQ.IEND ) THEN 196 Z( IBEGIN, IBEGIN ) = ONE 197 ISUPPZ( 2*IBEGIN-1 ) = IBEGIN 198 ISUPPZ( 2*IBEGIN ) = IBEGIN 199 IBEGIN = IEND + 1 200 GO TO 170 201 END IF 202 OLDIEN = IBEGIN - 1 203 IN = IEND - OLDIEN 204 OPS = OPS + REAL( 1 ) 205 RELTOL = MIN( 1.0E-2, ONE / REAL( IN ) ) 206 IM = IN 207 CALL SCOPY( IM, W( IBEGIN ), 1, WORK, 1 ) 208 OPS = OPS + REAL( IN-1 ) 209 DO 30 I = 1, IN - 1 210 WORK( INDGAP+I ) = WORK( I+1 ) - WORK( I ) 211 30 CONTINUE 212 WORK( INDGAP+IN ) = MAX( ABS( WORK( IN ) ), EPS ) 213 NDONE = 0 214* 215 NDEPTH = 0 216 LSBDPT = 1 217 NCLUS = 1 218 IWORK( IINDC1+1 ) = 1 219 IWORK( IINDC1+2 ) = IN 220* 221* While( NDONE.LT.IM ) do 222* 223 40 CONTINUE 224 IF( NDONE.LT.IM ) THEN 225 OLDNCL = NCLUS 226 NCLUS = 0 227 LSBDPT = 1 - LSBDPT 228 DO 150 I = 1, OLDNCL 229 IF( LSBDPT.EQ.0 ) THEN 230 OLDCLS = IINDC1 231 NEWCLS = IINDC2 232 ELSE 233 OLDCLS = IINDC2 234 NEWCLS = IINDC1 235 END IF 236* 237* If NDEPTH > 1, retrieve the relatively robust 238* representation (RRR) and perform limited bisection 239* (if necessary) to get approximate eigenvalues. 240* 241 J = OLDCLS + 2*I 242 OLDFST = IWORK( J-1 ) 243 OLDLST = IWORK( J ) 244 IF( NDEPTH.GT.0 ) THEN 245 J = OLDIEN + OLDFST 246 DO 45 K = 1, IN 247 D( IBEGIN+K-1 ) = REAL( Z( IBEGIN+K-1, 248 $ OLDIEN+OLDFST ) ) 249 L( IBEGIN+K-1 ) = REAL( Z( IBEGIN+K-1, 250 $ OLDIEN+OLDFST+1 ) ) 251 45 CONTINUE 252 SIGMA = L( IEND ) 253 END IF 254 K = IBEGIN 255 OPS = OPS + REAL( 2*(IN-1) ) 256 DO 50 J = 1, IN - 1 257 WORK( INDLD+J ) = D( K )*L( K ) 258 WORK( INDLLD+J ) = WORK( INDLD+J )*L( K ) 259 K = K + 1 260 50 CONTINUE 261 IF( NDEPTH.GT.0 ) THEN 262 CALL SLARRB( IN, D( IBEGIN ), L( IBEGIN ), 263 $ WORK( INDLD+1 ), WORK( INDLLD+1 ), 264 $ OLDFST, OLDLST, SIGMA, RELTOL, WORK, 265 $ WORK( INDGAP+1 ), WORK( INDERR ), 266 $ WORK( INDWRK ), IWORK( IINDWK ), IINFO ) 267 IF( IINFO.NE.0 ) THEN 268 INFO = 1 269 RETURN 270 END IF 271 END IF 272* 273* Classify eigenvalues of the current representation (RRR) 274* as (i) isolated, (ii) loosely clustered or (iii) tightly 275* clustered 276* 277 NEWFRS = OLDFST 278 DO 140 J = OLDFST, OLDLST 279 OPS = OPS + REAL( 1 ) 280 IF( J.EQ.OLDLST .OR. WORK( INDGAP+J ).GE.RELTOL* 281 $ ABS( WORK( J ) ) ) THEN 282 NEWLST = J 283 ELSE 284* 285* continue (to the next loop) 286* 287 OPS = OPS + REAL( 1 ) 288 RELGAP = WORK( INDGAP+J ) / ABS( WORK( J ) ) 289 IF( J.EQ.NEWFRS ) THEN 290 MINRGP = RELGAP 291 ELSE 292 MINRGP = MIN( MINRGP, RELGAP ) 293 END IF 294 GO TO 140 295 END IF 296 NEWSIZ = NEWLST - NEWFRS + 1 297 MAXITR = 10 298 NEWFTT = OLDIEN + NEWFRS 299 IF( NEWSIZ.GT.1 ) THEN 300 MGSCLS = NEWSIZ.LE.MGSSIZ .AND. MINRGP.GE.MGSTOL 301 IF( .NOT.MGSCLS ) THEN 302 DO 55 K = 1, IN 303 WORK( INDIN1+K-1 ) = REAL( Z( IBEGIN+K-1, 304 $ NEWFTT ) ) 305 WORK( INDIN2+K-1 ) = REAL( Z( IBEGIN+K-1, 306 $ NEWFTT+1 ) ) 307 55 CONTINUE 308 CALL SLARRF( IN, D( IBEGIN ), L( IBEGIN ), 309 $ WORK( INDLD+1 ), WORK( INDLLD+1 ), 310 $ NEWFRS, NEWLST, WORK, 311 $ WORK( INDIN1 ), WORK( INDIN2 ), 312 $ WORK( INDWRK ), IWORK( IINDWK ), 313 $ INFO ) 314 IF( INFO.EQ.0 ) THEN 315 NCLUS = NCLUS + 1 316 K = NEWCLS + 2*NCLUS 317 IWORK( K-1 ) = NEWFRS 318 IWORK( K ) = NEWLST 319 ELSE 320 INFO = 0 321 IF( MINRGP.GE.MGSTOL ) THEN 322 MGSCLS = .TRUE. 323 ELSE 324* 325* Call CSTEIN to process this tight cluster. 326* This happens only if MINRGP <= MGSTOL 327* and SLARRF returns INFO = 1. The latter 328* means that a new RRR to "break" the 329* cluster could not be found. 330* 331 WORK( INDWRK ) = D( IBEGIN ) 332 OPS = OPS + REAL( IN-1 ) 333 DO 60 K = 1, IN - 1 334 WORK( INDWRK+K ) = D( IBEGIN+K ) + 335 $ WORK( INDLLD+K ) 336 60 CONTINUE 337 DO 70 K = 1, NEWSIZ 338 IWORK( IINDWK+K-1 ) = 1 339 70 CONTINUE 340 DO 80 K = NEWFRS, NEWLST 341 ISUPPZ( 2*( IBEGIN+K )-3 ) = 1 342 ISUPPZ( 2*( IBEGIN+K )-2 ) = IN 343 80 CONTINUE 344 TEMP( 1 ) = IN 345 CALL CSTEIN( IN, WORK( INDWRK ), 346 $ WORK( INDLD+1 ), NEWSIZ, 347 $ WORK( NEWFRS ), 348 $ IWORK( IINDWK ), TEMP( 1 ), 349 $ Z( IBEGIN, NEWFTT ), LDZ, 350 $ WORK( INDWRK+IN ), 351 $ IWORK( IINDWK+IN ), 352 $ IWORK( IINDWK+2*IN ), IINFO ) 353 IF( IINFO.NE.0 ) THEN 354 INFO = 2 355 RETURN 356 END IF 357 NDONE = NDONE + NEWSIZ 358 END IF 359 END IF 360 END IF 361 ELSE 362 MGSCLS = .FALSE. 363 END IF 364 IF( NEWSIZ.EQ.1 .OR. MGSCLS ) THEN 365 KTOT = NEWFTT 366 DO 100 K = NEWFRS, NEWLST 367 ITER = 0 368 90 CONTINUE 369 LAMBDA = WORK( K ) 370 CALL CLAR1V( IN, 1, IN, LAMBDA, D( IBEGIN ), 371 $ L( IBEGIN ), WORK( INDLD+1 ), 372 $ WORK( INDLLD+1 ), 373 $ GERSCH( 2*OLDIEN+1 ), 374 $ Z( IBEGIN, KTOT ), ZTZ, MINGMA, 375 $ IWORK( IINDR+KTOT ), 376 $ ISUPPZ( 2*KTOT-1 ), 377 $ WORK( INDWRK ) ) 378 OPS = OPS + REAL( 4 ) 379 TMP1 = ONE / ZTZ 380 NRMINV = SQRT( TMP1 ) 381 RESID = ABS( MINGMA )*NRMINV 382 RQCORR = MINGMA*TMP1 383 IF( K.EQ.IN ) THEN 384 GAP = WORK( INDGAP+K-1 ) 385 ELSE IF( K.EQ.1 ) THEN 386 GAP = WORK( INDGAP+K ) 387 ELSE 388 GAP = MIN( WORK( INDGAP+K-1 ), 389 $ WORK( INDGAP+K ) ) 390 END IF 391 ITER = ITER + 1 392 OPS = OPS + REAL( 3 ) 393 IF( RESID.GT.TOL*GAP .AND. ABS( RQCORR ).GT. 394 $ FOUR*EPS*ABS( LAMBDA ) ) THEN 395 OPS = OPS + REAL( 1 ) 396 WORK( K ) = LAMBDA + RQCORR 397 IF( ITER.LT.MAXITR ) THEN 398 GO TO 90 399 END IF 400 END IF 401 IWORK( KTOT ) = 1 402 IF( NEWSIZ.EQ.1 ) 403 $ NDONE = NDONE + 1 404 OPS = OPS + REAL( 2*IN ) 405 CALL CSSCAL( IN, NRMINV, Z( IBEGIN, KTOT ), 1 ) 406 KTOT = KTOT + 1 407 100 CONTINUE 408 IF( NEWSIZ.GT.1 ) THEN 409 ITMP1 = ISUPPZ( 2*NEWFTT-1 ) 410 ITMP2 = ISUPPZ( 2*NEWFTT ) 411 KTOT = OLDIEN + NEWLST 412 DO 120 P = NEWFTT + 1, KTOT 413 DO 110 Q = NEWFTT, P - 1 414 OPS = OPS + REAL( 10*IN ) 415 TMP1 = -CDOTU( IN, Z( IBEGIN, P ), 1, 416 $ Z( IBEGIN, Q ), 1 ) 417 CALL CAXPY( IN, CMPLX( TMP1, ZERO ), 418 $ Z( IBEGIN, Q ), 1, 419 $ Z( IBEGIN, P ), 1 ) 420 110 CONTINUE 421 OPS = OPS + REAL( 8*IN+1 ) 422 TMP1 = ONE / SCNRM2( IN, Z( IBEGIN, P ), 1 ) 423 CALL CSSCAL( IN, TMP1, Z( IBEGIN, P ), 1 ) 424 ITMP1 = MIN( ITMP1, ISUPPZ( 2*P-1 ) ) 425 ITMP2 = MAX( ITMP2, ISUPPZ( 2*P ) ) 426 120 CONTINUE 427 DO 130 P = NEWFTT, KTOT 428 ISUPPZ( 2*P-1 ) = ITMP1 429 ISUPPZ( 2*P ) = ITMP2 430 130 CONTINUE 431 NDONE = NDONE + NEWSIZ 432 END IF 433 END IF 434 NEWFRS = J + 1 435 140 CONTINUE 436 150 CONTINUE 437 NDEPTH = NDEPTH + 1 438 GO TO 40 439 END IF 440 J = 2*IBEGIN 441 DO 160 I = IBEGIN, IEND 442 ISUPPZ( J-1 ) = ISUPPZ( J-1 ) + OLDIEN 443 ISUPPZ( J ) = ISUPPZ( J ) + OLDIEN 444 J = J + 2 445 160 CONTINUE 446 IBEGIN = IEND + 1 447 170 CONTINUE 448* 449 RETURN 450* 451* End of CLARRV 452* 453 END 454