Return to dlahrd.f CVS log

Up to [local] / rpl / lapack / lapack

Mise à jour de Lapack.

    1: *> \brief \b DLAHRD reduces the first nb columns of a general rectangular matrix A so that elements below the k-th subdiagonal are zero, and returns auxiliary matrices which are needed to apply the transformation to the unreduced part of A.
    2: *
    3: *  =========== DOCUMENTATION ===========
    4: *
    5: * Online html documentation available at
    6: *            http://www.netlib.org/lapack/explore-html/
    7: *
    8: *> \htmlonly
    9: *> Download DLAHRD + dependencies
   10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlahrd.f">
   11: *> [TGZ]</a>
   12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlahrd.f">
   13: *> [ZIP]</a>
   14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlahrd.f">
   15: *> [TXT]</a>
   16: *> \endhtmlonly
   17: *
   18: *  Definition:
   19: *  ===========
   20: *
   21: *       SUBROUTINE DLAHRD( N, K, NB, A, LDA, TAU, T, LDT, Y, LDY )
   22: *
   23: *       .. Scalar Arguments ..
   24: *       INTEGER            K, LDA, LDT, LDY, N, NB
   25: *       ..
   26: *       .. Array Arguments ..
   27: *       DOUBLE PRECISION   A( LDA, * ), T( LDT, NB ), TAU( NB ),
   28: *      $                   Y( LDY, NB )
   29: *       ..
   30: *
   31: *
   32: *> \par Purpose:
   33: *  =============
   34: *>
   35: *> \verbatim
   36: *>
   37: *> This routine is deprecated and has been replaced by routine DLAHR2.
   38: *>
   39: *> DLAHRD reduces the first NB columns of a real general n-by-(n-k+1)
   40: *> matrix A so that elements below the k-th subdiagonal are zero. The
   41: *> reduction is performed by an orthogonal similarity transformation
   42: *> Q**T * A * Q. The routine returns the matrices V and T which determine
   43: *> Q as a block reflector I - V*T*V**T, and also the matrix Y = A * V * T.
   44: *> \endverbatim
   45: *
   46: *  Arguments:
   47: *  ==========
   48: *
   49: *> \param[in] N
   50: *> \verbatim
   51: *>          N is INTEGER
   52: *>          The order of the matrix A.
   53: *> \endverbatim
   54: *>
   55: *> \param[in] K
   56: *> \verbatim
   57: *>          K is INTEGER
   58: *>          The offset for the reduction. Elements below the k-th
   59: *>          subdiagonal in the first NB columns are reduced to zero.
   60: *> \endverbatim
   61: *>
   62: *> \param[in] NB
   63: *> \verbatim
   64: *>          NB is INTEGER
   65: *>          The number of columns to be reduced.
   66: *> \endverbatim
   67: *>
   68: *> \param[in,out] A
   69: *> \verbatim
   70: *>          A is DOUBLE PRECISION array, dimension (LDA,N-K+1)
   71: *>          On entry, the n-by-(n-k+1) general matrix A.
   72: *>          On exit, the elements on and above the k-th subdiagonal in
   73: *>          the first NB columns are overwritten with the corresponding
   74: *>          elements of the reduced matrix; the elements below the k-th
   75: *>          subdiagonal, with the array TAU, represent the matrix Q as a
   76: *>          product of elementary reflectors. The other columns of A are
   77: *>          unchanged. See Further Details.
   78: *> \endverbatim
   79: *>
   80: *> \param[in] LDA
   81: *> \verbatim
   82: *>          LDA is INTEGER
   83: *>          The leading dimension of the array A.  LDA >= max(1,N).
   84: *> \endverbatim
   85: *>
   86: *> \param[out] TAU
   87: *> \verbatim
   88: *>          TAU is DOUBLE PRECISION array, dimension (NB)
   89: *>          The scalar factors of the elementary reflectors. See Further
   90: *>          Details.
   91: *> \endverbatim
   92: *>
   93: *> \param[out] T
   94: *> \verbatim
   95: *>          T is DOUBLE PRECISION array, dimension (LDT,NB)
   96: *>          The upper triangular matrix T.
   97: *> \endverbatim
   98: *>
   99: *> \param[in] LDT
  100: *> \verbatim
  101: *>          LDT is INTEGER
  102: *>          The leading dimension of the array T.  LDT >= NB.
  103: *> \endverbatim
  104: *>
  105: *> \param[out] Y
  106: *> \verbatim
  107: *>          Y is DOUBLE PRECISION array, dimension (LDY,NB)
  108: *>          The n-by-nb matrix Y.
  109: *> \endverbatim
  110: *>
  111: *> \param[in] LDY
  112: *> \verbatim
  113: *>          LDY is INTEGER
  114: *>          The leading dimension of the array Y. LDY >= N.
  115: *> \endverbatim
  116: *
  117: *  Authors:
  118: *  ========
  119: *
  120: *> \author Univ. of Tennessee
  121: *> \author Univ. of California Berkeley
  122: *> \author Univ. of Colorado Denver
  123: *> \author NAG Ltd.
  124: *
  125: *> \date December 2016
  126: *
  127: *> \ingroup doubleOTHERauxiliary
  128: *
  129: *> \par Further Details:
  130: *  =====================
  131: *>
  132: *> \verbatim
  133: *>
  134: *>  The matrix Q is represented as a product of nb elementary reflectors
  135: *>
  136: *>     Q = H(1) H(2) . . . H(nb).
  137: *>
  138: *>  Each H(i) has the form
  139: *>
  140: *>     H(i) = I - tau * v * v**T
  141: *>
  142: *>  where tau is a real scalar, and v is a real vector with
  143: *>  v(1:i+k-1) = 0, v(i+k) = 1; v(i+k+1:n) is stored on exit in
  144: *>  A(i+k+1:n,i), and tau in TAU(i).
  145: *>
  146: *>  The elements of the vectors v together form the (n-k+1)-by-nb matrix
  147: *>  V which is needed, with T and Y, to apply the transformation to the
  148: *>  unreduced part of the matrix, using an update of the form:
  149: *>  A := (I - V*T*V**T) * (A - Y*V**T).
  150: *>
  151: *>  The contents of A on exit are illustrated by the following example
  152: *>  with n = 7, k = 3 and nb = 2:
  153: *>
  154: *>     ( a   h   a   a   a )
  155: *>     ( a   h   a   a   a )
  156: *>     ( a   h   a   a   a )
  157: *>     ( h   h   a   a   a )
  158: *>     ( v1  h   a   a   a )
  159: *>     ( v1  v2  a   a   a )
  160: *>     ( v1  v2  a   a   a )
  161: *>
  162: *>  where a denotes an element of the original matrix A, h denotes a
  163: *>  modified element of the upper Hessenberg matrix H, and vi denotes an
  164: *>  element of the vector defining H(i).
  165: *> \endverbatim
  166: *>
  167: *  =====================================================================
  168:       SUBROUTINE DLAHRD( N, K, NB, A, LDA, TAU, T, LDT, Y, LDY )
  169: *
  170: *  -- LAPACK auxiliary routine (version 3.7.0) --
  171: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
  172: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  173: *     December 2016
  174: *
  175: *     .. Scalar Arguments ..
  176:       INTEGER            K, LDA, LDT, LDY, N, NB
  177: *     ..
  178: *     .. Array Arguments ..
  179:       DOUBLE PRECISION   A( LDA, * ), T( LDT, NB ), TAU( NB ),
  180:      $                   Y( LDY, NB )
  181: *     ..
  182: *
  183: *  =====================================================================
  184: *
  185: *     .. Parameters ..
  186:       DOUBLE PRECISION   ZERO, ONE
  187:       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
  188: *     ..
  189: *     .. Local Scalars ..
  190:       INTEGER            I
  191:       DOUBLE PRECISION   EI
  192: *     ..
  193: *     .. External Subroutines ..
  194:       EXTERNAL           DAXPY, DCOPY, DGEMV, DLARFG, DSCAL, DTRMV
  195: *     ..
  196: *     .. Intrinsic Functions ..
  197:       INTRINSIC          MIN
  198: *     ..
  199: *     .. Executable Statements ..
  200: *
  201: *     Quick return if possible
  202: *
  203:       IF( N.LE.1 )
  204:      $   RETURN
  205: *
  206:       DO 10 I = 1, NB
  207:          IF( I.GT.1 ) THEN
  208: *
  209: *           Update A(1:n,i)
  210: *
  211: *           Compute i-th column of A - Y * V**T
  212: *
  213:             CALL DGEMV( 'No transpose', N, I-1, -ONE, Y, LDY,
  214:      $                  A( K+I-1, 1 ), LDA, ONE, A( 1, I ), 1 )
  215: *
  216: *           Apply I - V * T**T * V**T to this column (call it b) from the
  217: *           left, using the last column of T as workspace
  218: *
  219: *           Let  V = ( V1 )   and   b = ( b1 )   (first I-1 rows)
  220: *                    ( V2 )             ( b2 )
  221: *
  222: *           where V1 is unit lower triangular
  223: *
  224: *           w := V1**T * b1
  225: *
  226:             CALL DCOPY( I-1, A( K+1, I ), 1, T( 1, NB ), 1 )
  227:             CALL DTRMV( 'Lower', 'Transpose', 'Unit', I-1, A( K+1, 1 ),
  228:      $                  LDA, T( 1, NB ), 1 )
  229: *
  230: *           w := w + V2**T *b2
  231: *
  232:             CALL DGEMV( 'Transpose', N-K-I+1, I-1, ONE, A( K+I, 1 ),
  233:      $                  LDA, A( K+I, I ), 1, ONE, T( 1, NB ), 1 )
  234: *
  235: *           w := T**T *w
  236: *
  237:             CALL DTRMV( 'Upper', 'Transpose', 'Non-unit', I-1, T, LDT,
  238:      $                  T( 1, NB ), 1 )
  239: *
  240: *           b2 := b2 - V2*w
  241: *
  242:             CALL DGEMV( 'No transpose', N-K-I+1, I-1, -ONE, A( K+I, 1 ),
  243:      $                  LDA, T( 1, NB ), 1, ONE, A( K+I, I ), 1 )
  244: *
  245: *           b1 := b1 - V1*w
  246: *
  247:             CALL DTRMV( 'Lower', 'No transpose', 'Unit', I-1,
  248:      $                  A( K+1, 1 ), LDA, T( 1, NB ), 1 )
  249:             CALL DAXPY( I-1, -ONE, T( 1, NB ), 1, A( K+1, I ), 1 )
  250: *
  251:             A( K+I-1, I-1 ) = EI
  252:          END IF
  253: *
  254: *        Generate the elementary reflector H(i) to annihilate
  255: *        A(k+i+1:n,i)
  256: *
  257:          CALL DLARFG( N-K-I+1, A( K+I, I ), A( MIN( K+I+1, N ), I ), 1,
  258:      $                TAU( I ) )
  259:          EI = A( K+I, I )
  260:          A( K+I, I ) = ONE
  261: *
  262: *        Compute  Y(1:n,i)
  263: *
  264:          CALL DGEMV( 'No transpose', N, N-K-I+1, ONE, A( 1, I+1 ), LDA,
  265:      $               A( K+I, I ), 1, ZERO, Y( 1, I ), 1 )
  266:          CALL DGEMV( 'Transpose', N-K-I+1, I-1, ONE, A( K+I, 1 ), LDA,
  267:      $               A( K+I, I ), 1, ZERO, T( 1, I ), 1 )
  268:          CALL DGEMV( 'No transpose', N, I-1, -ONE, Y, LDY, T( 1, I ), 1,
  269:      $               ONE, Y( 1, I ), 1 )
  270:          CALL DSCAL( N, TAU( I ), Y( 1, I ), 1 )
  271: *
  272: *        Compute T(1:i,i)
  273: *
  274:          CALL DSCAL( I-1, -TAU( I ), T( 1, I ), 1 )
  275:          CALL DTRMV( 'Upper', 'No transpose', 'Non-unit', I-1, T, LDT,
  276:      $               T( 1, I ), 1 )
  277:          T( I, I ) = TAU( I )
  278: *
  279:    10 CONTINUE
  280:       A( K+NB, NB ) = EI
  281: *
  282:       RETURN
  283: *
  284: *     End of DLAHRD
  285: *
  286:       END