Return to dlahrd.f CVS log

Up to [local] / rpl / lapack / lapack

Commit initial.

    1:       SUBROUTINE DLAHRD( N, K, NB, A, LDA, TAU, T, LDT, Y, LDY )
    2: *
    3: *  -- LAPACK auxiliary 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            K, LDA, LDT, LDY, N, NB
   10: *     ..
   11: *     .. Array Arguments ..
   12:       DOUBLE PRECISION   A( LDA, * ), T( LDT, NB ), TAU( NB ),
   13:      $                   Y( LDY, NB )
   14: *     ..
   15: *
   16: *  Purpose
   17: *  =======
   18: *
   19: *  DLAHRD reduces the first NB columns of a real general n-by-(n-k+1)
   20: *  matrix A so that elements below the k-th subdiagonal are zero. The
   21: *  reduction is performed by an orthogonal similarity transformation
   22: *  Q' * A * Q. The routine returns the matrices V and T which determine
   23: *  Q as a block reflector I - V*T*V', and also the matrix Y = A * V * T.
   24: *
   25: *  This is an OBSOLETE auxiliary routine. 
   26: *  This routine will be 'deprecated' in a  future release.
   27: *  Please use the new routine DLAHR2 instead.
   28: *
   29: *  Arguments
   30: *  =========
   31: *
   32: *  N       (input) INTEGER
   33: *          The order of the matrix A.
   34: *
   35: *  K       (input) INTEGER
   36: *          The offset for the reduction. Elements below the k-th
   37: *          subdiagonal in the first NB columns are reduced to zero.
   38: *
   39: *  NB      (input) INTEGER
   40: *          The number of columns to be reduced.
   41: *
   42: *  A       (input/output) DOUBLE PRECISION array, dimension (LDA,N-K+1)
   43: *          On entry, the n-by-(n-k+1) general matrix A.
   44: *          On exit, the elements on and above the k-th subdiagonal in
   45: *          the first NB columns are overwritten with the corresponding
   46: *          elements of the reduced matrix; the elements below the k-th
   47: *          subdiagonal, with the array TAU, represent the matrix Q as a
   48: *          product of elementary reflectors. The other columns of A are
   49: *          unchanged. See Further Details.
   50: *
   51: *  LDA     (input) INTEGER
   52: *          The leading dimension of the array A.  LDA >= max(1,N).
   53: *
   54: *  TAU     (output) DOUBLE PRECISION array, dimension (NB)
   55: *          The scalar factors of the elementary reflectors. See Further
   56: *          Details.
   57: *
   58: *  T       (output) DOUBLE PRECISION array, dimension (LDT,NB)
   59: *          The upper triangular matrix T.
   60: *
   61: *  LDT     (input) INTEGER
   62: *          The leading dimension of the array T.  LDT >= NB.
   63: *
   64: *  Y       (output) DOUBLE PRECISION array, dimension (LDY,NB)
   65: *          The n-by-nb matrix Y.
   66: *
   67: *  LDY     (input) INTEGER
   68: *          The leading dimension of the array Y. LDY >= N.
   69: *
   70: *  Further Details
   71: *  ===============
   72: *
   73: *  The matrix Q is represented as a product of nb elementary reflectors
   74: *
   75: *     Q = H(1) H(2) . . . H(nb).
   76: *
   77: *  Each H(i) has the form
   78: *
   79: *     H(i) = I - tau * v * v'
   80: *
   81: *  where tau is a real scalar, and v is a real vector with
   82: *  v(1:i+k-1) = 0, v(i+k) = 1; v(i+k+1:n) is stored on exit in
   83: *  A(i+k+1:n,i), and tau in TAU(i).
   84: *
   85: *  The elements of the vectors v together form the (n-k+1)-by-nb matrix
   86: *  V which is needed, with T and Y, to apply the transformation to the
   87: *  unreduced part of the matrix, using an update of the form:
   88: *  A := (I - V*T*V') * (A - Y*V').
   89: *
   90: *  The contents of A on exit are illustrated by the following example
   91: *  with n = 7, k = 3 and nb = 2:
   92: *
   93: *     ( a   h   a   a   a )
   94: *     ( a   h   a   a   a )
   95: *     ( a   h   a   a   a )
   96: *     ( h   h   a   a   a )
   97: *     ( v1  h   a   a   a )
   98: *     ( v1  v2  a   a   a )
   99: *     ( v1  v2  a   a   a )
  100: *
  101: *  where a denotes an element of the original matrix A, h denotes a
  102: *  modified element of the upper Hessenberg matrix H, and vi denotes an
  103: *  element of the vector defining H(i).
  104: *
  105: *  =====================================================================
  106: *
  107: *     .. Parameters ..
  108:       DOUBLE PRECISION   ZERO, ONE
  109:       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
  110: *     ..
  111: *     .. Local Scalars ..
  112:       INTEGER            I
  113:       DOUBLE PRECISION   EI
  114: *     ..
  115: *     .. External Subroutines ..
  116:       EXTERNAL           DAXPY, DCOPY, DGEMV, DLARFG, DSCAL, DTRMV
  117: *     ..
  118: *     .. Intrinsic Functions ..
  119:       INTRINSIC          MIN
  120: *     ..
  121: *     .. Executable Statements ..
  122: *
  123: *     Quick return if possible
  124: *
  125:       IF( N.LE.1 )
  126:      $   RETURN
  127: *
  128:       DO 10 I = 1, NB
  129:          IF( I.GT.1 ) THEN
  130: *
  131: *           Update A(1:n,i)
  132: *
  133: *           Compute i-th column of A - Y * V'
  134: *
  135:             CALL DGEMV( 'No transpose', N, I-1, -ONE, Y, LDY,
  136:      $                  A( K+I-1, 1 ), LDA, ONE, A( 1, I ), 1 )
  137: *
  138: *           Apply I - V * T' * V' to this column (call it b) from the
  139: *           left, using the last column of T as workspace
  140: *
  141: *           Let  V = ( V1 )   and   b = ( b1 )   (first I-1 rows)
  142: *                    ( V2 )             ( b2 )
  143: *
  144: *           where V1 is unit lower triangular
  145: *
  146: *           w := V1' * b1
  147: *
  148:             CALL DCOPY( I-1, A( K+1, I ), 1, T( 1, NB ), 1 )
  149:             CALL DTRMV( 'Lower', 'Transpose', 'Unit', I-1, A( K+1, 1 ),
  150:      $                  LDA, T( 1, NB ), 1 )
  151: *
  152: *           w := w + V2'*b2
  153: *
  154:             CALL DGEMV( 'Transpose', N-K-I+1, I-1, ONE, A( K+I, 1 ),
  155:      $                  LDA, A( K+I, I ), 1, ONE, T( 1, NB ), 1 )
  156: *
  157: *           w := T'*w
  158: *
  159:             CALL DTRMV( 'Upper', 'Transpose', 'Non-unit', I-1, T, LDT,
  160:      $                  T( 1, NB ), 1 )
  161: *
  162: *           b2 := b2 - V2*w
  163: *
  164:             CALL DGEMV( 'No transpose', N-K-I+1, I-1, -ONE, A( K+I, 1 ),
  165:      $                  LDA, T( 1, NB ), 1, ONE, A( K+I, I ), 1 )
  166: *
  167: *           b1 := b1 - V1*w
  168: *
  169:             CALL DTRMV( 'Lower', 'No transpose', 'Unit', I-1,
  170:      $                  A( K+1, 1 ), LDA, T( 1, NB ), 1 )
  171:             CALL DAXPY( I-1, -ONE, T( 1, NB ), 1, A( K+1, I ), 1 )
  172: *
  173:             A( K+I-1, I-1 ) = EI
  174:          END IF
  175: *
  176: *        Generate the elementary reflector H(i) to annihilate
  177: *        A(k+i+1:n,i)
  178: *
  179:          CALL DLARFG( N-K-I+1, A( K+I, I ), A( MIN( K+I+1, N ), I ), 1,
  180:      $                TAU( I ) )
  181:          EI = A( K+I, I )
  182:          A( K+I, I ) = ONE
  183: *
  184: *        Compute  Y(1:n,i)
  185: *
  186:          CALL DGEMV( 'No transpose', N, N-K-I+1, ONE, A( 1, I+1 ), LDA,
  187:      $               A( K+I, I ), 1, ZERO, Y( 1, I ), 1 )
  188:          CALL DGEMV( 'Transpose', N-K-I+1, I-1, ONE, A( K+I, 1 ), LDA,
  189:      $               A( K+I, I ), 1, ZERO, T( 1, I ), 1 )
  190:          CALL DGEMV( 'No transpose', N, I-1, -ONE, Y, LDY, T( 1, I ), 1,
  191:      $               ONE, Y( 1, I ), 1 )
  192:          CALL DSCAL( N, TAU( I ), Y( 1, I ), 1 )
  193: *
  194: *        Compute T(1:i,i)
  195: *
  196:          CALL DSCAL( I-1, -TAU( I ), T( 1, I ), 1 )
  197:          CALL DTRMV( 'Upper', 'No transpose', 'Non-unit', I-1, T, LDT,
  198:      $               T( 1, I ), 1 )
  199:          T( I, I ) = TAU( I )
  200: *
  201:    10 CONTINUE
  202:       A( K+NB, NB ) = EI
  203: *
  204:       RETURN
  205: *
  206: *     End of DLAHRD
  207: *
  208:       END