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