Return to zlaein.f CVS log

Up to [local] / rpl / lapack / lapack

Commit initial.

    1:       SUBROUTINE ZLAEIN( RIGHTV, NOINIT, N, H, LDH, W, V, B, LDB, RWORK,
    2:      $                   EPS3, SMLNUM, INFO )
    3: *
    4: *  -- LAPACK auxiliary routine (version 3.2) --
    5: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
    6: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
    7: *     November 2006
    8: *
    9: *     .. Scalar Arguments ..
   10:       LOGICAL            NOINIT, RIGHTV
   11:       INTEGER            INFO, LDB, LDH, N
   12:       DOUBLE PRECISION   EPS3, SMLNUM
   13:       COMPLEX*16         W
   14: *     ..
   15: *     .. Array Arguments ..
   16:       DOUBLE PRECISION   RWORK( * )
   17:       COMPLEX*16         B( LDB, * ), H( LDH, * ), V( * )
   18: *     ..
   19: *
   20: *  Purpose
   21: *  =======
   22: *
   23: *  ZLAEIN uses inverse iteration to find a right or left eigenvector
   24: *  corresponding to the eigenvalue W of a complex upper Hessenberg
   25: *  matrix H.
   26: *
   27: *  Arguments
   28: *  =========
   29: *
   30: *  RIGHTV   (input) LOGICAL
   31: *          = .TRUE. : compute right eigenvector;
   32: *          = .FALSE.: compute left eigenvector.
   33: *
   34: *  NOINIT   (input) LOGICAL
   35: *          = .TRUE. : no initial vector supplied in V
   36: *          = .FALSE.: initial vector supplied in V.
   37: *
   38: *  N       (input) INTEGER
   39: *          The order of the matrix H.  N >= 0.
   40: *
   41: *  H       (input) COMPLEX*16 array, dimension (LDH,N)
   42: *          The upper Hessenberg matrix H.
   43: *
   44: *  LDH     (input) INTEGER
   45: *          The leading dimension of the array H.  LDH >= max(1,N).
   46: *
   47: *  W       (input) COMPLEX*16
   48: *          The eigenvalue of H whose corresponding right or left
   49: *          eigenvector is to be computed.
   50: *
   51: *  V       (input/output) COMPLEX*16 array, dimension (N)
   52: *          On entry, if NOINIT = .FALSE., V must contain a starting
   53: *          vector for inverse iteration; otherwise V need not be set.
   54: *          On exit, V contains the computed eigenvector, normalized so
   55: *          that the component of largest magnitude has magnitude 1; here
   56: *          the magnitude of a complex number (x,y) is taken to be
   57: *          |x| + |y|.
   58: *
   59: *  B       (workspace) COMPLEX*16 array, dimension (LDB,N)
   60: *
   61: *  LDB     (input) INTEGER
   62: *          The leading dimension of the array B.  LDB >= max(1,N).
   63: *
   64: *  RWORK   (workspace) DOUBLE PRECISION array, dimension (N)
   65: *
   66: *  EPS3    (input) DOUBLE PRECISION
   67: *          A small machine-dependent value which is used to perturb
   68: *          close eigenvalues, and to replace zero pivots.
   69: *
   70: *  SMLNUM  (input) DOUBLE PRECISION
   71: *          A machine-dependent value close to the underflow threshold.
   72: *
   73: *  INFO    (output) INTEGER
   74: *          = 0:  successful exit
   75: *          = 1:  inverse iteration did not converge; V is set to the
   76: *                last iterate.
   77: *
   78: *  =====================================================================
   79: *
   80: *     .. Parameters ..
   81:       DOUBLE PRECISION   ONE, TENTH
   82:       PARAMETER          ( ONE = 1.0D+0, TENTH = 1.0D-1 )
   83:       COMPLEX*16         ZERO
   84:       PARAMETER          ( ZERO = ( 0.0D+0, 0.0D+0 ) )
   85: *     ..
   86: *     .. Local Scalars ..
   87:       CHARACTER          NORMIN, TRANS
   88:       INTEGER            I, IERR, ITS, J
   89:       DOUBLE PRECISION   GROWTO, NRMSML, ROOTN, RTEMP, SCALE, VNORM
   90:       COMPLEX*16         CDUM, EI, EJ, TEMP, X
   91: *     ..
   92: *     .. External Functions ..
   93:       INTEGER            IZAMAX
   94:       DOUBLE PRECISION   DZASUM, DZNRM2
   95:       COMPLEX*16         ZLADIV
   96:       EXTERNAL           IZAMAX, DZASUM, DZNRM2, ZLADIV
   97: *     ..
   98: *     .. External Subroutines ..
   99:       EXTERNAL           ZDSCAL, ZLATRS
  100: *     ..
  101: *     .. Intrinsic Functions ..
  102:       INTRINSIC          ABS, DBLE, DIMAG, MAX, SQRT
  103: *     ..
  104: *     .. Statement Functions ..
  105:       DOUBLE PRECISION   CABS1
  106: *     ..
  107: *     .. Statement Function definitions ..
  108:       CABS1( CDUM ) = ABS( DBLE( CDUM ) ) + ABS( DIMAG( CDUM ) )
  109: *     ..
  110: *     .. Executable Statements ..
  111: *
  112:       INFO = 0
  113: *
  114: *     GROWTO is the threshold used in the acceptance test for an
  115: *     eigenvector.
  116: *
  117:       ROOTN = SQRT( DBLE( N ) )
  118:       GROWTO = TENTH / ROOTN
  119:       NRMSML = MAX( ONE, EPS3*ROOTN )*SMLNUM
  120: *
  121: *     Form B = H - W*I (except that the subdiagonal elements are not
  122: *     stored).
  123: *
  124:       DO 20 J = 1, N
  125:          DO 10 I = 1, J - 1
  126:             B( I, J ) = H( I, J )
  127:    10    CONTINUE
  128:          B( J, J ) = H( J, J ) - W
  129:    20 CONTINUE
  130: *
  131:       IF( NOINIT ) THEN
  132: *
  133: *        Initialize V.
  134: *
  135:          DO 30 I = 1, N
  136:             V( I ) = EPS3
  137:    30    CONTINUE
  138:       ELSE
  139: *
  140: *        Scale supplied initial vector.
  141: *
  142:          VNORM = DZNRM2( N, V, 1 )
  143:          CALL ZDSCAL( N, ( EPS3*ROOTN ) / MAX( VNORM, NRMSML ), V, 1 )
  144:       END IF
  145: *
  146:       IF( RIGHTV ) THEN
  147: *
  148: *        LU decomposition with partial pivoting of B, replacing zero
  149: *        pivots by EPS3.
  150: *
  151:          DO 60 I = 1, N - 1
  152:             EI = H( I+1, I )
  153:             IF( CABS1( B( I, I ) ).LT.CABS1( EI ) ) THEN
  154: *
  155: *              Interchange rows and eliminate.
  156: *
  157:                X = ZLADIV( B( I, I ), EI )
  158:                B( I, I ) = EI
  159:                DO 40 J = I + 1, N
  160:                   TEMP = B( I+1, J )
  161:                   B( I+1, J ) = B( I, J ) - X*TEMP
  162:                   B( I, J ) = TEMP
  163:    40          CONTINUE
  164:             ELSE
  165: *
  166: *              Eliminate without interchange.
  167: *
  168:                IF( B( I, I ).EQ.ZERO )
  169:      $            B( I, I ) = EPS3
  170:                X = ZLADIV( EI, B( I, I ) )
  171:                IF( X.NE.ZERO ) THEN
  172:                   DO 50 J = I + 1, N
  173:                      B( I+1, J ) = B( I+1, J ) - X*B( I, J )
  174:    50             CONTINUE
  175:                END IF
  176:             END IF
  177:    60    CONTINUE
  178:          IF( B( N, N ).EQ.ZERO )
  179:      $      B( N, N ) = EPS3
  180: *
  181:          TRANS = 'N'
  182: *
  183:       ELSE
  184: *
  185: *        UL decomposition with partial pivoting of B, replacing zero
  186: *        pivots by EPS3.
  187: *
  188:          DO 90 J = N, 2, -1
  189:             EJ = H( J, J-1 )
  190:             IF( CABS1( B( J, J ) ).LT.CABS1( EJ ) ) THEN
  191: *
  192: *              Interchange columns and eliminate.
  193: *
  194:                X = ZLADIV( B( J, J ), EJ )
  195:                B( J, J ) = EJ
  196:                DO 70 I = 1, J - 1
  197:                   TEMP = B( I, J-1 )
  198:                   B( I, J-1 ) = B( I, J ) - X*TEMP
  199:                   B( I, J ) = TEMP
  200:    70          CONTINUE
  201:             ELSE
  202: *
  203: *              Eliminate without interchange.
  204: *
  205:                IF( B( J, J ).EQ.ZERO )
  206:      $            B( J, J ) = EPS3
  207:                X = ZLADIV( EJ, B( J, J ) )
  208:                IF( X.NE.ZERO ) THEN
  209:                   DO 80 I = 1, J - 1
  210:                      B( I, J-1 ) = B( I, J-1 ) - X*B( I, J )
  211:    80             CONTINUE
  212:                END IF
  213:             END IF
  214:    90    CONTINUE
  215:          IF( B( 1, 1 ).EQ.ZERO )
  216:      $      B( 1, 1 ) = EPS3
  217: *
  218:          TRANS = 'C'
  219: *
  220:       END IF
  221: *
  222:       NORMIN = 'N'
  223:       DO 110 ITS = 1, N
  224: *
  225: *        Solve U*x = scale*v for a right eigenvector
  226: *          or U'*x = scale*v for a left eigenvector,
  227: *        overwriting x on v.
  228: *
  229:          CALL ZLATRS( 'Upper', TRANS, 'Nonunit', NORMIN, N, B, LDB, V,
  230:      $                SCALE, RWORK, IERR )
  231:          NORMIN = 'Y'
  232: *
  233: *        Test for sufficient growth in the norm of v.
  234: *
  235:          VNORM = DZASUM( N, V, 1 )
  236:          IF( VNORM.GE.GROWTO*SCALE )
  237:      $      GO TO 120
  238: *
  239: *        Choose new orthogonal starting vector and try again.
  240: *
  241:          RTEMP = EPS3 / ( ROOTN+ONE )
  242:          V( 1 ) = EPS3
  243:          DO 100 I = 2, N
  244:             V( I ) = RTEMP
  245:   100    CONTINUE
  246:          V( N-ITS+1 ) = V( N-ITS+1 ) - EPS3*ROOTN
  247:   110 CONTINUE
  248: *
  249: *     Failure to find eigenvector in N iterations.
  250: *
  251:       INFO = 1
  252: *
  253:   120 CONTINUE
  254: *
  255: *     Normalize eigenvector.
  256: *
  257:       I = IZAMAX( N, V, 1 )
  258:       CALL ZDSCAL( N, ONE / CABS1( V( I ) ), V, 1 )
  259: *
  260:       RETURN
  261: *
  262: *     End of ZLAEIN
  263: *
  264:       END