Annotation of rpl/lapack/lapack/dlarzb.f, revision 1.9

1.9     ! bertrand    1: *> \brief \b DLARZB
        !             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 November 2011
        !           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: *  =====================================================================
1.1       bertrand  183:       SUBROUTINE DLARZB( SIDE, TRANS, DIRECT, STOREV, M, N, K, L, V,
                    184:      $                   LDV, T, LDT, C, LDC, WORK, LDWORK )
                    185: *
1.9     ! bertrand  186: *  -- LAPACK computational routine (version 3.4.0) --
1.1       bertrand  187: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
                    188: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.9     ! bertrand  189: *     November 2011
1.1       bertrand  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: *
1.8       bertrand  245: *        Form  H * C  or  H**T * C
1.1       bertrand  246: *
1.8       bertrand  247: *        W( 1:n, 1:k ) = C( 1:k, 1:n )**T
1.1       bertrand  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 ) + ...
1.8       bertrand  254: *                        C( m-l+1:m, 1:n )**T * V( 1:k, 1:l )**T
1.1       bertrand  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: *
1.8       bertrand  260: *        W( 1:n, 1:k ) = W( 1:n, 1:k ) * T**T  or  W( 1:m, 1:k ) * T
1.1       bertrand  261: *
                    262:          CALL DTRMM( 'Right', 'Lower', TRANST, 'Non-unit', N, K, ONE, T,
                    263:      $               LDT, WORK, LDWORK )
                    264: *
1.8       bertrand  265: *        C( 1:k, 1:n ) = C( 1:k, 1:n ) - W( 1:n, 1:k )**T
1.1       bertrand  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 ) - ...
1.8       bertrand  274: *                            V( 1:k, 1:l )**T * W( 1:n, 1:k )**T
1.1       bertrand  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: *
1.8       bertrand  282: *        Form  C * H  or  C * H**T
1.1       bertrand  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 ) + ...
1.8       bertrand  291: *                        C( 1:m, n-l+1:n ) * V( 1:k, 1:l )**T
1.1       bertrand  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: *
1.8       bertrand  297: *        W( 1:m, 1:k ) = W( 1:m, 1:k ) * T  or  W( 1:m, 1:k ) * T**T
1.1       bertrand  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>