LAPACK: SRC/dstein.f Source File

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