LAPACK: SRC/cstein.f Source File

Go to the documentation of this file.
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