LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
cstein.f
Go to the documentation of this file.
00001 *> \brief \b CSTEIN
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *> \htmlonly
00009 *> Download CSTEIN + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/cstein.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/cstein.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/cstein.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       SUBROUTINE CSTEIN( 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( * )
00031 *       COMPLEX            Z( LDZ, * )
00032 *       ..
00033 *  
00034 *
00035 *> \par Purpose:
00036 *  =============
00037 *>
00038 *> \verbatim
00039 *>
00040 *> CSTEIN 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 CUNMTR or CUPMTR 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 REAL 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 REAL 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 REAL 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 SSTEBZ 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 SSTEBZ 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 SSTEBZ is expected here. )
00110 *> \endverbatim
00111 *>
00112 *> \param[out] Z
00113 *> \verbatim
00114 *>          Z is COMPLEX 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 REAL 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 complexOTHERcomputational
00180 *
00181 *  =====================================================================
00182       SUBROUTINE CSTEIN( 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       REAL               D( * ), E( * ), W( * ), WORK( * )
00197       COMPLEX            Z( LDZ, * )
00198 *     ..
00199 *
00200 * =====================================================================
00201 *
00202 *     .. Parameters ..
00203       COMPLEX            CZERO, CONE
00204       PARAMETER          ( CZERO = ( 0.0E+0, 0.0E+0 ),
00205      $                   CONE = ( 1.0E+0, 0.0E+0 ) )
00206       REAL               ZERO, ONE, TEN, ODM3, ODM1
00207       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0, TEN = 1.0E+1,
00208      $                   ODM3 = 1.0E-3, ODM1 = 1.0E-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       REAL               CTR, EPS, EPS1, NRM, ONENRM, ORTOL, PERTOL,
00217      $                   SCL, SEP, STPCRT, TOL, XJ, XJM
00218 *     ..
00219 *     .. Local Arrays ..
00220       INTEGER            ISEED( 4 )
00221 *     ..
00222 *     .. External Functions ..
00223       INTEGER            ISAMAX
00224       REAL               SASUM, SLAMCH, SNRM2
00225       EXTERNAL           ISAMAX, SASUM, SLAMCH, SNRM2
00226 *     ..
00227 *     .. External Subroutines ..
00228       EXTERNAL           SCOPY, SLAGTF, SLAGTS, SLARNV, SSCAL, XERBLA
00229 *     ..
00230 *     .. Intrinsic Functions ..
00231       INTRINSIC          ABS, CMPLX, MAX, REAL, 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( 'CSTEIN', -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 = SLAMCH( 'Precision' )
00280 *
00281 *     Initialize seed for random number generator SLARNV.
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          STPCRT = 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 SLARNV( 2, ISEED, BLKSIZ, WORK( INDRV1+1 ) )
00361 *
00362 *           Copy the matrix T so it won't be destroyed in factorization.
00363 *
00364             CALL SCOPY( BLKSIZ, D( B1 ), 1, WORK( INDRV4+1 ), 1 )
00365             CALL SCOPY( BLKSIZ-1, E( B1 ), 1, WORK( INDRV2+2 ), 1 )
00366             CALL SCOPY( 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 SLAGTF( 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      $            SASUM( BLKSIZ, WORK( INDRV1+1 ), 1 )
00387             CALL SSCAL( BLKSIZ, SCL, WORK( INDRV1+1 ), 1 )
00388 *
00389 *           Solve the system LU = Pb.
00390 *
00391             CALL SLAGTS( -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                   CTR = ZERO
00405                   DO 80 JR = 1, BLKSIZ
00406                      CTR = CTR + WORK( INDRV1+JR )*
00407      $                     REAL( Z( B1-1+JR, I ) )
00408    80             CONTINUE
00409                   DO 90 JR = 1, BLKSIZ
00410                      WORK( INDRV1+JR ) = WORK( INDRV1+JR ) -
00411      $                                   CTR*REAL( 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 = ISAMAX( 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.STPCRT )
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 / SNRM2( BLKSIZ, WORK( INDRV1+1 ), 1 )
00444             JMAX = ISAMAX( BLKSIZ, WORK( INDRV1+1 ), 1 )
00445             IF( WORK( INDRV1+JMAX ).LT.ZERO )
00446      $         SCL = -SCL
00447             CALL SSCAL( 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 ) = CMPLX( 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 CSTEIN
00467 *
00468       END
 All Files Functions