LAPACK 3.3.1
Linear Algebra PACKage

clarzb.f

Go to the documentation of this file.
00001       SUBROUTINE CLARZB( SIDE, TRANS, DIRECT, STOREV, M, N, K, L, V,
00002      $                   LDV, T, LDT, C, LDC, WORK, LDWORK )
00003 *
00004 *  -- LAPACK routine (version 3.3.1) --
00005 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00006 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00007 *  -- April 2011                                                      --
00008 *
00009 *     .. Scalar Arguments ..
00010       CHARACTER          DIRECT, SIDE, STOREV, TRANS
00011       INTEGER            K, L, LDC, LDT, LDV, LDWORK, M, N
00012 *     ..
00013 *     .. Array Arguments ..
00014       COMPLEX            C( LDC, * ), T( LDT, * ), V( LDV, * ),
00015      $                   WORK( LDWORK, * )
00016 *     ..
00017 *
00018 *  Purpose
00019 *  =======
00020 *
00021 *  CLARZB applies a complex block reflector H or its transpose H**H
00022 *  to a complex distributed M-by-N  C from the left or the right.
00023 *
00024 *  Currently, only STOREV = 'R' and DIRECT = 'B' are supported.
00025 *
00026 *  Arguments
00027 *  =========
00028 *
00029 *  SIDE    (input) CHARACTER*1
00030 *          = 'L': apply H or H**H from the Left
00031 *          = 'R': apply H or H**H from the Right
00032 *
00033 *  TRANS   (input) CHARACTER*1
00034 *          = 'N': apply H (No transpose)
00035 *          = 'C': apply H**H (Conjugate transpose)
00036 *
00037 *  DIRECT  (input) CHARACTER*1
00038 *          Indicates how H is formed from a product of elementary
00039 *          reflectors
00040 *          = 'F': H = H(1) H(2) . . . H(k) (Forward, not supported yet)
00041 *          = 'B': H = H(k) . . . H(2) H(1) (Backward)
00042 *
00043 *  STOREV  (input) CHARACTER*1
00044 *          Indicates how the vectors which define the elementary
00045 *          reflectors are stored:
00046 *          = 'C': Columnwise                        (not supported yet)
00047 *          = 'R': Rowwise
00048 *
00049 *  M       (input) INTEGER
00050 *          The number of rows of the matrix C.
00051 *
00052 *  N       (input) INTEGER
00053 *          The number of columns of the matrix C.
00054 *
00055 *  K       (input) INTEGER
00056 *          The order of the matrix T (= the number of elementary
00057 *          reflectors whose product defines the block reflector).
00058 *
00059 *  L       (input) INTEGER
00060 *          The number of columns of the matrix V containing the
00061 *          meaningful part of the Householder reflectors.
00062 *          If SIDE = 'L', M >= L >= 0, if SIDE = 'R', N >= L >= 0.
00063 *
00064 *  V       (input) COMPLEX array, dimension (LDV,NV).
00065 *          If STOREV = 'C', NV = K; if STOREV = 'R', NV = L.
00066 *
00067 *  LDV     (input) INTEGER
00068 *          The leading dimension of the array V.
00069 *          If STOREV = 'C', LDV >= L; if STOREV = 'R', LDV >= K.
00070 *
00071 *  T       (input) COMPLEX array, dimension (LDT,K)
00072 *          The triangular K-by-K matrix T in the representation of the
00073 *          block reflector.
00074 *
00075 *  LDT     (input) INTEGER
00076 *          The leading dimension of the array T. LDT >= K.
00077 *
00078 *  C       (input/output) COMPLEX array, dimension (LDC,N)
00079 *          On entry, the M-by-N matrix C.
00080 *          On exit, C is overwritten by H*C or H**H*C or C*H or C*H**H.
00081 *
00082 *  LDC     (input) INTEGER
00083 *          The leading dimension of the array C. LDC >= max(1,M).
00084 *
00085 *  WORK    (workspace) COMPLEX array, dimension (LDWORK,K)
00086 *
00087 *  LDWORK  (input) INTEGER
00088 *          The leading dimension of the array WORK.
00089 *          If SIDE = 'L', LDWORK >= max(1,N);
00090 *          if SIDE = 'R', LDWORK >= max(1,M).
00091 *
00092 *  Further Details
00093 *  ===============
00094 *
00095 *  Based on contributions by
00096 *    A. Petitet, Computer Science Dept., Univ. of Tenn., Knoxville, USA
00097 *
00098 *  =====================================================================
00099 *
00100 *     .. Parameters ..
00101       COMPLEX            ONE
00102       PARAMETER          ( ONE = ( 1.0E+0, 0.0E+0 ) )
00103 *     ..
00104 *     .. Local Scalars ..
00105       CHARACTER          TRANST
00106       INTEGER            I, INFO, J
00107 *     ..
00108 *     .. External Functions ..
00109       LOGICAL            LSAME
00110       EXTERNAL           LSAME
00111 *     ..
00112 *     .. External Subroutines ..
00113       EXTERNAL           CCOPY, CGEMM, CLACGV, CTRMM, XERBLA
00114 *     ..
00115 *     .. Executable Statements ..
00116 *
00117 *     Quick return if possible
00118 *
00119       IF( M.LE.0 .OR. N.LE.0 )
00120      $   RETURN
00121 *
00122 *     Check for currently supported options
00123 *
00124       INFO = 0
00125       IF( .NOT.LSAME( DIRECT, 'B' ) ) THEN
00126          INFO = -3
00127       ELSE IF( .NOT.LSAME( STOREV, 'R' ) ) THEN
00128          INFO = -4
00129       END IF
00130       IF( INFO.NE.0 ) THEN
00131          CALL XERBLA( 'CLARZB', -INFO )
00132          RETURN
00133       END IF
00134 *
00135       IF( LSAME( TRANS, 'N' ) ) THEN
00136          TRANST = 'C'
00137       ELSE
00138          TRANST = 'N'
00139       END IF
00140 *
00141       IF( LSAME( SIDE, 'L' ) ) THEN
00142 *
00143 *        Form  H * C  or  H**H * C
00144 *
00145 *        W( 1:n, 1:k ) = C( 1:k, 1:n )**H
00146 *
00147          DO 10 J = 1, K
00148             CALL CCOPY( N, C( J, 1 ), LDC, WORK( 1, J ), 1 )
00149    10    CONTINUE
00150 *
00151 *        W( 1:n, 1:k ) = W( 1:n, 1:k ) + ...
00152 *                        C( m-l+1:m, 1:n )**H * V( 1:k, 1:l )**T
00153 *
00154          IF( L.GT.0 )
00155      $      CALL CGEMM( 'Transpose', 'Conjugate transpose', N, K, L,
00156      $                  ONE, C( M-L+1, 1 ), LDC, V, LDV, ONE, WORK,
00157      $                  LDWORK )
00158 *
00159 *        W( 1:n, 1:k ) = W( 1:n, 1:k ) * T**T  or  W( 1:m, 1:k ) * T
00160 *
00161          CALL CTRMM( 'Right', 'Lower', TRANST, 'Non-unit', N, K, ONE, T,
00162      $               LDT, WORK, LDWORK )
00163 *
00164 *        C( 1:k, 1:n ) = C( 1:k, 1:n ) - W( 1:n, 1:k )**H
00165 *
00166          DO 30 J = 1, N
00167             DO 20 I = 1, K
00168                C( I, J ) = C( I, J ) - WORK( J, I )
00169    20       CONTINUE
00170    30    CONTINUE
00171 *
00172 *        C( m-l+1:m, 1:n ) = C( m-l+1:m, 1:n ) - ...
00173 *                            V( 1:k, 1:l )**H * W( 1:n, 1:k )**H
00174 *
00175          IF( L.GT.0 )
00176      $      CALL CGEMM( 'Transpose', 'Transpose', L, N, K, -ONE, V, LDV,
00177      $                  WORK, LDWORK, ONE, C( M-L+1, 1 ), LDC )
00178 *
00179       ELSE IF( LSAME( SIDE, 'R' ) ) THEN
00180 *
00181 *        Form  C * H  or  C * H**H
00182 *
00183 *        W( 1:m, 1:k ) = C( 1:m, 1:k )
00184 *
00185          DO 40 J = 1, K
00186             CALL CCOPY( M, C( 1, J ), 1, WORK( 1, J ), 1 )
00187    40    CONTINUE
00188 *
00189 *        W( 1:m, 1:k ) = W( 1:m, 1:k ) + ...
00190 *                        C( 1:m, n-l+1:n ) * V( 1:k, 1:l )**H
00191 *
00192          IF( L.GT.0 )
00193      $      CALL CGEMM( 'No transpose', 'Transpose', M, K, L, ONE,
00194      $                  C( 1, N-L+1 ), LDC, V, LDV, ONE, WORK, LDWORK )
00195 *
00196 *        W( 1:m, 1:k ) = W( 1:m, 1:k ) * conjg( T )  or
00197 *                        W( 1:m, 1:k ) * T**H
00198 *
00199          DO 50 J = 1, K
00200             CALL CLACGV( K-J+1, T( J, J ), 1 )
00201    50    CONTINUE
00202          CALL CTRMM( 'Right', 'Lower', TRANS, 'Non-unit', M, K, ONE, T,
00203      $               LDT, WORK, LDWORK )
00204          DO 60 J = 1, K
00205             CALL CLACGV( K-J+1, T( J, J ), 1 )
00206    60    CONTINUE
00207 *
00208 *        C( 1:m, 1:k ) = C( 1:m, 1:k ) - W( 1:m, 1:k )
00209 *
00210          DO 80 J = 1, K
00211             DO 70 I = 1, M
00212                C( I, J ) = C( I, J ) - WORK( I, J )
00213    70       CONTINUE
00214    80    CONTINUE
00215 *
00216 *        C( 1:m, n-l+1:n ) = C( 1:m, n-l+1:n ) - ...
00217 *                            W( 1:m, 1:k ) * conjg( V( 1:k, 1:l ) )
00218 *
00219          DO 90 J = 1, L
00220             CALL CLACGV( K, V( 1, J ), 1 )
00221    90    CONTINUE
00222          IF( L.GT.0 )
00223      $      CALL CGEMM( 'No transpose', 'No transpose', M, L, K, -ONE,
00224      $                  WORK, LDWORK, V, LDV, ONE, C( 1, N-L+1 ), LDC )
00225          DO 100 J = 1, L
00226             CALL CLACGV( K, V( 1, J ), 1 )
00227   100    CONTINUE
00228 *
00229       END IF
00230 *
00231       RETURN
00232 *
00233 *     End of CLARZB
00234 *
00235       END
 All Files Functions