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