File:  [local] / rpl / lapack / lapack / dlarzb.f
Revision 1.15: download - view: text, annotated - select for diffs - revision graph
Sat Aug 27 15:34:31 2016 UTC (7 years, 8 months ago) by bertrand
Branches: MAIN
CVS tags: rpl-4_1_25, HEAD
Cohérence Lapack.

    1: *> \brief \b DLARZB applies a block reflector or its transpose to a general matrix.
    2: *
    3: *  =========== DOCUMENTATION ===========
    4: *
    5: * Online html documentation available at 
    6: *            http://www.netlib.org/lapack/explore-html/ 
    7: *
    8: *> \htmlonly
    9: *> Download DLARZB + dependencies 
   10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlarzb.f"> 
   11: *> [TGZ]</a> 
   12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlarzb.f"> 
   13: *> [ZIP]</a> 
   14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlarzb.f"> 
   15: *> [TXT]</a>
   16: *> \endhtmlonly 
   17: *
   18: *  Definition:
   19: *  ===========
   20: *
   21: *       SUBROUTINE DLARZB( SIDE, TRANS, DIRECT, STOREV, M, N, K, L, V,
   22: *                          LDV, T, LDT, C, LDC, WORK, LDWORK )
   23:    24: *       .. Scalar Arguments ..
   25: *       CHARACTER          DIRECT, SIDE, STOREV, TRANS
   26: *       INTEGER            K, L, LDC, LDT, LDV, LDWORK, M, N
   27: *       ..
   28: *       .. Array Arguments ..
   29: *       DOUBLE PRECISION   C( LDC, * ), T( LDT, * ), V( LDV, * ),
   30: *      $                   WORK( LDWORK, * )
   31: *       ..
   32: *  
   33: *
   34: *> \par Purpose:
   35: *  =============
   36: *>
   37: *> \verbatim
   38: *>
   39: *> DLARZB applies a real block reflector H or its transpose H**T to
   40: *> a real distributed M-by-N  C from the left or the right.
   41: *>
   42: *> Currently, only STOREV = 'R' and DIRECT = 'B' are supported.
   43: *> \endverbatim
   44: *
   45: *  Arguments:
   46: *  ==========
   47: *
   48: *> \param[in] SIDE
   49: *> \verbatim
   50: *>          SIDE is CHARACTER*1
   51: *>          = 'L': apply H or H**T from the Left
   52: *>          = 'R': apply H or H**T from the Right
   53: *> \endverbatim
   54: *>
   55: *> \param[in] TRANS
   56: *> \verbatim
   57: *>          TRANS is CHARACTER*1
   58: *>          = 'N': apply H (No transpose)
   59: *>          = 'C': apply H**T (Transpose)
   60: *> \endverbatim
   61: *>
   62: *> \param[in] DIRECT
   63: *> \verbatim
   64: *>          DIRECT is CHARACTER*1
   65: *>          Indicates how H is formed from a product of elementary
   66: *>          reflectors
   67: *>          = 'F': H = H(1) H(2) . . . H(k) (Forward, not supported yet)
   68: *>          = 'B': H = H(k) . . . H(2) H(1) (Backward)
   69: *> \endverbatim
   70: *>
   71: *> \param[in] STOREV
   72: *> \verbatim
   73: *>          STOREV is CHARACTER*1
   74: *>          Indicates how the vectors which define the elementary
   75: *>          reflectors are stored:
   76: *>          = 'C': Columnwise                        (not supported yet)
   77: *>          = 'R': Rowwise
   78: *> \endverbatim
   79: *>
   80: *> \param[in] M
   81: *> \verbatim
   82: *>          M is INTEGER
   83: *>          The number of rows of the matrix C.
   84: *> \endverbatim
   85: *>
   86: *> \param[in] N
   87: *> \verbatim
   88: *>          N is INTEGER
   89: *>          The number of columns of the matrix C.
   90: *> \endverbatim
   91: *>
   92: *> \param[in] K
   93: *> \verbatim
   94: *>          K is INTEGER
   95: *>          The order of the matrix T (= the number of elementary
   96: *>          reflectors whose product defines the block reflector).
   97: *> \endverbatim
   98: *>
   99: *> \param[in] L
  100: *> \verbatim
  101: *>          L is INTEGER
  102: *>          The number of columns of the matrix V containing the
  103: *>          meaningful part of the Householder reflectors.
  104: *>          If SIDE = 'L', M >= L >= 0, if SIDE = 'R', N >= L >= 0.
  105: *> \endverbatim
  106: *>
  107: *> \param[in] V
  108: *> \verbatim
  109: *>          V is DOUBLE PRECISION array, dimension (LDV,NV).
  110: *>          If STOREV = 'C', NV = K; if STOREV = 'R', NV = L.
  111: *> \endverbatim
  112: *>
  113: *> \param[in] LDV
  114: *> \verbatim
  115: *>          LDV is INTEGER
  116: *>          The leading dimension of the array V.
  117: *>          If STOREV = 'C', LDV >= L; if STOREV = 'R', LDV >= K.
  118: *> \endverbatim
  119: *>
  120: *> \param[in] T
  121: *> \verbatim
  122: *>          T is DOUBLE PRECISION array, dimension (LDT,K)
  123: *>          The triangular K-by-K matrix T in the representation of the
  124: *>          block reflector.
  125: *> \endverbatim
  126: *>
  127: *> \param[in] LDT
  128: *> \verbatim
  129: *>          LDT is INTEGER
  130: *>          The leading dimension of the array T. LDT >= K.
  131: *> \endverbatim
  132: *>
  133: *> \param[in,out] C
  134: *> \verbatim
  135: *>          C is DOUBLE PRECISION array, dimension (LDC,N)
  136: *>          On entry, the M-by-N matrix C.
  137: *>          On exit, C is overwritten by H*C or H**T*C or C*H or C*H**T.
  138: *> \endverbatim
  139: *>
  140: *> \param[in] LDC
  141: *> \verbatim
  142: *>          LDC is INTEGER
  143: *>          The leading dimension of the array C. LDC >= max(1,M).
  144: *> \endverbatim
  145: *>
  146: *> \param[out] WORK
  147: *> \verbatim
  148: *>          WORK is DOUBLE PRECISION array, dimension (LDWORK,K)
  149: *> \endverbatim
  150: *>
  151: *> \param[in] LDWORK
  152: *> \verbatim
  153: *>          LDWORK is INTEGER
  154: *>          The leading dimension of the array WORK.
  155: *>          If SIDE = 'L', LDWORK >= max(1,N);
  156: *>          if SIDE = 'R', LDWORK >= max(1,M).
  157: *> \endverbatim
  158: *
  159: *  Authors:
  160: *  ========
  161: *
  162: *> \author Univ. of Tennessee 
  163: *> \author Univ. of California Berkeley 
  164: *> \author Univ. of Colorado Denver 
  165: *> \author NAG Ltd. 
  166: *
  167: *> \date September 2012
  168: *
  169: *> \ingroup doubleOTHERcomputational
  170: *
  171: *> \par Contributors:
  172: *  ==================
  173: *>
  174: *>    A. Petitet, Computer Science Dept., Univ. of Tenn., Knoxville, USA
  175: *
  176: *> \par Further Details:
  177: *  =====================
  178: *>
  179: *> \verbatim
  180: *> \endverbatim
  181: *>
  182: *  =====================================================================
  183:       SUBROUTINE DLARZB( SIDE, TRANS, DIRECT, STOREV, M, N, K, L, V,
  184:      $                   LDV, T, LDT, C, LDC, WORK, LDWORK )
  185: *
  186: *  -- LAPACK computational routine (version 3.4.2) --
  187: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
  188: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  189: *     September 2012
  190: *
  191: *     .. Scalar Arguments ..
  192:       CHARACTER          DIRECT, SIDE, STOREV, TRANS
  193:       INTEGER            K, L, LDC, LDT, LDV, LDWORK, M, N
  194: *     ..
  195: *     .. Array Arguments ..
  196:       DOUBLE PRECISION   C( LDC, * ), T( LDT, * ), V( LDV, * ),
  197:      $                   WORK( LDWORK, * )
  198: *     ..
  199: *
  200: *  =====================================================================
  201: *
  202: *     .. Parameters ..
  203:       DOUBLE PRECISION   ONE
  204:       PARAMETER          ( ONE = 1.0D+0 )
  205: *     ..
  206: *     .. Local Scalars ..
  207:       CHARACTER          TRANST
  208:       INTEGER            I, INFO, J
  209: *     ..
  210: *     .. External Functions ..
  211:       LOGICAL            LSAME
  212:       EXTERNAL           LSAME
  213: *     ..
  214: *     .. External Subroutines ..
  215:       EXTERNAL           DCOPY, DGEMM, DTRMM, XERBLA
  216: *     ..
  217: *     .. Executable Statements ..
  218: *
  219: *     Quick return if possible
  220: *
  221:       IF( M.LE.0 .OR. N.LE.0 )
  222:      $   RETURN
  223: *
  224: *     Check for currently supported options
  225: *
  226:       INFO = 0
  227:       IF( .NOT.LSAME( DIRECT, 'B' ) ) THEN
  228:          INFO = -3
  229:       ELSE IF( .NOT.LSAME( STOREV, 'R' ) ) THEN
  230:          INFO = -4
  231:       END IF
  232:       IF( INFO.NE.0 ) THEN
  233:          CALL XERBLA( 'DLARZB', -INFO )
  234:          RETURN
  235:       END IF
  236: *
  237:       IF( LSAME( TRANS, 'N' ) ) THEN
  238:          TRANST = 'T'
  239:       ELSE
  240:          TRANST = 'N'
  241:       END IF
  242: *
  243:       IF( LSAME( SIDE, 'L' ) ) THEN
  244: *
  245: *        Form  H * C  or  H**T * C
  246: *
  247: *        W( 1:n, 1:k ) = C( 1:k, 1:n )**T
  248: *
  249:          DO 10 J = 1, K
  250:             CALL DCOPY( N, C( J, 1 ), LDC, WORK( 1, J ), 1 )
  251:    10    CONTINUE
  252: *
  253: *        W( 1:n, 1:k ) = W( 1:n, 1:k ) + ...
  254: *                        C( m-l+1:m, 1:n )**T * V( 1:k, 1:l )**T
  255: *
  256:          IF( L.GT.0 )
  257:      $      CALL DGEMM( 'Transpose', 'Transpose', N, K, L, ONE,
  258:      $                  C( M-L+1, 1 ), LDC, V, LDV, ONE, WORK, LDWORK )
  259: *
  260: *        W( 1:n, 1:k ) = W( 1:n, 1:k ) * T**T  or  W( 1:m, 1:k ) * T
  261: *
  262:          CALL DTRMM( 'Right', 'Lower', TRANST, 'Non-unit', N, K, ONE, T,
  263:      $               LDT, WORK, LDWORK )
  264: *
  265: *        C( 1:k, 1:n ) = C( 1:k, 1:n ) - W( 1:n, 1:k )**T
  266: *
  267:          DO 30 J = 1, N
  268:             DO 20 I = 1, K
  269:                C( I, J ) = C( I, J ) - WORK( J, I )
  270:    20       CONTINUE
  271:    30    CONTINUE
  272: *
  273: *        C( m-l+1:m, 1:n ) = C( m-l+1:m, 1:n ) - ...
  274: *                            V( 1:k, 1:l )**T * W( 1:n, 1:k )**T
  275: *
  276:          IF( L.GT.0 )
  277:      $      CALL DGEMM( 'Transpose', 'Transpose', L, N, K, -ONE, V, LDV,
  278:      $                  WORK, LDWORK, ONE, C( M-L+1, 1 ), LDC )
  279: *
  280:       ELSE IF( LSAME( SIDE, 'R' ) ) THEN
  281: *
  282: *        Form  C * H  or  C * H**T
  283: *
  284: *        W( 1:m, 1:k ) = C( 1:m, 1:k )
  285: *
  286:          DO 40 J = 1, K
  287:             CALL DCOPY( M, C( 1, J ), 1, WORK( 1, J ), 1 )
  288:    40    CONTINUE
  289: *
  290: *        W( 1:m, 1:k ) = W( 1:m, 1:k ) + ...
  291: *                        C( 1:m, n-l+1:n ) * V( 1:k, 1:l )**T
  292: *
  293:          IF( L.GT.0 )
  294:      $      CALL DGEMM( 'No transpose', 'Transpose', M, K, L, ONE,
  295:      $                  C( 1, N-L+1 ), LDC, V, LDV, ONE, WORK, LDWORK )
  296: *
  297: *        W( 1:m, 1:k ) = W( 1:m, 1:k ) * T  or  W( 1:m, 1:k ) * T**T
  298: *
  299:          CALL DTRMM( 'Right', 'Lower', TRANS, 'Non-unit', M, K, ONE, T,
  300:      $               LDT, WORK, LDWORK )
  301: *
  302: *        C( 1:m, 1:k ) = C( 1:m, 1:k ) - W( 1:m, 1:k )
  303: *
  304:          DO 60 J = 1, K
  305:             DO 50 I = 1, M
  306:                C( I, J ) = C( I, J ) - WORK( I, J )
  307:    50       CONTINUE
  308:    60    CONTINUE
  309: *
  310: *        C( 1:m, n-l+1:n ) = C( 1:m, n-l+1:n ) - ...
  311: *                            W( 1:m, 1:k ) * V( 1:k, 1:l )
  312: *
  313:          IF( L.GT.0 )
  314:      $      CALL DGEMM( 'No transpose', 'No transpose', M, L, K, -ONE,
  315:      $                  WORK, LDWORK, V, LDV, ONE, C( 1, N-L+1 ), LDC )
  316: *
  317:       END IF
  318: *
  319:       RETURN
  320: *
  321: *     End of DLARZB
  322: *
  323:       END

CVSweb interface <joel.bertrand@systella.fr>