dstein.f

*> \brief \b DSTEIN
*
*  =========== DOCUMENTATION ===========
*
* Online html documentation available at
*            http://www.netlib.org/lapack/explore-html/
*
*> \htmlonly
*> Download DSTEIN + dependencies
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dstein.f">
*> [TGZ]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dstein.f">
*> [ZIP]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dstein.f">
*> [TXT]</a>
*> \endhtmlonly
*
*  Definition:
*  ===========
*
*       SUBROUTINE DSTEIN( N, D, E, M, W, IBLOCK, ISPLIT, Z, LDZ, WORK,
*                          IWORK, IFAIL, INFO )
*
*       .. Scalar Arguments ..
*       INTEGER            INFO, LDZ, M, N
*       ..
*       .. Array Arguments ..
*       INTEGER            IBLOCK( * ), IFAIL( * ), ISPLIT( * ),
*      $                   IWORK( * )
*       DOUBLE PRECISION   D( * ), E( * ), W( * ), WORK( * ), Z( LDZ, * )
*       ..
*
*
*> \par Purpose:
*  =============
*>
*> \verbatim
*>
*> DSTEIN computes the eigenvectors of a real symmetric tridiagonal
*> matrix T corresponding to specified eigenvalues, using inverse
*> iteration.
*>
*> The maximum number of iterations allowed for each eigenvector is
*> specified by an internal parameter MAXITS (currently set to 5).
*> \endverbatim
*
*  Arguments:
*  ==========
*
*> \param[in] N
*> \verbatim
*>          N is INTEGER
*>          The order of the matrix.  N >= 0.
*> \endverbatim
*>
*> \param[in] D
*> \verbatim
*>          D is DOUBLE PRECISION array, dimension (N)
*>          The n diagonal elements of the tridiagonal matrix T.
*> \endverbatim
*>
*> \param[in] E
*> \verbatim
*>          E is DOUBLE PRECISION array, dimension (N-1)
*>          The (n-1) subdiagonal elements of the tridiagonal matrix
*>          T, in elements 1 to N-1.
*> \endverbatim
*>
*> \param[in] M
*> \verbatim
*>          M is INTEGER
*>          The number of eigenvectors to be found.  0 <= M <= N.
*> \endverbatim
*>
*> \param[in] W
*> \verbatim
*>          W is DOUBLE PRECISION array, dimension (N)
*>          The first M elements of W contain the eigenvalues for
*>          which eigenvectors are to be computed.  The eigenvalues
*>          should be grouped by split-off block and ordered from
*>          smallest to largest within the block.  ( The output array
*>          W from DSTEBZ with ORDER = 'B' is expected here. )
*> \endverbatim
*>
*> \param[in] IBLOCK
*> \verbatim
*>          IBLOCK is INTEGER array, dimension (N)
*>          The submatrix indices associated with the corresponding
*>          eigenvalues in W; IBLOCK(i)=1 if eigenvalue W(i) belongs to
*>          the first submatrix from the top, =2 if W(i) belongs to
*>          the second submatrix, etc.  ( The output array IBLOCK
*>          from DSTEBZ is expected here. )
*> \endverbatim
*>
*> \param[in] ISPLIT
*> \verbatim
*>          ISPLIT is INTEGER array, dimension (N)
*>          The splitting points, at which T breaks up into submatrices.
*>          The first submatrix consists of rows/columns 1 to
*>          ISPLIT( 1 ), the second of rows/columns ISPLIT( 1 )+1
*>          through ISPLIT( 2 ), etc.
*>          ( The output array ISPLIT from DSTEBZ is expected here. )
*> \endverbatim
*>
*> \param[out] Z
*> \verbatim
*>          Z is DOUBLE PRECISION array, dimension (LDZ, M)
*>          The computed eigenvectors.  The eigenvector associated
*>          with the eigenvalue W(i) is stored in the i-th column of
*>          Z.  Any vector which fails to converge is set to its current
*>          iterate after MAXITS iterations.
*> \endverbatim
*>
*> \param[in] LDZ
*> \verbatim
*>          LDZ is INTEGER
*>          The leading dimension of the array Z.  LDZ >= max(1,N).
*> \endverbatim
*>
*> \param[out] WORK
*> \verbatim
*>          WORK is DOUBLE PRECISION array, dimension (5*N)
*> \endverbatim
*>
*> \param[out] IWORK
*> \verbatim
*>          IWORK is INTEGER array, dimension (N)
*> \endverbatim
*>
*> \param[out] IFAIL
*> \verbatim
*>          IFAIL is INTEGER array, dimension (M)
*>          On normal exit, all elements of IFAIL are zero.
*>          If one or more eigenvectors fail to converge after
*>          MAXITS iterations, then their indices are stored in
*>          array IFAIL.
*> \endverbatim
*>
*> \param[out] INFO
*> \verbatim
*>          INFO is INTEGER
*>          = 0: successful exit.
*>          < 0: if INFO = -i, the i-th argument had an illegal value
*>          > 0: if INFO = i, then i eigenvectors failed to converge
*>               in MAXITS iterations.  Their indices are stored in
*>               array IFAIL.
*> \endverbatim
*
*> \par Internal Parameters:
*  =========================
*>
*> \verbatim
*>  MAXITS  INTEGER, default = 5
*>          The maximum number of iterations performed.
*>
*>  EXTRA   INTEGER, default = 2
*>          The number of iterations performed after norm growth
*>          criterion is satisfied, should be at least 1.
*> \endverbatim
*
*  Authors:
*  ========
*
*> \author Univ. of Tennessee
*> \author Univ. of California Berkeley
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
*> \ingroup stein
*
*  =====================================================================
      SUBROUTINE dstein( N, D, E, M, W, IBLOCK, ISPLIT, Z, LDZ, WORK,
     $                   IWORK, IFAIL, INFO )
*
*  -- LAPACK computational routine --
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*
*     .. Scalar Arguments ..
      INTEGER            INFO, LDZ, M, N
*     ..
*     .. Array Arguments ..
      INTEGER            IBLOCK( * ), IFAIL( * ), ISPLIT( * ),
     $                   iwork( * )
      DOUBLE PRECISION   D( * ), E( * ), W( * ), WORK( * ), Z( ldz, * )
*     ..
*
*  =====================================================================
*
*     .. Parameters ..
      DOUBLE PRECISION   ZERO, ONE, TEN, ODM3, ODM1
      parameter( zero = 0.0d+0, one = 1.0d+0, ten = 1.0d+1,
     $                   odm3 = 1.0d-3, odm1 = 1.0d-1 )
      INTEGER            MAXITS, EXTRA
      parameter( maxits = 5, extra = 2 )
*     ..
*     .. Local Scalars ..
      INTEGER            B1, BLKSIZ, BN, GPIND, I, IINFO, INDRV1,
     $                   indrv2, indrv3, indrv4, indrv5, its, j, j1,
     $                   jblk, jmax, nblk, nrmchk
      DOUBLE PRECISION   DTPCRT, EPS, EPS1, NRM, ONENRM, ORTOL, PERTOL,
     $                   scl, sep, tol, xj, xjm, ztr
*     ..
*     .. Local Arrays ..
      INTEGER            ISEED( 4 )
*     ..
*     .. External Functions ..
      INTEGER            IDAMAX
      DOUBLE PRECISION   DDOT, DLAMCH, DNRM2
      EXTERNAL           idamax, ddot, dlamch, dnrm2
*     ..
*     .. External Subroutines ..
      EXTERNAL           daxpy, dcopy, dlagtf, dlagts, dlarnv, dscal,
     $                   xerbla
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          abs, max, sqrt
*     ..
*     .. Executable Statements ..
*
*     Test the input parameters.
*
      info = 0
      DO 10 i = 1, m
         ifail( i ) = 0
   10 CONTINUE
*
      IF( n.LT.0 ) THEN
         info = -1
      ELSE IF( m.LT.0 .OR. m.GT.n ) THEN
         info = -4
      ELSE IF( ldz.LT.max( 1, n ) ) THEN
         info = -9
      ELSE
         DO 20 j = 2, m
            IF( iblock( j ).LT.iblock( j-1 ) ) THEN
               info = -6
               GO TO 30
            END IF
            IF( iblock( j ).EQ.iblock( j-1 ) .AND. w( j ).LT.w( j-1 ) )
     $           THEN
               info = -5
               GO TO 30
            END IF
   20    CONTINUE
   30    CONTINUE
      END IF
*
      IF( info.NE.0 ) THEN
         CALL xerbla( 'DSTEIN', -info )
         RETURN
      END IF
*
*     Quick return if possible
*
      IF( n.EQ.0 .OR. m.EQ.0 ) THEN
         RETURN
      ELSE IF( n.EQ.1 ) THEN
         z( 1, 1 ) = one
         RETURN
      END IF
*
*     Get machine constants.
*
      eps = dlamch( 'Precision' )
*
*     Initialize seed for random number generator DLARNV.
*
      DO 40 i = 1, 4
         iseed( i ) = 1
   40 CONTINUE
*
*     Initialize pointers.
*
      indrv1 = 0
      indrv2 = indrv1 + n
      indrv3 = indrv2 + n
      indrv4 = indrv3 + n
      indrv5 = indrv4 + n
*
*     Compute eigenvectors of matrix blocks.
*
      j1 = 1
      DO 160 nblk = 1, iblock( m )
*
*        Find starting and ending indices of block nblk.
*
         IF( nblk.EQ.1 ) THEN
            b1 = 1
         ELSE
            b1 = isplit( nblk-1 ) + 1
         END IF
         bn = isplit( nblk )
         blksiz = bn - b1 + 1
         IF( blksiz.EQ.1 )
     $      GO TO 60
         gpind = j1
*
*        Compute reorthogonalization criterion and stopping criterion.
*
         onenrm = abs( d( b1 ) ) + abs( e( b1 ) )
         onenrm = max( onenrm, abs( d( bn ) )+abs( e( bn-1 ) ) )
         DO 50 i = b1 + 1, bn - 1
            onenrm = max( onenrm, abs( d( i ) )+abs( e( i-1 ) )+
     $               abs( e( i ) ) )
   50    CONTINUE
         ortol = odm3*onenrm
*
         dtpcrt = sqrt( odm1 / blksiz )
*
*        Loop through eigenvalues of block nblk.
*
   60    CONTINUE
         jblk = 0
         DO 150 j = j1, m
            IF( iblock( j ).NE.nblk ) THEN
               j1 = j
               GO TO 160
            END IF
            jblk = jblk + 1
            xj = w( j )
*
*           Skip all the work if the block size is one.
*
            IF( blksiz.EQ.1 ) THEN
               work( indrv1+1 ) = one
               GO TO 120
            END IF
*
*           If eigenvalues j and j-1 are too close, add a relatively
*           small perturbation.
*
            IF( jblk.GT.1 ) THEN
               eps1 = abs( eps*xj )
               pertol = ten*eps1
               sep = xj - xjm
               IF( sep.LT.pertol )
     $            xj = xjm + pertol
            END IF
*
            its = 0
            nrmchk = 0
*
*           Get random starting vector.
*
            CALL dlarnv( 2, iseed, blksiz, work( indrv1+1 ) )
*
*           Copy the matrix T so it won't be destroyed in factorization.
*
            CALL dcopy( blksiz, d( b1 ), 1, work( indrv4+1 ), 1 )
            CALL dcopy( blksiz-1, e( b1 ), 1, work( indrv2+2 ), 1 )
            CALL dcopy( blksiz-1, e( b1 ), 1, work( indrv3+1 ), 1 )
*
*           Compute LU factors with partial pivoting  ( PT = LU )
*
            tol = zero
            CALL dlagtf( blksiz, work( indrv4+1 ), xj, work( indrv2+2 ),
     $                   work( indrv3+1 ), tol, work( indrv5+1 ), iwork,
     $                   iinfo )
*
*           Update iteration count.
*
   70       CONTINUE
            its = its + 1
            IF( its.GT.maxits )
     $         GO TO 100
*
*           Normalize and scale the righthand side vector Pb.
*
            jmax = idamax( blksiz, work( indrv1+1 ), 1 )
            scl = blksiz*onenrm*max( eps,
     $            abs( work( indrv4+blksiz ) ) ) /
     $            abs( work( indrv1+jmax ) )
            CALL dscal( blksiz, scl, work( indrv1+1 ), 1 )
*
*           Solve the system LU = Pb.
*
            CALL dlagts( -1, blksiz, work( indrv4+1 ), work( indrv2+2 ),
     $                   work( indrv3+1 ), work( indrv5+1 ), iwork,
     $                   work( indrv1+1 ), tol, iinfo )
*
*           Reorthogonalize by modified Gram-Schmidt if eigenvalues are
*           close enough.
*
            IF( jblk.EQ.1 )
     $         GO TO 90
            IF( abs( xj-xjm ).GT.ortol )
     $         gpind = j
            IF( gpind.NE.j ) THEN
               DO 80 i = gpind, j - 1
                  ztr = -ddot( blksiz, work( indrv1+1 ), 1, z( b1, i ),
     $                  1 )
                  CALL daxpy( blksiz, ztr, z( b1, i ), 1,
     $                        work( indrv1+1 ), 1 )
   80          CONTINUE
            END IF
*
*           Check the infinity norm of the iterate.
*
   90       CONTINUE
            jmax = idamax( blksiz, work( indrv1+1 ), 1 )
            nrm = abs( work( indrv1+jmax ) )
*
*           Continue for additional iterations after norm reaches
*           stopping criterion.
*
            IF( nrm.LT.dtpcrt )
     $         GO TO 70
            nrmchk = nrmchk + 1
            IF( nrmchk.LT.extra+1 )
     $         GO TO 70
*
            GO TO 110
*
*           If stopping criterion was not satisfied, update info and
*           store eigenvector number in array ifail.
*
  100       CONTINUE
            info = info + 1
            ifail( info ) = j
*
*           Accept iterate as jth eigenvector.
*
  110       CONTINUE
            scl = one / dnrm2( blksiz, work( indrv1+1 ), 1 )
            jmax = idamax( blksiz, work( indrv1+1 ), 1 )
            IF( work( indrv1+jmax ).LT.zero )
     $         scl = -scl
            CALL dscal( blksiz, scl, work( indrv1+1 ), 1 )
  120       CONTINUE
            DO 130 i = 1, n
               z( i, j ) = zero
  130       CONTINUE
            DO 140 i = 1, blksiz
               z( b1+i-1, j ) = work( indrv1+i )
  140       CONTINUE
*
*           Save the shift to check eigenvalue spacing at next
*           iteration.
*
            xjm = xj
*
  150    CONTINUE
  160 CONTINUE
*
      RETURN
*
*     End of DSTEIN
*
      END