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