1 SUBROUTINE DSTEIN( N, D, E, M, W, IBLOCK, ISPLIT, Z, LDZ, WORK, 2 $ IWORK, IFAIL, INFO ) 3* 4* -- LAPACK routine (instrumented to count operations, version 3.0) -- 5* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., 6* Courant Institute, Argonne National Lab, and Rice University 7* September 30, 1994 8* 9* .. Scalar Arguments .. 10 INTEGER INFO, LDZ, M, N 11* .. 12* .. Array Arguments .. 13 INTEGER IBLOCK( * ), IFAIL( * ), ISPLIT( * ), 14 $ IWORK( * ) 15 DOUBLE PRECISION D( * ), E( * ), W( * ), WORK( * ), Z( LDZ, * ) 16* .. 17* Common block to return operation count and iteration count 18* ITCNT is initialized to 0, OPS is only incremented 19* .. Common blocks .. 20 COMMON / LATIME / OPS, ITCNT 21* .. 22* .. Scalars in Common .. 23 DOUBLE PRECISION ITCNT, OPS 24* .. 25* 26* Purpose 27* ======= 28* 29* DSTEIN computes the eigenvectors of a real symmetric tridiagonal 30* matrix T corresponding to specified eigenvalues, using inverse 31* iteration. 32* 33* The maximum number of iterations allowed for each eigenvector is 34* specified by an internal parameter MAXITS (currently set to 5). 35* 36* Arguments 37* ========= 38* 39* N (input) INTEGER 40* The order of the matrix. N >= 0. 41* 42* D (input) DOUBLE PRECISION array, dimension (N) 43* The n diagonal elements of the tridiagonal matrix T. 44* 45* E (input) DOUBLE PRECISION array, dimension (N) 46* The (n-1) subdiagonal elements of the tridiagonal matrix 47* T, in elements 1 to N-1. E(N) need not be set. 48* 49* M (input) INTEGER 50* The number of eigenvectors to be found. 0 <= M <= N. 51* 52* W (input) DOUBLE PRECISION array, dimension (N) 53* The first M elements of W contain the eigenvalues for 54* which eigenvectors are to be computed. The eigenvalues 55* should be grouped by split-off block and ordered from 56* smallest to largest within the block. ( The output array 57* W from DSTEBZ with ORDER = 'B' is expected here. ) 58* 59* IBLOCK (input) INTEGER array, dimension (N) 60* The submatrix indices associated with the corresponding 61* eigenvalues in W; IBLOCK(i)=1 if eigenvalue W(i) belongs to 62* the first submatrix from the top, =2 if W(i) belongs to 63* the second submatrix, etc. ( The output array IBLOCK 64* from DSTEBZ is expected here. ) 65* 66* ISPLIT (input) INTEGER array, dimension (N) 67* The splitting points, at which T breaks up into submatrices. 68* The first submatrix consists of rows/columns 1 to 69* ISPLIT( 1 ), the second of rows/columns ISPLIT( 1 )+1 70* through ISPLIT( 2 ), etc. 71* ( The output array ISPLIT from DSTEBZ is expected here. ) 72* 73* Z (output) DOUBLE PRECISION array, dimension (LDZ, M) 74* The computed eigenvectors. The eigenvector associated 75* with the eigenvalue W(i) is stored in the i-th column of 76* Z. Any vector which fails to converge is set to its current 77* iterate after MAXITS iterations. 78* 79* LDZ (input) INTEGER 80* The leading dimension of the array Z. LDZ >= max(1,N). 81* 82* WORK (workspace) DOUBLE PRECISION array, dimension (5*N) 83* 84* IWORK (workspace) INTEGER array, dimension (N) 85* 86* IFAIL (output) INTEGER array, dimension (M) 87* On normal exit, all elements of IFAIL are zero. 88* If one or more eigenvectors fail to converge after 89* MAXITS iterations, then their indices are stored in 90* array IFAIL. 91* 92* INFO (output) INTEGER 93* = 0: successful exit. 94* < 0: if INFO = -i, the i-th argument had an illegal value 95* > 0: if INFO = i, then i eigenvectors failed to converge 96* in MAXITS iterations. Their indices are stored in 97* array IFAIL. 98* 99* Internal Parameters 100* =================== 101* 102* MAXITS INTEGER, default = 5 103* The maximum number of iterations performed. 104* 105* EXTRA INTEGER, default = 2 106* The number of iterations performed after norm growth 107* criterion is satisfied, should be at least 1. 108* 109* ===================================================================== 110* 111* .. Parameters .. 112 DOUBLE PRECISION ZERO, ONE, TEN, ODM3, ODM1 113 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0, TEN = 1.0D+1, 114 $ ODM3 = 1.0D-3, ODM1 = 1.0D-1 ) 115 INTEGER MAXITS, EXTRA 116 PARAMETER ( MAXITS = 5, EXTRA = 2 ) 117* .. 118* .. Local Scalars .. 119 INTEGER B1, BLKSIZ, BN, GPIND, I, IINFO, INDRV1, 120 $ INDRV2, INDRV3, INDRV4, INDRV5, ITS, J, J1, 121 $ JBLK, JMAX, NBLK, NRMCHK 122 DOUBLE PRECISION DTPCRT, EPS, EPS1, NRM, ONENRM, ORTOL, PERTOL, 123 $ SCL, SEP, TOL, XJ, XJM, ZTR 124* .. 125* .. Local Arrays .. 126 INTEGER ISEED( 4 ) 127* .. 128* .. External Functions .. 129 INTEGER IDAMAX 130 DOUBLE PRECISION DASUM, DDOT, DLAMCH, DNRM2 131 EXTERNAL IDAMAX, DASUM, DDOT, DLAMCH, DNRM2 132* .. 133* .. External Subroutines .. 134 EXTERNAL DAXPY, DCOPY, DLAGTF, DLAGTS, DLARNV, DSCAL, 135 $ XERBLA 136* .. 137* .. Intrinsic Functions .. 138 INTRINSIC ABS, MAX, SQRT 139* .. 140* .. Executable Statements .. 141* 142* Test the input parameters. 143* 144 INFO = 0 145 DO 10 I = 1, M 146 IFAIL( I ) = 0 147 10 CONTINUE 148* 149 IF( N.LT.0 ) THEN 150 INFO = -1 151 ELSE IF( M.LT.0 .OR. M.GT.N ) THEN 152 INFO = -4 153 ELSE IF( LDZ.LT.MAX( 1, N ) ) THEN 154 INFO = -9 155 ELSE 156 DO 20 J = 2, M 157 IF( IBLOCK( J ).LT.IBLOCK( J-1 ) ) THEN 158 INFO = -6 159 GO TO 30 160 END IF 161 IF( IBLOCK( J ).EQ.IBLOCK( J-1 ) .AND. W( J ).LT.W( J-1 ) ) 162 $ THEN 163 INFO = -5 164 GO TO 30 165 END IF 166 20 CONTINUE 167 30 CONTINUE 168 END IF 169* 170 IF( INFO.NE.0 ) THEN 171 CALL XERBLA( 'DSTEIN', -INFO ) 172 RETURN 173 END IF 174* 175* Initialize iteration count. 176* 177 ITCNT = 0 178* 179* Quick return if possible 180* 181 IF( N.EQ.0 .OR. M.EQ.0 ) THEN 182 RETURN 183 ELSE IF( N.EQ.1 ) THEN 184 Z( 1, 1 ) = ONE 185 RETURN 186 END IF 187* 188* Get machine constants. 189* 190 EPS = DLAMCH( 'Precision' ) 191* 192* Initialize seed for random number generator DLARNV. 193* 194 DO 40 I = 1, 4 195 ISEED( I ) = 1 196 40 CONTINUE 197* 198* Initialize pointers. 199* 200 INDRV1 = 0 201 INDRV2 = INDRV1 + N 202 INDRV3 = INDRV2 + N 203 INDRV4 = INDRV3 + N 204 INDRV5 = INDRV4 + N 205* 206* Compute eigenvectors of matrix blocks. 207* 208 J1 = 1 209 DO 160 NBLK = 1, IBLOCK( M ) 210* 211* Find starting and ending indices of block nblk. 212* 213 IF( NBLK.EQ.1 ) THEN 214 B1 = 1 215 ELSE 216 B1 = ISPLIT( NBLK-1 ) + 1 217 END IF 218 BN = ISPLIT( NBLK ) 219 BLKSIZ = BN - B1 + 1 220 IF( BLKSIZ.EQ.1 ) 221 $ GO TO 60 222 GPIND = B1 223* 224* Compute reorthogonalization criterion and stopping criterion. 225* 226 ONENRM = ABS( D( B1 ) ) + ABS( E( B1 ) ) 227 ONENRM = MAX( ONENRM, ABS( D( BN ) )+ABS( E( BN-1 ) ) ) 228 DO 50 I = B1 + 1, BN - 1 229 ONENRM = MAX( ONENRM, ABS( D( I ) )+ABS( E( I-1 ) )+ 230 $ ABS( E( I ) ) ) 231 50 CONTINUE 232 ORTOL = ODM3*ONENRM 233* 234 DTPCRT = SQRT( ODM1 / BLKSIZ ) 235* 236* Increment opcount for computing criteria. 237* 238 OPS = OPS + ( BN-B1 )*2 + 3 239* 240* Loop through eigenvalues of block nblk. 241* 242 60 CONTINUE 243 JBLK = 0 244 DO 150 J = J1, M 245 IF( IBLOCK( J ).NE.NBLK ) THEN 246 J1 = J 247 GO TO 160 248 END IF 249 JBLK = JBLK + 1 250 XJ = W( J ) 251* 252* Skip all the work if the block size is one. 253* 254 IF( BLKSIZ.EQ.1 ) THEN 255 WORK( INDRV1+1 ) = ONE 256 GO TO 120 257 END IF 258* 259* If eigenvalues j and j-1 are too close, add a relatively 260* small perturbation. 261* 262 IF( JBLK.GT.1 ) THEN 263 EPS1 = ABS( EPS*XJ ) 264 PERTOL = TEN*EPS1 265 SEP = XJ - XJM 266 IF( SEP.LT.PERTOL ) 267 $ XJ = XJM + PERTOL 268 END IF 269* 270 ITS = 0 271 NRMCHK = 0 272* 273* Get random starting vector. 274* 275 CALL DLARNV( 2, ISEED, BLKSIZ, WORK( INDRV1+1 ) ) 276* 277* Increment opcount for getting random starting vector. 278* ( DLARND(2,.) requires 9 flops. ) 279* 280 OPS = OPS + BLKSIZ*9 281* 282* Copy the matrix T so it won't be destroyed in factorization. 283* 284 CALL DCOPY( BLKSIZ, D( B1 ), 1, WORK( INDRV4+1 ), 1 ) 285 CALL DCOPY( BLKSIZ-1, E( B1 ), 1, WORK( INDRV2+2 ), 1 ) 286 CALL DCOPY( BLKSIZ-1, E( B1 ), 1, WORK( INDRV3+1 ), 1 ) 287* 288* Compute LU factors with partial pivoting ( PT = LU ) 289* 290 TOL = ZERO 291 CALL DLAGTF( BLKSIZ, WORK( INDRV4+1 ), XJ, WORK( INDRV2+2 ), 292 $ WORK( INDRV3+1 ), TOL, WORK( INDRV5+1 ), IWORK, 293 $ IINFO ) 294* 295* Increment opcount for computing LU factors. 296* ( DLAGTF(BLKSIZ,...) requires about 8*BLKSIZ flops. ) 297* 298 OPS = OPS + 8*BLKSIZ 299* 300* Update iteration count. 301* 302 70 CONTINUE 303 ITS = ITS + 1 304 IF( ITS.GT.MAXITS ) 305 $ GO TO 100 306* 307* Normalize and scale the righthand side vector Pb. 308* 309 SCL = BLKSIZ*ONENRM*MAX( EPS, 310 $ ABS( WORK( INDRV4+BLKSIZ ) ) ) / 311 $ DASUM( BLKSIZ, WORK( INDRV1+1 ), 1 ) 312 CALL DSCAL( BLKSIZ, SCL, WORK( INDRV1+1 ), 1 ) 313* 314* Solve the system LU = Pb. 315* 316 CALL DLAGTS( -1, BLKSIZ, WORK( INDRV4+1 ), WORK( INDRV2+2 ), 317 $ WORK( INDRV3+1 ), WORK( INDRV5+1 ), IWORK, 318 $ WORK( INDRV1+1 ), TOL, IINFO ) 319* 320* Increment opcount for scaling and solving linear system. 321* ( DLAGTS(-1,BLKSIZ,...) requires about 8*BLKSIZ flops. ) 322* 323 OPS = OPS + 3 + 10*BLKSIZ 324* 325* Reorthogonalize by modified Gram-Schmidt if eigenvalues are 326* close enough. 327* 328 IF( JBLK.EQ.1 ) 329 $ GO TO 90 330 IF( ABS( XJ-XJM ).GT.ORTOL ) 331 $ GPIND = J 332 IF( GPIND.NE.J ) THEN 333 DO 80 I = GPIND, J - 1 334 ZTR = -DDOT( BLKSIZ, WORK( INDRV1+1 ), 1, Z( B1, I ), 335 $ 1 ) 336 CALL DAXPY( BLKSIZ, ZTR, Z( B1, I ), 1, 337 $ WORK( INDRV1+1 ), 1 ) 338 80 CONTINUE 339* 340* Increment opcount for reorthogonalizing. 341* 342 OPS = OPS + ( J-GPIND )*BLKSIZ*4 343* 344 END IF 345* 346* Check the infinity norm of the iterate. 347* 348 90 CONTINUE 349 JMAX = IDAMAX( BLKSIZ, WORK( INDRV1+1 ), 1 ) 350 NRM = ABS( WORK( INDRV1+JMAX ) ) 351* 352* Continue for additional iterations after norm reaches 353* stopping criterion. 354* 355 IF( NRM.LT.DTPCRT ) 356 $ GO TO 70 357 NRMCHK = NRMCHK + 1 358 IF( NRMCHK.LT.EXTRA+1 ) 359 $ GO TO 70 360* 361 GO TO 110 362* 363* If stopping criterion was not satisfied, update info and 364* store eigenvector number in array ifail. 365* 366 100 CONTINUE 367 INFO = INFO + 1 368 IFAIL( INFO ) = J 369* 370* Accept iterate as jth eigenvector. 371* 372 110 CONTINUE 373 SCL = ONE / DNRM2( BLKSIZ, WORK( INDRV1+1 ), 1 ) 374 JMAX = IDAMAX( BLKSIZ, WORK( INDRV1+1 ), 1 ) 375 IF( WORK( INDRV1+JMAX ).LT.ZERO ) 376 $ SCL = -SCL 377 CALL DSCAL( BLKSIZ, SCL, WORK( INDRV1+1 ), 1 ) 378* 379* Increment opcount for scaling. 380* 381 OPS = OPS + 3*BLKSIZ 382* 383 120 CONTINUE 384 DO 130 I = 1, N 385 Z( I, J ) = ZERO 386 130 CONTINUE 387 DO 140 I = 1, BLKSIZ 388 Z( B1+I-1, J ) = WORK( INDRV1+I ) 389 140 CONTINUE 390* 391* Save the shift to check eigenvalue spacing at next 392* iteration. 393* 394 XJM = XJ 395* 396 150 CONTINUE 397 160 CONTINUE 398* 399 RETURN 400* 401* End of DSTEIN 402* 403 END 404