Return to zgehd2.f CVS log

Up to [local] / rpl / lapack / lapack

Initial revision

    1:       SUBROUTINE ZGEHD2( N, ILO, IHI, A, LDA, TAU, WORK, INFO )
    2: *
    3: *  -- LAPACK routine (version 3.2) --
    4: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
    5: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
    6: *     November 2006
    7: *
    8: *     .. Scalar Arguments ..
    9:       INTEGER            IHI, ILO, INFO, LDA, N
   10: *     ..
   11: *     .. Array Arguments ..
   12:       COMPLEX*16         A( LDA, * ), TAU( * ), WORK( * )
   13: *     ..
   14: *
   15: *  Purpose
   16: *  =======
   17: *
   18: *  ZGEHD2 reduces a complex general matrix A to upper Hessenberg form H
   19: *  by a unitary similarity transformation:  Q' * A * Q = H .
   20: *
   21: *  Arguments
   22: *  =========
   23: *
   24: *  N       (input) INTEGER
   25: *          The order of the matrix A.  N >= 0.
   26: *
   27: *  ILO     (input) INTEGER
   28: *  IHI     (input) INTEGER
   29: *          It is assumed that A is already upper triangular in rows
   30: *          and columns 1:ILO-1 and IHI+1:N. ILO and IHI are normally
   31: *          set by a previous call to ZGEBAL; otherwise they should be
   32: *          set to 1 and N respectively. See Further Details.
   33: *          1 <= ILO <= IHI <= max(1,N).
   34: *
   35: *  A       (input/output) COMPLEX*16 array, dimension (LDA,N)
   36: *          On entry, the n by n general matrix to be reduced.
   37: *          On exit, the upper triangle and the first subdiagonal of A
   38: *          are overwritten with the upper Hessenberg matrix H, and the
   39: *          elements below the first subdiagonal, with the array TAU,
   40: *          represent the unitary matrix Q as a product of elementary
   41: *          reflectors. See Further Details.
   42: *
   43: *  LDA     (input) INTEGER
   44: *          The leading dimension of the array A.  LDA >= max(1,N).
   45: *
   46: *  TAU     (output) COMPLEX*16 array, dimension (N-1)
   47: *          The scalar factors of the elementary reflectors (see Further
   48: *          Details).
   49: *
   50: *  WORK    (workspace) COMPLEX*16 array, dimension (N)
   51: *
   52: *  INFO    (output) INTEGER
   53: *          = 0:  successful exit
   54: *          < 0:  if INFO = -i, the i-th argument had an illegal value.
   55: *
   56: *  Further Details
   57: *  ===============
   58: *
   59: *  The matrix Q is represented as a product of (ihi-ilo) elementary
   60: *  reflectors
   61: *
   62: *     Q = H(ilo) H(ilo+1) . . . H(ihi-1).
   63: *
   64: *  Each H(i) has the form
   65: *
   66: *     H(i) = I - tau * v * v'
   67: *
   68: *  where tau is a complex scalar, and v is a complex vector with
   69: *  v(1:i) = 0, v(i+1) = 1 and v(ihi+1:n) = 0; v(i+2:ihi) is stored on
   70: *  exit in A(i+2:ihi,i), and tau in TAU(i).
   71: *
   72: *  The contents of A are illustrated by the following example, with
   73: *  n = 7, ilo = 2 and ihi = 6:
   74: *
   75: *  on entry,                        on exit,
   76: *
   77: *  ( a   a   a   a   a   a   a )    (  a   a   h   h   h   h   a )
   78: *  (     a   a   a   a   a   a )    (      a   h   h   h   h   a )
   79: *  (     a   a   a   a   a   a )    (      h   h   h   h   h   h )
   80: *  (     a   a   a   a   a   a )    (      v2  h   h   h   h   h )
   81: *  (     a   a   a   a   a   a )    (      v2  v3  h   h   h   h )
   82: *  (     a   a   a   a   a   a )    (      v2  v3  v4  h   h   h )
   83: *  (                         a )    (                          a )
   84: *
   85: *  where a denotes an element of the original matrix A, h denotes a
   86: *  modified element of the upper Hessenberg matrix H, and vi denotes an
   87: *  element of the vector defining H(i).
   88: *
   89: *  =====================================================================
   90: *
   91: *     .. Parameters ..
   92:       COMPLEX*16         ONE
   93:       PARAMETER          ( ONE = ( 1.0D+0, 0.0D+0 ) )
   94: *     ..
   95: *     .. Local Scalars ..
   96:       INTEGER            I
   97:       COMPLEX*16         ALPHA
   98: *     ..
   99: *     .. External Subroutines ..
  100:       EXTERNAL           XERBLA, ZLARF, ZLARFG
  101: *     ..
  102: *     .. Intrinsic Functions ..
  103:       INTRINSIC          DCONJG, MAX, MIN
  104: *     ..
  105: *     .. Executable Statements ..
  106: *
  107: *     Test the input parameters
  108: *
  109:       INFO = 0
  110:       IF( N.LT.0 ) THEN
  111:          INFO = -1
  112:       ELSE IF( ILO.LT.1 .OR. ILO.GT.MAX( 1, N ) ) THEN
  113:          INFO = -2
  114:       ELSE IF( IHI.LT.MIN( ILO, N ) .OR. IHI.GT.N ) THEN
  115:          INFO = -3
  116:       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
  117:          INFO = -5
  118:       END IF
  119:       IF( INFO.NE.0 ) THEN
  120:          CALL XERBLA( 'ZGEHD2', -INFO )
  121:          RETURN
  122:       END IF
  123: *
  124:       DO 10 I = ILO, IHI - 1
  125: *
  126: *        Compute elementary reflector H(i) to annihilate A(i+2:ihi,i)
  127: *
  128:          ALPHA = A( I+1, I )
  129:          CALL ZLARFG( IHI-I, ALPHA, A( MIN( I+2, N ), I ), 1, TAU( I ) )
  130:          A( I+1, I ) = ONE
  131: *
  132: *        Apply H(i) to A(1:ihi,i+1:ihi) from the right
  133: *
  134:          CALL ZLARF( 'Right', IHI, IHI-I, A( I+1, I ), 1, TAU( I ),
  135:      $               A( 1, I+1 ), LDA, WORK )
  136: *
  137: *        Apply H(i)' to A(i+1:ihi,i+1:n) from the left
  138: *
  139:          CALL ZLARF( 'Left', IHI-I, N-I, A( I+1, I ), 1,
  140:      $               DCONJG( TAU( I ) ), A( I+1, I+1 ), LDA, WORK )
  141: *
  142:          A( I+1, I ) = ALPHA
  143:    10 CONTINUE
  144: *
  145:       RETURN
  146: *
  147: *     End of ZGEHD2
  148: *
  149:       END