001:       SUBROUTINE ZLARZB( SIDE, TRANS, DIRECT, STOREV, M, N, K, L, V,
002:      $                   LDV, T, LDT, C, LDC, WORK, LDWORK )
003: *
004: *  -- LAPACK routine (version 3.2) --
005: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
006: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
007: *     November 2006
008: *
009: *     .. Scalar Arguments ..
010:       CHARACTER          DIRECT, SIDE, STOREV, TRANS
011:       INTEGER            K, L, LDC, LDT, LDV, LDWORK, M, N
012: *     ..
013: *     .. Array Arguments ..
014:       COMPLEX*16         C( LDC, * ), T( LDT, * ), V( LDV, * ),
015:      $                   WORK( LDWORK, * )
016: *     ..
017: *
018: *  Purpose
019: *  =======
020: *
021: *  ZLARZB applies a complex block reflector H or its transpose H**H
022: *  to a complex distributed M-by-N  C from the left or the right.
023: *
024: *  Currently, only STOREV = 'R' and DIRECT = 'B' are supported.
025: *
026: *  Arguments
027: *  =========
028: *
029: *  SIDE    (input) CHARACTER*1
030: *          = 'L': apply H or H' from the Left
031: *          = 'R': apply H or H' from the Right
032: *
033: *  TRANS   (input) CHARACTER*1
034: *          = 'N': apply H (No transpose)
035: *          = 'C': apply H' (Conjugate transpose)
036: *
037: *  DIRECT  (input) CHARACTER*1
038: *          Indicates how H is formed from a product of elementary
039: *          reflectors
040: *          = 'F': H = H(1) H(2) . . . H(k) (Forward, not supported yet)
041: *          = 'B': H = H(k) . . . H(2) H(1) (Backward)
042: *
043: *  STOREV  (input) CHARACTER*1
044: *          Indicates how the vectors which define the elementary
045: *          reflectors are stored:
046: *          = 'C': Columnwise                        (not supported yet)
047: *          = 'R': Rowwise
048: *
049: *  M       (input) INTEGER
050: *          The number of rows of the matrix C.
051: *
052: *  N       (input) INTEGER
053: *          The number of columns of the matrix C.
054: *
055: *  K       (input) INTEGER
056: *          The order of the matrix T (= the number of elementary
057: *          reflectors whose product defines the block reflector).
058: *
059: *  L       (input) INTEGER
060: *          The number of columns of the matrix V containing the
061: *          meaningful part of the Householder reflectors.
062: *          If SIDE = 'L', M >= L >= 0, if SIDE = 'R', N >= L >= 0.
063: *
064: *  V       (input) COMPLEX*16 array, dimension (LDV,NV).
065: *          If STOREV = 'C', NV = K; if STOREV = 'R', NV = L.
066: *
067: *  LDV     (input) INTEGER
068: *          The leading dimension of the array V.
069: *          If STOREV = 'C', LDV >= L; if STOREV = 'R', LDV >= K.
070: *
071: *  T       (input) COMPLEX*16 array, dimension (LDT,K)
072: *          The triangular K-by-K matrix T in the representation of the
073: *          block reflector.
074: *
075: *  LDT     (input) INTEGER
076: *          The leading dimension of the array T. LDT >= K.
077: *
078: *  C       (input/output) COMPLEX*16 array, dimension (LDC,N)
079: *          On entry, the M-by-N matrix C.
080: *          On exit, C is overwritten by H*C or H'*C or C*H or C*H'.
081: *
082: *  LDC     (input) INTEGER
083: *          The leading dimension of the array C. LDC >= max(1,M).
084: *
085: *  WORK    (workspace) COMPLEX*16 array, dimension (LDWORK,K)
086: *
087: *  LDWORK  (input) INTEGER
088: *          The leading dimension of the array WORK.
089: *          If SIDE = 'L', LDWORK >= max(1,N);
090: *          if SIDE = 'R', LDWORK >= max(1,M).
091: *
092: *  Further Details
093: *  ===============
094: *
095: *  Based on contributions by
096: *    A. Petitet, Computer Science Dept., Univ. of Tenn., Knoxville, USA
097: *
098: *  =====================================================================
099: *
100: *     .. Parameters ..
101:       COMPLEX*16         ONE
102:       PARAMETER          ( ONE = ( 1.0D+0, 0.0D+0 ) )
103: *     ..
104: *     .. Local Scalars ..
105:       CHARACTER          TRANST
106:       INTEGER            I, INFO, J
107: *     ..
108: *     .. External Functions ..
109:       LOGICAL            LSAME
110:       EXTERNAL           LSAME
111: *     ..
112: *     .. External Subroutines ..
113:       EXTERNAL           XERBLA, ZCOPY, ZGEMM, ZLACGV, ZTRMM
114: *     ..
115: *     .. Executable Statements ..
116: *
117: *     Quick return if possible
118: *
119:       IF( M.LE.0 .OR. N.LE.0 )
120:      $   RETURN
121: *
122: *     Check for currently supported options
123: *
124:       INFO = 0
125:       IF( .NOT.LSAME( DIRECT, 'B' ) ) THEN
126:          INFO = -3
127:       ELSE IF( .NOT.LSAME( STOREV, 'R' ) ) THEN
128:          INFO = -4
129:       END IF
130:       IF( INFO.NE.0 ) THEN
131:          CALL XERBLA( 'ZLARZB', -INFO )
132:          RETURN
133:       END IF
134: *
135:       IF( LSAME( TRANS, 'N' ) ) THEN
136:          TRANST = 'C'
137:       ELSE
138:          TRANST = 'N'
139:       END IF
140: *
141:       IF( LSAME( SIDE, 'L' ) ) THEN
142: *
143: *        Form  H * C  or  H' * C
144: *
145: *        W( 1:n, 1:k ) = conjg( C( 1:k, 1:n )' )
146: *
147:          DO 10 J = 1, K
148:             CALL ZCOPY( N, C( J, 1 ), LDC, WORK( 1, J ), 1 )
149:    10    CONTINUE
150: *
151: *        W( 1:n, 1:k ) = W( 1:n, 1:k ) + ...
152: *                        conjg( C( m-l+1:m, 1:n )' ) * V( 1:k, 1:l )'
153: *
154:          IF( L.GT.0 )
155:      $      CALL ZGEMM( 'Transpose', 'Conjugate transpose', N, K, L,
156:      $                  ONE, C( M-L+1, 1 ), LDC, V, LDV, ONE, WORK,
157:      $                  LDWORK )
158: *
159: *        W( 1:n, 1:k ) = W( 1:n, 1:k ) * T'  or  W( 1:m, 1:k ) * T
160: *
161:          CALL ZTRMM( 'Right', 'Lower', TRANST, 'Non-unit', N, K, ONE, T,
162:      $               LDT, WORK, LDWORK )
163: *
164: *        C( 1:k, 1:n ) = C( 1:k, 1:n ) - conjg( W( 1:n, 1:k )' )
165: *
166:          DO 30 J = 1, N
167:             DO 20 I = 1, K
168:                C( I, J ) = C( I, J ) - WORK( J, I )
169:    20       CONTINUE
170:    30    CONTINUE
171: *
172: *        C( m-l+1:m, 1:n ) = C( m-l+1:m, 1:n ) - ...
173: *                    conjg( V( 1:k, 1:l )' ) * conjg( W( 1:n, 1:k )' )
174: *
175:          IF( L.GT.0 )
176:      $      CALL ZGEMM( 'Transpose', 'Transpose', L, N, K, -ONE, V, LDV,
177:      $                  WORK, LDWORK, ONE, C( M-L+1, 1 ), LDC )
178: *
179:       ELSE IF( LSAME( SIDE, 'R' ) ) THEN
180: *
181: *        Form  C * H  or  C * H'
182: *
183: *        W( 1:m, 1:k ) = C( 1:m, 1:k )
184: *
185:          DO 40 J = 1, K
186:             CALL ZCOPY( M, C( 1, J ), 1, WORK( 1, J ), 1 )
187:    40    CONTINUE
188: *
189: *        W( 1:m, 1:k ) = W( 1:m, 1:k ) + ...
190: *                        C( 1:m, n-l+1:n ) * conjg( V( 1:k, 1:l )' )
191: *
192:          IF( L.GT.0 )
193:      $      CALL ZGEMM( 'No transpose', 'Transpose', M, K, L, ONE,
194:      $                  C( 1, N-L+1 ), LDC, V, LDV, ONE, WORK, LDWORK )
195: *
196: *        W( 1:m, 1:k ) = W( 1:m, 1:k ) * conjg( T )  or
197: *                        W( 1:m, 1:k ) * conjg( T' )
198: *
199:          DO 50 J = 1, K
200:             CALL ZLACGV( K-J+1, T( J, J ), 1 )
201:    50    CONTINUE
202:          CALL ZTRMM( 'Right', 'Lower', TRANS, 'Non-unit', M, K, ONE, T,
203:      $               LDT, WORK, LDWORK )
204:          DO 60 J = 1, K
205:             CALL ZLACGV( K-J+1, T( J, J ), 1 )
206:    60    CONTINUE
207: *
208: *        C( 1:m, 1:k ) = C( 1:m, 1:k ) - W( 1:m, 1:k )
209: *
210:          DO 80 J = 1, K
211:             DO 70 I = 1, M
212:                C( I, J ) = C( I, J ) - WORK( I, J )
213:    70       CONTINUE
214:    80    CONTINUE
215: *
216: *        C( 1:m, n-l+1:n ) = C( 1:m, n-l+1:n ) - ...
217: *                            W( 1:m, 1:k ) * conjg( V( 1:k, 1:l ) )
218: *
219:          DO 90 J = 1, L
220:             CALL ZLACGV( K, V( 1, J ), 1 )
221:    90    CONTINUE
222:          IF( L.GT.0 )
223:      $      CALL ZGEMM( 'No transpose', 'No transpose', M, L, K, -ONE,
224:      $                  WORK, LDWORK, V, LDV, ONE, C( 1, N-L+1 ), LDC )
225:          DO 100 J = 1, L
226:             CALL ZLACGV( K, V( 1, J ), 1 )
227:   100    CONTINUE
228: *
229:       END IF
230: *
231:       RETURN
232: *
233: *     End of ZLARZB
234: *
235:       END
236: