|
LAPACK 3.3.0
|
|
LAPACK 3.3.0
|
00001 SUBROUTINE DSTEIN( N, D, E, M, W, IBLOCK, ISPLIT, Z, LDZ, WORK, 00002 $ IWORK, IFAIL, INFO ) 00003 * 00004 * -- LAPACK routine (version 3.2) -- 00005 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00006 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00007 * November 2006 00008 * 00009 * .. Scalar Arguments .. 00010 INTEGER INFO, LDZ, M, N 00011 * .. 00012 * .. Array Arguments .. 00013 INTEGER IBLOCK( * ), IFAIL( * ), ISPLIT( * ), 00014 $ IWORK( * ) 00015 DOUBLE PRECISION D( * ), E( * ), W( * ), WORK( * ), Z( LDZ, * ) 00016 * .. 00017 * 00018 * Purpose 00019 * ======= 00020 * 00021 * DSTEIN computes the eigenvectors of a real symmetric tridiagonal 00022 * matrix T corresponding to specified eigenvalues, using inverse 00023 * iteration. 00024 * 00025 * The maximum number of iterations allowed for each eigenvector is 00026 * specified by an internal parameter MAXITS (currently set to 5). 00027 * 00028 * Arguments 00029 * ========= 00030 * 00031 * N (input) INTEGER 00032 * The order of the matrix. N >= 0. 00033 * 00034 * D (input) DOUBLE PRECISION array, dimension (N) 00035 * The n diagonal elements of the tridiagonal matrix T. 00036 * 00037 * E (input) DOUBLE PRECISION array, dimension (N-1) 00038 * The (n-1) subdiagonal elements of the tridiagonal matrix 00039 * T, in elements 1 to N-1. 00040 * 00041 * M (input) INTEGER 00042 * The number of eigenvectors to be found. 0 <= M <= N. 00043 * 00044 * W (input) DOUBLE PRECISION array, dimension (N) 00045 * The first M elements of W contain the eigenvalues for 00046 * which eigenvectors are to be computed. The eigenvalues 00047 * should be grouped by split-off block and ordered from 00048 * smallest to largest within the block. ( The output array 00049 * W from DSTEBZ with ORDER = 'B' is expected here. ) 00050 * 00051 * IBLOCK (input) INTEGER array, dimension (N) 00052 * The submatrix indices associated with the corresponding 00053 * eigenvalues in W; IBLOCK(i)=1 if eigenvalue W(i) belongs to 00054 * the first submatrix from the top, =2 if W(i) belongs to 00055 * the second submatrix, etc. ( The output array IBLOCK 00056 * from DSTEBZ is expected here. ) 00057 * 00058 * ISPLIT (input) INTEGER array, dimension (N) 00059 * The splitting points, at which T breaks up into submatrices. 00060 * The first submatrix consists of rows/columns 1 to 00061 * ISPLIT( 1 ), the second of rows/columns ISPLIT( 1 )+1 00062 * through ISPLIT( 2 ), etc. 00063 * ( The output array ISPLIT from DSTEBZ is expected here. ) 00064 * 00065 * Z (output) DOUBLE PRECISION array, dimension (LDZ, M) 00066 * The computed eigenvectors. The eigenvector associated 00067 * with the eigenvalue W(i) is stored in the i-th column of 00068 * Z. Any vector which fails to converge is set to its current 00069 * iterate after MAXITS iterations. 00070 * 00071 * LDZ (input) INTEGER 00072 * The leading dimension of the array Z. LDZ >= max(1,N). 00073 * 00074 * WORK (workspace) DOUBLE PRECISION array, dimension (5*N) 00075 * 00076 * IWORK (workspace) INTEGER array, dimension (N) 00077 * 00078 * IFAIL (output) INTEGER array, dimension (M) 00079 * On normal exit, all elements of IFAIL are zero. 00080 * If one or more eigenvectors fail to converge after 00081 * MAXITS iterations, then their indices are stored in 00082 * array IFAIL. 00083 * 00084 * INFO (output) INTEGER 00085 * = 0: successful exit. 00086 * < 0: if INFO = -i, the i-th argument had an illegal value 00087 * > 0: if INFO = i, then i eigenvectors failed to converge 00088 * in MAXITS iterations. Their indices are stored in 00089 * array IFAIL. 00090 * 00091 * Internal Parameters 00092 * =================== 00093 * 00094 * MAXITS INTEGER, default = 5 00095 * The maximum number of iterations performed. 00096 * 00097 * EXTRA INTEGER, default = 2 00098 * The number of iterations performed after norm growth 00099 * criterion is satisfied, should be at least 1. 00100 * 00101 * ===================================================================== 00102 * 00103 * .. Parameters .. 00104 DOUBLE PRECISION ZERO, ONE, TEN, ODM3, ODM1 00105 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0, TEN = 1.0D+1, 00106 $ ODM3 = 1.0D-3, ODM1 = 1.0D-1 ) 00107 INTEGER MAXITS, EXTRA 00108 PARAMETER ( MAXITS = 5, EXTRA = 2 ) 00109 * .. 00110 * .. Local Scalars .. 00111 INTEGER B1, BLKSIZ, BN, GPIND, I, IINFO, INDRV1, 00112 $ INDRV2, INDRV3, INDRV4, INDRV5, ITS, J, J1, 00113 $ JBLK, JMAX, NBLK, NRMCHK 00114 DOUBLE PRECISION DTPCRT, EPS, EPS1, NRM, ONENRM, ORTOL, PERTOL, 00115 $ SCL, SEP, TOL, XJ, XJM, ZTR 00116 * .. 00117 * .. Local Arrays .. 00118 INTEGER ISEED( 4 ) 00119 * .. 00120 * .. External Functions .. 00121 INTEGER IDAMAX 00122 DOUBLE PRECISION DASUM, DDOT, DLAMCH, DNRM2 00123 EXTERNAL IDAMAX, DASUM, DDOT, DLAMCH, DNRM2 00124 * .. 00125 * .. External Subroutines .. 00126 EXTERNAL DAXPY, DCOPY, DLAGTF, DLAGTS, DLARNV, DSCAL, 00127 $ XERBLA 00128 * .. 00129 * .. Intrinsic Functions .. 00130 INTRINSIC ABS, MAX, SQRT 00131 * .. 00132 * .. Executable Statements .. 00133 * 00134 * Test the input parameters. 00135 * 00136 INFO = 0 00137 DO 10 I = 1, M 00138 IFAIL( I ) = 0 00139 10 CONTINUE 00140 * 00141 IF( N.LT.0 ) THEN 00142 INFO = -1 00143 ELSE IF( M.LT.0 .OR. M.GT.N ) THEN 00144 INFO = -4 00145 ELSE IF( LDZ.LT.MAX( 1, N ) ) THEN 00146 INFO = -9 00147 ELSE 00148 DO 20 J = 2, M 00149 IF( IBLOCK( J ).LT.IBLOCK( J-1 ) ) THEN 00150 INFO = -6 00151 GO TO 30 00152 END IF 00153 IF( IBLOCK( J ).EQ.IBLOCK( J-1 ) .AND. W( J ).LT.W( J-1 ) ) 00154 $ THEN 00155 INFO = -5 00156 GO TO 30 00157 END IF 00158 20 CONTINUE 00159 30 CONTINUE 00160 END IF 00161 * 00162 IF( INFO.NE.0 ) THEN 00163 CALL XERBLA( 'DSTEIN', -INFO ) 00164 RETURN 00165 END IF 00166 * 00167 * Quick return if possible 00168 * 00169 IF( N.EQ.0 .OR. M.EQ.0 ) THEN 00170 RETURN 00171 ELSE IF( N.EQ.1 ) THEN 00172 Z( 1, 1 ) = ONE 00173 RETURN 00174 END IF 00175 * 00176 * Get machine constants. 00177 * 00178 EPS = DLAMCH( 'Precision' ) 00179 * 00180 * Initialize seed for random number generator DLARNV. 00181 * 00182 DO 40 I = 1, 4 00183 ISEED( I ) = 1 00184 40 CONTINUE 00185 * 00186 * Initialize pointers. 00187 * 00188 INDRV1 = 0 00189 INDRV2 = INDRV1 + N 00190 INDRV3 = INDRV2 + N 00191 INDRV4 = INDRV3 + N 00192 INDRV5 = INDRV4 + N 00193 * 00194 * Compute eigenvectors of matrix blocks. 00195 * 00196 J1 = 1 00197 DO 160 NBLK = 1, IBLOCK( M ) 00198 * 00199 * Find starting and ending indices of block nblk. 00200 * 00201 IF( NBLK.EQ.1 ) THEN 00202 B1 = 1 00203 ELSE 00204 B1 = ISPLIT( NBLK-1 ) + 1 00205 END IF 00206 BN = ISPLIT( NBLK ) 00207 BLKSIZ = BN - B1 + 1 00208 IF( BLKSIZ.EQ.1 ) 00209 $ GO TO 60 00210 GPIND = B1 00211 * 00212 * Compute reorthogonalization criterion and stopping criterion. 00213 * 00214 ONENRM = ABS( D( B1 ) ) + ABS( E( B1 ) ) 00215 ONENRM = MAX( ONENRM, ABS( D( BN ) )+ABS( E( BN-1 ) ) ) 00216 DO 50 I = B1 + 1, BN - 1 00217 ONENRM = MAX( ONENRM, ABS( D( I ) )+ABS( E( I-1 ) )+ 00218 $ ABS( E( I ) ) ) 00219 50 CONTINUE 00220 ORTOL = ODM3*ONENRM 00221 * 00222 DTPCRT = SQRT( ODM1 / BLKSIZ ) 00223 * 00224 * Loop through eigenvalues of block nblk. 00225 * 00226 60 CONTINUE 00227 JBLK = 0 00228 DO 150 J = J1, M 00229 IF( IBLOCK( J ).NE.NBLK ) THEN 00230 J1 = J 00231 GO TO 160 00232 END IF 00233 JBLK = JBLK + 1 00234 XJ = W( J ) 00235 * 00236 * Skip all the work if the block size is one. 00237 * 00238 IF( BLKSIZ.EQ.1 ) THEN 00239 WORK( INDRV1+1 ) = ONE 00240 GO TO 120 00241 END IF 00242 * 00243 * If eigenvalues j and j-1 are too close, add a relatively 00244 * small perturbation. 00245 * 00246 IF( JBLK.GT.1 ) THEN 00247 EPS1 = ABS( EPS*XJ ) 00248 PERTOL = TEN*EPS1 00249 SEP = XJ - XJM 00250 IF( SEP.LT.PERTOL ) 00251 $ XJ = XJM + PERTOL 00252 END IF 00253 * 00254 ITS = 0 00255 NRMCHK = 0 00256 * 00257 * Get random starting vector. 00258 * 00259 CALL DLARNV( 2, ISEED, BLKSIZ, WORK( INDRV1+1 ) ) 00260 * 00261 * Copy the matrix T so it won't be destroyed in factorization. 00262 * 00263 CALL DCOPY( BLKSIZ, D( B1 ), 1, WORK( INDRV4+1 ), 1 ) 00264 CALL DCOPY( BLKSIZ-1, E( B1 ), 1, WORK( INDRV2+2 ), 1 ) 00265 CALL DCOPY( BLKSIZ-1, E( B1 ), 1, WORK( INDRV3+1 ), 1 ) 00266 * 00267 * Compute LU factors with partial pivoting ( PT = LU ) 00268 * 00269 TOL = ZERO 00270 CALL DLAGTF( BLKSIZ, WORK( INDRV4+1 ), XJ, WORK( INDRV2+2 ), 00271 $ WORK( INDRV3+1 ), TOL, WORK( INDRV5+1 ), IWORK, 00272 $ IINFO ) 00273 * 00274 * Update iteration count. 00275 * 00276 70 CONTINUE 00277 ITS = ITS + 1 00278 IF( ITS.GT.MAXITS ) 00279 $ GO TO 100 00280 * 00281 * Normalize and scale the righthand side vector Pb. 00282 * 00283 SCL = BLKSIZ*ONENRM*MAX( EPS, 00284 $ ABS( WORK( INDRV4+BLKSIZ ) ) ) / 00285 $ DASUM( BLKSIZ, WORK( INDRV1+1 ), 1 ) 00286 CALL DSCAL( BLKSIZ, SCL, WORK( INDRV1+1 ), 1 ) 00287 * 00288 * Solve the system LU = Pb. 00289 * 00290 CALL DLAGTS( -1, BLKSIZ, WORK( INDRV4+1 ), WORK( INDRV2+2 ), 00291 $ WORK( INDRV3+1 ), WORK( INDRV5+1 ), IWORK, 00292 $ WORK( INDRV1+1 ), TOL, IINFO ) 00293 * 00294 * Reorthogonalize by modified Gram-Schmidt if eigenvalues are 00295 * close enough. 00296 * 00297 IF( JBLK.EQ.1 ) 00298 $ GO TO 90 00299 IF( ABS( XJ-XJM ).GT.ORTOL ) 00300 $ GPIND = J 00301 IF( GPIND.NE.J ) THEN 00302 DO 80 I = GPIND, J - 1 00303 ZTR = -DDOT( BLKSIZ, WORK( INDRV1+1 ), 1, Z( B1, I ), 00304 $ 1 ) 00305 CALL DAXPY( BLKSIZ, ZTR, Z( B1, I ), 1, 00306 $ WORK( INDRV1+1 ), 1 ) 00307 80 CONTINUE 00308 END IF 00309 * 00310 * Check the infinity norm of the iterate. 00311 * 00312 90 CONTINUE 00313 JMAX = IDAMAX( BLKSIZ, WORK( INDRV1+1 ), 1 ) 00314 NRM = ABS( WORK( INDRV1+JMAX ) ) 00315 * 00316 * Continue for additional iterations after norm reaches 00317 * stopping criterion. 00318 * 00319 IF( NRM.LT.DTPCRT ) 00320 $ GO TO 70 00321 NRMCHK = NRMCHK + 1 00322 IF( NRMCHK.LT.EXTRA+1 ) 00323 $ GO TO 70 00324 * 00325 GO TO 110 00326 * 00327 * If stopping criterion was not satisfied, update info and 00328 * store eigenvector number in array ifail. 00329 * 00330 100 CONTINUE 00331 INFO = INFO + 1 00332 IFAIL( INFO ) = J 00333 * 00334 * Accept iterate as jth eigenvector. 00335 * 00336 110 CONTINUE 00337 SCL = ONE / DNRM2( BLKSIZ, WORK( INDRV1+1 ), 1 ) 00338 JMAX = IDAMAX( BLKSIZ, WORK( INDRV1+1 ), 1 ) 00339 IF( WORK( INDRV1+JMAX ).LT.ZERO ) 00340 $ SCL = -SCL 00341 CALL DSCAL( BLKSIZ, SCL, WORK( INDRV1+1 ), 1 ) 00342 120 CONTINUE 00343 DO 130 I = 1, N 00344 Z( I, J ) = ZERO 00345 130 CONTINUE 00346 DO 140 I = 1, BLKSIZ 00347 Z( B1+I-1, J ) = WORK( INDRV1+I ) 00348 140 CONTINUE 00349 * 00350 * Save the shift to check eigenvalue spacing at next 00351 * iteration. 00352 * 00353 XJM = XJ 00354 * 00355 150 CONTINUE 00356 160 CONTINUE 00357 * 00358 RETURN 00359 * 00360 * End of DSTEIN 00361 * 00362 END
1.7.3 
