001:       SUBROUTINE SLARZB( 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:       REAL               C( LDC, * ), T( LDT, * ), V( LDV, * ),
015:      $                   WORK( LDWORK, * )
016: *     ..
017: *
018: *  Purpose
019: *  =======
020: *
021: *  SLARZB applies a real block reflector H or its transpose H**T to
022: *  a real 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' (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) REAL 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) REAL 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) REAL 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) REAL 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:       REAL               ONE
102:       PARAMETER          ( ONE = 1.0E+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           SCOPY, SGEMM, STRMM, XERBLA
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( 'SLARZB', -INFO )
132:          RETURN
133:       END IF
134: *
135:       IF( LSAME( TRANS, 'N' ) ) THEN
136:          TRANST = 'T'
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 ) = C( 1:k, 1:n )'
146: *
147:          DO 10 J = 1, K
148:             CALL SCOPY( 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: *                        C( m-l+1:m, 1:n )' * V( 1:k, 1:l )'
153: *
154:          IF( L.GT.0 )
155:      $      CALL SGEMM( 'Transpose', 'Transpose', N, K, L, ONE,
156:      $                  C( M-L+1, 1 ), LDC, V, LDV, ONE, WORK, LDWORK )
157: *
158: *        W( 1:n, 1:k ) = W( 1:n, 1:k ) * T'  or  W( 1:m, 1:k ) * T
159: *
160:          CALL STRMM( 'Right', 'Lower', TRANST, 'Non-unit', N, K, ONE, T,
161:      $               LDT, WORK, LDWORK )
162: *
163: *        C( 1:k, 1:n ) = C( 1:k, 1:n ) - W( 1:n, 1:k )'
164: *
165:          DO 30 J = 1, N
166:             DO 20 I = 1, K
167:                C( I, J ) = C( I, J ) - WORK( J, I )
168:    20       CONTINUE
169:    30    CONTINUE
170: *
171: *        C( m-l+1:m, 1:n ) = C( m-l+1:m, 1:n ) - ...
172: *                            V( 1:k, 1:l )' * W( 1:n, 1:k )'
173: *
174:          IF( L.GT.0 )
175:      $      CALL SGEMM( 'Transpose', 'Transpose', L, N, K, -ONE, V, LDV,
176:      $                  WORK, LDWORK, ONE, C( M-L+1, 1 ), LDC )
177: *
178:       ELSE IF( LSAME( SIDE, 'R' ) ) THEN
179: *
180: *        Form  C * H  or  C * H'
181: *
182: *        W( 1:m, 1:k ) = C( 1:m, 1:k )
183: *
184:          DO 40 J = 1, K
185:             CALL SCOPY( M, C( 1, J ), 1, WORK( 1, J ), 1 )
186:    40    CONTINUE
187: *
188: *        W( 1:m, 1:k ) = W( 1:m, 1:k ) + ...
189: *                        C( 1:m, n-l+1:n ) * V( 1:k, 1:l )'
190: *
191:          IF( L.GT.0 )
192:      $      CALL SGEMM( 'No transpose', 'Transpose', M, K, L, ONE,
193:      $                  C( 1, N-L+1 ), LDC, V, LDV, ONE, WORK, LDWORK )
194: *
195: *        W( 1:m, 1:k ) = W( 1:m, 1:k ) * T  or  W( 1:m, 1:k ) * T'
196: *
197:          CALL STRMM( 'Right', 'Lower', TRANS, 'Non-unit', M, K, ONE, T,
198:      $               LDT, WORK, LDWORK )
199: *
200: *        C( 1:m, 1:k ) = C( 1:m, 1:k ) - W( 1:m, 1:k )
201: *
202:          DO 60 J = 1, K
203:             DO 50 I = 1, M
204:                C( I, J ) = C( I, J ) - WORK( I, J )
205:    50       CONTINUE
206:    60    CONTINUE
207: *
208: *        C( 1:m, n-l+1:n ) = C( 1:m, n-l+1:n ) - ...
209: *                            W( 1:m, 1:k ) * V( 1:k, 1:l )
210: *
211:          IF( L.GT.0 )
212:      $      CALL SGEMM( 'No transpose', 'No transpose', M, L, K, -ONE,
213:      $                  WORK, LDWORK, V, LDV, ONE, C( 1, N-L+1 ), LDC )
214: *
215:       END IF
216: *
217:       RETURN
218: *
219: *     End of SLARZB
220: *
221:       END
222: