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