![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b SSTEIN 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download SSTEIN + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/sstein.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/sstein.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/sstein.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE SSTEIN( N, D, E, M, W, IBLOCK, ISPLIT, Z, LDZ, WORK, 00022 * IWORK, IFAIL, INFO ) 00023 * 00024 * .. Scalar Arguments .. 00025 * INTEGER INFO, LDZ, M, N 00026 * .. 00027 * .. Array Arguments .. 00028 * INTEGER IBLOCK( * ), IFAIL( * ), ISPLIT( * ), 00029 * $ IWORK( * ) 00030 * REAL D( * ), E( * ), W( * ), WORK( * ), Z( LDZ, * ) 00031 * .. 00032 * 00033 * 00034 *> \par Purpose: 00035 * ============= 00036 *> 00037 *> \verbatim 00038 *> 00039 *> SSTEIN computes the eigenvectors of a real symmetric tridiagonal 00040 *> matrix T corresponding to specified eigenvalues, using inverse 00041 *> iteration. 00042 *> 00043 *> The maximum number of iterations allowed for each eigenvector is 00044 *> specified by an internal parameter MAXITS (currently set to 5). 00045 *> \endverbatim 00046 * 00047 * Arguments: 00048 * ========== 00049 * 00050 *> \param[in] N 00051 *> \verbatim 00052 *> N is INTEGER 00053 *> The order of the matrix. N >= 0. 00054 *> \endverbatim 00055 *> 00056 *> \param[in] D 00057 *> \verbatim 00058 *> D is REAL array, dimension (N) 00059 *> The n diagonal elements of the tridiagonal matrix T. 00060 *> \endverbatim 00061 *> 00062 *> \param[in] E 00063 *> \verbatim 00064 *> E is REAL array, dimension (N-1) 00065 *> The (n-1) subdiagonal elements of the tridiagonal matrix 00066 *> T, in elements 1 to N-1. 00067 *> \endverbatim 00068 *> 00069 *> \param[in] M 00070 *> \verbatim 00071 *> M is INTEGER 00072 *> The number of eigenvectors to be found. 0 <= M <= N. 00073 *> \endverbatim 00074 *> 00075 *> \param[in] W 00076 *> \verbatim 00077 *> W is REAL array, dimension (N) 00078 *> The first M elements of W contain the eigenvalues for 00079 *> which eigenvectors are to be computed. The eigenvalues 00080 *> should be grouped by split-off block and ordered from 00081 *> smallest to largest within the block. ( The output array 00082 *> W from SSTEBZ with ORDER = 'B' is expected here. ) 00083 *> \endverbatim 00084 *> 00085 *> \param[in] IBLOCK 00086 *> \verbatim 00087 *> IBLOCK is INTEGER array, dimension (N) 00088 *> The submatrix indices associated with the corresponding 00089 *> eigenvalues in W; IBLOCK(i)=1 if eigenvalue W(i) belongs to 00090 *> the first submatrix from the top, =2 if W(i) belongs to 00091 *> the second submatrix, etc. ( The output array IBLOCK 00092 *> from SSTEBZ is expected here. ) 00093 *> \endverbatim 00094 *> 00095 *> \param[in] ISPLIT 00096 *> \verbatim 00097 *> ISPLIT is INTEGER array, dimension (N) 00098 *> The splitting points, at which T breaks up into submatrices. 00099 *> The first submatrix consists of rows/columns 1 to 00100 *> ISPLIT( 1 ), the second of rows/columns ISPLIT( 1 )+1 00101 *> through ISPLIT( 2 ), etc. 00102 *> ( The output array ISPLIT from SSTEBZ is expected here. ) 00103 *> \endverbatim 00104 *> 00105 *> \param[out] Z 00106 *> \verbatim 00107 *> Z is REAL array, dimension (LDZ, M) 00108 *> The computed eigenvectors. The eigenvector associated 00109 *> with the eigenvalue W(i) is stored in the i-th column of 00110 *> Z. Any vector which fails to converge is set to its current 00111 *> iterate after MAXITS iterations. 00112 *> \endverbatim 00113 *> 00114 *> \param[in] LDZ 00115 *> \verbatim 00116 *> LDZ is INTEGER 00117 *> The leading dimension of the array Z. LDZ >= max(1,N). 00118 *> \endverbatim 00119 *> 00120 *> \param[out] WORK 00121 *> \verbatim 00122 *> WORK is REAL array, dimension (5*N) 00123 *> \endverbatim 00124 *> 00125 *> \param[out] IWORK 00126 *> \verbatim 00127 *> IWORK is INTEGER array, dimension (N) 00128 *> \endverbatim 00129 *> 00130 *> \param[out] IFAIL 00131 *> \verbatim 00132 *> IFAIL is INTEGER array, dimension (M) 00133 *> On normal exit, all elements of IFAIL are zero. 00134 *> If one or more eigenvectors fail to converge after 00135 *> MAXITS iterations, then their indices are stored in 00136 *> array IFAIL. 00137 *> \endverbatim 00138 *> 00139 *> \param[out] INFO 00140 *> \verbatim 00141 *> INFO is INTEGER 00142 *> = 0: successful exit. 00143 *> < 0: if INFO = -i, the i-th argument had an illegal value 00144 *> > 0: if INFO = i, then i eigenvectors failed to converge 00145 *> in MAXITS iterations. Their indices are stored in 00146 *> array IFAIL. 00147 *> \endverbatim 00148 * 00149 *> \par Internal Parameters: 00150 * ========================= 00151 *> 00152 *> \verbatim 00153 *> MAXITS INTEGER, default = 5 00154 *> The maximum number of iterations performed. 00155 *> 00156 *> EXTRA INTEGER, default = 2 00157 *> The number of iterations performed after norm growth 00158 *> criterion is satisfied, should be at least 1. 00159 *> \endverbatim 00160 * 00161 * Authors: 00162 * ======== 00163 * 00164 *> \author Univ. of Tennessee 00165 *> \author Univ. of California Berkeley 00166 *> \author Univ. of Colorado Denver 00167 *> \author NAG Ltd. 00168 * 00169 *> \date November 2011 00170 * 00171 *> \ingroup realOTHERcomputational 00172 * 00173 * ===================================================================== 00174 SUBROUTINE SSTEIN( N, D, E, M, W, IBLOCK, ISPLIT, Z, LDZ, WORK, 00175 $ IWORK, IFAIL, INFO ) 00176 * 00177 * -- LAPACK computational routine (version 3.4.0) -- 00178 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00179 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00180 * November 2011 00181 * 00182 * .. Scalar Arguments .. 00183 INTEGER INFO, LDZ, M, N 00184 * .. 00185 * .. Array Arguments .. 00186 INTEGER IBLOCK( * ), IFAIL( * ), ISPLIT( * ), 00187 $ IWORK( * ) 00188 REAL D( * ), E( * ), W( * ), WORK( * ), Z( LDZ, * ) 00189 * .. 00190 * 00191 * ===================================================================== 00192 * 00193 * .. Parameters .. 00194 REAL ZERO, ONE, TEN, ODM3, ODM1 00195 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0, TEN = 1.0E+1, 00196 $ ODM3 = 1.0E-3, ODM1 = 1.0E-1 ) 00197 INTEGER MAXITS, EXTRA 00198 PARAMETER ( MAXITS = 5, EXTRA = 2 ) 00199 * .. 00200 * .. Local Scalars .. 00201 INTEGER B1, BLKSIZ, BN, GPIND, I, IINFO, INDRV1, 00202 $ INDRV2, INDRV3, INDRV4, INDRV5, ITS, J, J1, 00203 $ JBLK, JMAX, NBLK, NRMCHK 00204 REAL CTR, EPS, EPS1, NRM, ONENRM, ORTOL, PERTOL, 00205 $ SCL, SEP, STPCRT, TOL, XJ, XJM 00206 * .. 00207 * .. Local Arrays .. 00208 INTEGER ISEED( 4 ) 00209 * .. 00210 * .. External Functions .. 00211 INTEGER ISAMAX 00212 REAL SASUM, SDOT, SLAMCH, SNRM2 00213 EXTERNAL ISAMAX, SASUM, SDOT, SLAMCH, SNRM2 00214 * .. 00215 * .. External Subroutines .. 00216 EXTERNAL SAXPY, SCOPY, SLAGTF, SLAGTS, SLARNV, SSCAL, 00217 $ XERBLA 00218 * .. 00219 * .. Intrinsic Functions .. 00220 INTRINSIC ABS, MAX, SQRT 00221 * .. 00222 * .. Executable Statements .. 00223 * 00224 * Test the input parameters. 00225 * 00226 INFO = 0 00227 DO 10 I = 1, M 00228 IFAIL( I ) = 0 00229 10 CONTINUE 00230 * 00231 IF( N.LT.0 ) THEN 00232 INFO = -1 00233 ELSE IF( M.LT.0 .OR. M.GT.N ) THEN 00234 INFO = -4 00235 ELSE IF( LDZ.LT.MAX( 1, N ) ) THEN 00236 INFO = -9 00237 ELSE 00238 DO 20 J = 2, M 00239 IF( IBLOCK( J ).LT.IBLOCK( J-1 ) ) THEN 00240 INFO = -6 00241 GO TO 30 00242 END IF 00243 IF( IBLOCK( J ).EQ.IBLOCK( J-1 ) .AND. W( J ).LT.W( J-1 ) ) 00244 $ THEN 00245 INFO = -5 00246 GO TO 30 00247 END IF 00248 20 CONTINUE 00249 30 CONTINUE 00250 END IF 00251 * 00252 IF( INFO.NE.0 ) THEN 00253 CALL XERBLA( 'SSTEIN', -INFO ) 00254 RETURN 00255 END IF 00256 * 00257 * Quick return if possible 00258 * 00259 IF( N.EQ.0 .OR. M.EQ.0 ) THEN 00260 RETURN 00261 ELSE IF( N.EQ.1 ) THEN 00262 Z( 1, 1 ) = ONE 00263 RETURN 00264 END IF 00265 * 00266 * Get machine constants. 00267 * 00268 EPS = SLAMCH( 'Precision' ) 00269 * 00270 * Initialize seed for random number generator SLARNV. 00271 * 00272 DO 40 I = 1, 4 00273 ISEED( I ) = 1 00274 40 CONTINUE 00275 * 00276 * Initialize pointers. 00277 * 00278 INDRV1 = 0 00279 INDRV2 = INDRV1 + N 00280 INDRV3 = INDRV2 + N 00281 INDRV4 = INDRV3 + N 00282 INDRV5 = INDRV4 + N 00283 * 00284 * Compute eigenvectors of matrix blocks. 00285 * 00286 J1 = 1 00287 DO 160 NBLK = 1, IBLOCK( M ) 00288 * 00289 * Find starting and ending indices of block nblk. 00290 * 00291 IF( NBLK.EQ.1 ) THEN 00292 B1 = 1 00293 ELSE 00294 B1 = ISPLIT( NBLK-1 ) + 1 00295 END IF 00296 BN = ISPLIT( NBLK ) 00297 BLKSIZ = BN - B1 + 1 00298 IF( BLKSIZ.EQ.1 ) 00299 $ GO TO 60 00300 GPIND = B1 00301 * 00302 * Compute reorthogonalization criterion and stopping criterion. 00303 * 00304 ONENRM = ABS( D( B1 ) ) + ABS( E( B1 ) ) 00305 ONENRM = MAX( ONENRM, ABS( D( BN ) )+ABS( E( BN-1 ) ) ) 00306 DO 50 I = B1 + 1, BN - 1 00307 ONENRM = MAX( ONENRM, ABS( D( I ) )+ABS( E( I-1 ) )+ 00308 $ ABS( E( I ) ) ) 00309 50 CONTINUE 00310 ORTOL = ODM3*ONENRM 00311 * 00312 STPCRT = SQRT( ODM1 / BLKSIZ ) 00313 * 00314 * Loop through eigenvalues of block nblk. 00315 * 00316 60 CONTINUE 00317 JBLK = 0 00318 DO 150 J = J1, M 00319 IF( IBLOCK( J ).NE.NBLK ) THEN 00320 J1 = J 00321 GO TO 160 00322 END IF 00323 JBLK = JBLK + 1 00324 XJ = W( J ) 00325 * 00326 * Skip all the work if the block size is one. 00327 * 00328 IF( BLKSIZ.EQ.1 ) THEN 00329 WORK( INDRV1+1 ) = ONE 00330 GO TO 120 00331 END IF 00332 * 00333 * If eigenvalues j and j-1 are too close, add a relatively 00334 * small perturbation. 00335 * 00336 IF( JBLK.GT.1 ) THEN 00337 EPS1 = ABS( EPS*XJ ) 00338 PERTOL = TEN*EPS1 00339 SEP = XJ - XJM 00340 IF( SEP.LT.PERTOL ) 00341 $ XJ = XJM + PERTOL 00342 END IF 00343 * 00344 ITS = 0 00345 NRMCHK = 0 00346 * 00347 * Get random starting vector. 00348 * 00349 CALL SLARNV( 2, ISEED, BLKSIZ, WORK( INDRV1+1 ) ) 00350 * 00351 * Copy the matrix T so it won't be destroyed in factorization. 00352 * 00353 CALL SCOPY( BLKSIZ, D( B1 ), 1, WORK( INDRV4+1 ), 1 ) 00354 CALL SCOPY( BLKSIZ-1, E( B1 ), 1, WORK( INDRV2+2 ), 1 ) 00355 CALL SCOPY( BLKSIZ-1, E( B1 ), 1, WORK( INDRV3+1 ), 1 ) 00356 * 00357 * Compute LU factors with partial pivoting ( PT = LU ) 00358 * 00359 TOL = ZERO 00360 CALL SLAGTF( BLKSIZ, WORK( INDRV4+1 ), XJ, WORK( INDRV2+2 ), 00361 $ WORK( INDRV3+1 ), TOL, WORK( INDRV5+1 ), IWORK, 00362 $ IINFO ) 00363 * 00364 * Update iteration count. 00365 * 00366 70 CONTINUE 00367 ITS = ITS + 1 00368 IF( ITS.GT.MAXITS ) 00369 $ GO TO 100 00370 * 00371 * Normalize and scale the righthand side vector Pb. 00372 * 00373 SCL = BLKSIZ*ONENRM*MAX( EPS, 00374 $ ABS( WORK( INDRV4+BLKSIZ ) ) ) / 00375 $ SASUM( BLKSIZ, WORK( INDRV1+1 ), 1 ) 00376 CALL SSCAL( BLKSIZ, SCL, WORK( INDRV1+1 ), 1 ) 00377 * 00378 * Solve the system LU = Pb. 00379 * 00380 CALL SLAGTS( -1, BLKSIZ, WORK( INDRV4+1 ), WORK( INDRV2+2 ), 00381 $ WORK( INDRV3+1 ), WORK( INDRV5+1 ), IWORK, 00382 $ WORK( INDRV1+1 ), TOL, IINFO ) 00383 * 00384 * Reorthogonalize by modified Gram-Schmidt if eigenvalues are 00385 * close enough. 00386 * 00387 IF( JBLK.EQ.1 ) 00388 $ GO TO 90 00389 IF( ABS( XJ-XJM ).GT.ORTOL ) 00390 $ GPIND = J 00391 IF( GPIND.NE.J ) THEN 00392 DO 80 I = GPIND, J - 1 00393 CTR = -SDOT( BLKSIZ, WORK( INDRV1+1 ), 1, Z( B1, I ), 00394 $ 1 ) 00395 CALL SAXPY( BLKSIZ, CTR, Z( B1, I ), 1, 00396 $ WORK( INDRV1+1 ), 1 ) 00397 80 CONTINUE 00398 END IF 00399 * 00400 * Check the infinity norm of the iterate. 00401 * 00402 90 CONTINUE 00403 JMAX = ISAMAX( BLKSIZ, WORK( INDRV1+1 ), 1 ) 00404 NRM = ABS( WORK( INDRV1+JMAX ) ) 00405 * 00406 * Continue for additional iterations after norm reaches 00407 * stopping criterion. 00408 * 00409 IF( NRM.LT.STPCRT ) 00410 $ GO TO 70 00411 NRMCHK = NRMCHK + 1 00412 IF( NRMCHK.LT.EXTRA+1 ) 00413 $ GO TO 70 00414 * 00415 GO TO 110 00416 * 00417 * If stopping criterion was not satisfied, update info and 00418 * store eigenvector number in array ifail. 00419 * 00420 100 CONTINUE 00421 INFO = INFO + 1 00422 IFAIL( INFO ) = J 00423 * 00424 * Accept iterate as jth eigenvector. 00425 * 00426 110 CONTINUE 00427 SCL = ONE / SNRM2( BLKSIZ, WORK( INDRV1+1 ), 1 ) 00428 JMAX = ISAMAX( BLKSIZ, WORK( INDRV1+1 ), 1 ) 00429 IF( WORK( INDRV1+JMAX ).LT.ZERO ) 00430 $ SCL = -SCL 00431 CALL SSCAL( BLKSIZ, SCL, WORK( INDRV1+1 ), 1 ) 00432 120 CONTINUE 00433 DO 130 I = 1, N 00434 Z( I, J ) = ZERO 00435 130 CONTINUE 00436 DO 140 I = 1, BLKSIZ 00437 Z( B1+I-1, J ) = WORK( INDRV1+I ) 00438 140 CONTINUE 00439 * 00440 * Save the shift to check eigenvalue spacing at next 00441 * iteration. 00442 * 00443 XJM = XJ 00444 * 00445 150 CONTINUE 00446 160 CONTINUE 00447 * 00448 RETURN 00449 * 00450 * End of SSTEIN 00451 * 00452 END