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