Return to dlahr2.f CVS log

Up to [local] / rpl / lapack / lapack

Cohérence

    1:       SUBROUTINE DLAHR2( N, K, NB, A, LDA, TAU, T, LDT, Y, LDY )
    2: *
    3: *  -- LAPACK auxiliary routine (version 3.2.1)                        --
    4: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
    5: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
    6: *  -- April 2009                                                      --
    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: *  DLAHR2 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 auxiliary routine called by DGEHRD.
   26: *
   27: *  Arguments
   28: *  =========
   29: *
   30: *  N       (input) INTEGER
   31: *          The order of the matrix A.
   32: *
   33: *  K       (input) INTEGER
   34: *          The offset for the reduction. Elements below the k-th
   35: *          subdiagonal in the first NB columns are reduced to zero.
   36: *          K < N.
   37: *
   38: *  NB      (input) INTEGER
   39: *          The number of columns to be reduced.
   40: *
   41: *  A       (input/output) DOUBLE PRECISION array, dimension (LDA,N-K+1)
   42: *          On entry, the n-by-(n-k+1) general matrix A.
   43: *          On exit, the elements on and above the k-th subdiagonal in
   44: *          the first NB columns are overwritten with the corresponding
   45: *          elements of the reduced matrix; the elements below the k-th
   46: *          subdiagonal, with the array TAU, represent the matrix Q as a
   47: *          product of elementary reflectors. The other columns of A are
   48: *          unchanged. See Further Details.
   49: *
   50: *  LDA     (input) INTEGER
   51: *          The leading dimension of the array A.  LDA >= max(1,N).
   52: *
   53: *  TAU     (output) DOUBLE PRECISION array, dimension (NB)
   54: *          The scalar factors of the elementary reflectors. See Further
   55: *          Details.
   56: *
   57: *  T       (output) DOUBLE PRECISION array, dimension (LDT,NB)
   58: *          The upper triangular matrix T.
   59: *
   60: *  LDT     (input) INTEGER
   61: *          The leading dimension of the array T.  LDT >= NB.
   62: *
   63: *  Y       (output) DOUBLE PRECISION array, dimension (LDY,NB)
   64: *          The n-by-nb matrix Y.
   65: *
   66: *  LDY     (input) INTEGER
   67: *          The leading dimension of the array Y. LDY >= N.
   68: *
   69: *  Further Details
   70: *  ===============
   71: *
   72: *  The matrix Q is represented as a product of nb elementary reflectors
   73: *
   74: *     Q = H(1) H(2) . . . H(nb).
   75: *
   76: *  Each H(i) has the form
   77: *
   78: *     H(i) = I - tau * v * v'
   79: *
   80: *  where tau is a real scalar, and v is a real vector with
   81: *  v(1:i+k-1) = 0, v(i+k) = 1; v(i+k+1:n) is stored on exit in
   82: *  A(i+k+1:n,i), and tau in TAU(i).
   83: *
   84: *  The elements of the vectors v together form the (n-k+1)-by-nb matrix
   85: *  V which is needed, with T and Y, to apply the transformation to the
   86: *  unreduced part of the matrix, using an update of the form:
   87: *  A := (I - V*T*V') * (A - Y*V').
   88: *
   89: *  The contents of A on exit are illustrated by the following example
   90: *  with n = 7, k = 3 and nb = 2:
   91: *
   92: *     ( a   a   a   a   a )
   93: *     ( a   a   a   a   a )
   94: *     ( a   a   a   a   a )
   95: *     ( h   h   a   a   a )
   96: *     ( v1  h   a   a   a )
   97: *     ( v1  v2  a   a   a )
   98: *     ( v1  v2  a   a   a )
   99: *
  100: *  where a denotes an element of the original matrix A, h denotes a
  101: *  modified element of the upper Hessenberg matrix H, and vi denotes an
  102: *  element of the vector defining H(i).
  103: *
  104: *  This subroutine is a slight modification of LAPACK-3.0's DLAHRD
  105: *  incorporating improvements proposed by Quintana-Orti and Van de
  106: *  Gejin. Note that the entries of A(1:K,2:NB) differ from those
  107: *  returned by the original LAPACK-3.0's DLAHRD routine. (This
  108: *  subroutine is not backward compatible with LAPACK-3.0's DLAHRD.)
  109: *
  110: *  References
  111: *  ==========
  112: *
  113: *  Gregorio Quintana-Orti and Robert van de Geijn, "Improving the
  114: *  performance of reduction to Hessenberg form," ACM Transactions on
  115: *  Mathematical Software, 32(2):180-194, June 2006.
  116: *
  117: *  =====================================================================
  118: *
  119: *     .. Parameters ..
  120:       DOUBLE PRECISION  ZERO, ONE
  121:       PARAMETER          ( ZERO = 0.0D+0, 
  122:      $                     ONE = 1.0D+0 )
  123: *     ..
  124: *     .. Local Scalars ..
  125:       INTEGER            I
  126:       DOUBLE PRECISION  EI
  127: *     ..
  128: *     .. External Subroutines ..
  129:       EXTERNAL           DAXPY, DCOPY, DGEMM, DGEMV, DLACPY,
  130:      $                   DLARFG, DSCAL, DTRMM, DTRMV
  131: *     ..
  132: *     .. Intrinsic Functions ..
  133:       INTRINSIC          MIN
  134: *     ..
  135: *     .. Executable Statements ..
  136: *
  137: *     Quick return if possible
  138: *
  139:       IF( N.LE.1 )
  140:      $   RETURN
  141: *
  142:       DO 10 I = 1, NB
  143:          IF( I.GT.1 ) THEN
  144: *
  145: *           Update A(K+1:N,I)
  146: *
  147: *           Update I-th column of A - Y * V'
  148: *
  149:             CALL DGEMV( 'NO TRANSPOSE', N-K, I-1, -ONE, Y(K+1,1), LDY,
  150:      $                  A( K+I-1, 1 ), LDA, ONE, A( K+1, I ), 1 )
  151: *
  152: *           Apply I - V * T' * V' to this column (call it b) from the
  153: *           left, using the last column of T as workspace
  154: *
  155: *           Let  V = ( V1 )   and   b = ( b1 )   (first I-1 rows)
  156: *                    ( V2 )             ( b2 )
  157: *
  158: *           where V1 is unit lower triangular
  159: *
  160: *           w := V1' * b1
  161: *
  162:             CALL DCOPY( I-1, A( K+1, I ), 1, T( 1, NB ), 1 )
  163:             CALL DTRMV( 'Lower', 'Transpose', 'UNIT', 
  164:      $                  I-1, A( K+1, 1 ),
  165:      $                  LDA, T( 1, NB ), 1 )
  166: *
  167: *           w := w + V2'*b2
  168: *
  169:             CALL DGEMV( 'Transpose', N-K-I+1, I-1, 
  170:      $                  ONE, A( K+I, 1 ),
  171:      $                  LDA, A( K+I, I ), 1, ONE, T( 1, NB ), 1 )
  172: *
  173: *           w := T'*w
  174: *
  175:             CALL DTRMV( 'Upper', 'Transpose', 'NON-UNIT', 
  176:      $                  I-1, T, LDT,
  177:      $                  T( 1, NB ), 1 )
  178: *
  179: *           b2 := b2 - V2*w
  180: *
  181:             CALL DGEMV( 'NO TRANSPOSE', N-K-I+1, I-1, -ONE, 
  182:      $                  A( K+I, 1 ),
  183:      $                  LDA, T( 1, NB ), 1, ONE, A( K+I, I ), 1 )
  184: *
  185: *           b1 := b1 - V1*w
  186: *
  187:             CALL DTRMV( 'Lower', 'NO TRANSPOSE', 
  188:      $                  'UNIT', I-1,
  189:      $                  A( K+1, 1 ), LDA, T( 1, NB ), 1 )
  190:             CALL DAXPY( I-1, -ONE, T( 1, NB ), 1, A( K+1, I ), 1 )
  191: *
  192:             A( K+I-1, I-1 ) = EI
  193:          END IF
  194: *
  195: *        Generate the elementary reflector H(I) to annihilate
  196: *        A(K+I+1:N,I)
  197: *
  198:          CALL DLARFG( N-K-I+1, A( K+I, I ), A( MIN( K+I+1, N ), I ), 1,
  199:      $                TAU( I ) )
  200:          EI = A( K+I, I )
  201:          A( K+I, I ) = ONE
  202: *
  203: *        Compute  Y(K+1:N,I)
  204: *
  205:          CALL DGEMV( 'NO TRANSPOSE', N-K, N-K-I+1, 
  206:      $               ONE, A( K+1, I+1 ),
  207:      $               LDA, A( K+I, I ), 1, ZERO, Y( K+1, I ), 1 )
  208:          CALL DGEMV( 'Transpose', N-K-I+1, I-1, 
  209:      $               ONE, A( K+I, 1 ), LDA,
  210:      $               A( K+I, I ), 1, ZERO, T( 1, I ), 1 )
  211:          CALL DGEMV( 'NO TRANSPOSE', N-K, I-1, -ONE, 
  212:      $               Y( K+1, 1 ), LDY,
  213:      $               T( 1, I ), 1, ONE, Y( K+1, I ), 1 )
  214:          CALL DSCAL( N-K, TAU( I ), Y( K+1, I ), 1 )
  215: *
  216: *        Compute T(1:I,I)
  217: *
  218:          CALL DSCAL( I-1, -TAU( I ), T( 1, I ), 1 )
  219:          CALL DTRMV( 'Upper', 'No Transpose', 'NON-UNIT', 
  220:      $               I-1, T, LDT,
  221:      $               T( 1, I ), 1 )
  222:          T( I, I ) = TAU( I )
  223: *
  224:    10 CONTINUE
  225:       A( K+NB, NB ) = EI
  226: *
  227: *     Compute Y(1:K,1:NB)
  228: *
  229:       CALL DLACPY( 'ALL', K, NB, A( 1, 2 ), LDA, Y, LDY )
  230:       CALL DTRMM( 'RIGHT', 'Lower', 'NO TRANSPOSE', 
  231:      $            'UNIT', K, NB,
  232:      $            ONE, A( K+1, 1 ), LDA, Y, LDY )
  233:       IF( N.GT.K+NB )
  234:      $   CALL DGEMM( 'NO TRANSPOSE', 'NO TRANSPOSE', K, 
  235:      $               NB, N-K-NB, ONE,
  236:      $               A( 1, 2+NB ), LDA, A( K+1+NB, 1 ), LDA, ONE, Y,
  237:      $               LDY )
  238:       CALL DTRMM( 'RIGHT', 'Upper', 'NO TRANSPOSE', 
  239:      $            'NON-UNIT', K, NB,
  240:      $            ONE, T, LDT, Y, LDY )
  241: *
  242:       RETURN
  243: *
  244: *     End of DLAHR2
  245: *
  246:       END