LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
dstein.f
Go to the documentation of this file.
00001 *> \brief \b DSTEIN
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *> \htmlonly
00009 *> Download DSTEIN + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dstein.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dstein.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dstein.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       SUBROUTINE DSTEIN( 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( * ), Z( LDZ, * )
00031 *       ..
00032 *  
00033 *
00034 *> \par Purpose:
00035 *  =============
00036 *>
00037 *> \verbatim
00038 *>
00039 *> DSTEIN 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DSTEBZ 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 DSTEBZ 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 DSTEBZ is expected here. )
00103 *> \endverbatim
00104 *>
00105 *> \param[out] Z
00106 *> \verbatim
00107 *>          Z is DOUBLE PRECISION 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 DOUBLE PRECISION 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 doubleOTHERcomputational
00172 *
00173 *  =====================================================================
00174       SUBROUTINE DSTEIN( 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       DOUBLE PRECISION   D( * ), E( * ), W( * ), WORK( * ), Z( LDZ, * )
00189 *     ..
00190 *
00191 *  =====================================================================
00192 *
00193 *     .. Parameters ..
00194       DOUBLE PRECISION   ZERO, ONE, TEN, ODM3, ODM1
00195       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0, TEN = 1.0D+1,
00196      $                   ODM3 = 1.0D-3, ODM1 = 1.0D-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       DOUBLE PRECISION   DTPCRT, EPS, EPS1, NRM, ONENRM, ORTOL, PERTOL,
00205      $                   SCL, SEP, TOL, XJ, XJM, ZTR
00206 *     ..
00207 *     .. Local Arrays ..
00208       INTEGER            ISEED( 4 )
00209 *     ..
00210 *     .. External Functions ..
00211       INTEGER            IDAMAX
00212       DOUBLE PRECISION   DASUM, DDOT, DLAMCH, DNRM2
00213       EXTERNAL           IDAMAX, DASUM, DDOT, DLAMCH, DNRM2
00214 *     ..
00215 *     .. External Subroutines ..
00216       EXTERNAL           DAXPY, DCOPY, DLAGTF, DLAGTS, DLARNV, DSCAL,
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( 'DSTEIN', -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 = DLAMCH( 'Precision' )
00269 *
00270 *     Initialize seed for random number generator DLARNV.
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          DTPCRT = 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 DLARNV( 2, ISEED, BLKSIZ, WORK( INDRV1+1 ) )
00350 *
00351 *           Copy the matrix T so it won't be destroyed in factorization.
00352 *
00353             CALL DCOPY( BLKSIZ, D( B1 ), 1, WORK( INDRV4+1 ), 1 )
00354             CALL DCOPY( BLKSIZ-1, E( B1 ), 1, WORK( INDRV2+2 ), 1 )
00355             CALL DCOPY( 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 DLAGTF( 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      $            DASUM( BLKSIZ, WORK( INDRV1+1 ), 1 )
00376             CALL DSCAL( BLKSIZ, SCL, WORK( INDRV1+1 ), 1 )
00377 *
00378 *           Solve the system LU = Pb.
00379 *
00380             CALL DLAGTS( -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                   ZTR = -DDOT( BLKSIZ, WORK( INDRV1+1 ), 1, Z( B1, I ),
00394      $                  1 )
00395                   CALL DAXPY( BLKSIZ, ZTR, 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 = IDAMAX( 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.DTPCRT )
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 / DNRM2( BLKSIZ, WORK( INDRV1+1 ), 1 )
00428             JMAX = IDAMAX( BLKSIZ, WORK( INDRV1+1 ), 1 )
00429             IF( WORK( INDRV1+JMAX ).LT.ZERO )
00430      $         SCL = -SCL
00431             CALL DSCAL( 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 DSTEIN
00451 *
00452       END
 All Files Functions