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