001:       SUBROUTINE CORE_SSSRFB( SIDE, TRANS, DIRECT, STOREV, M1, M2, NN,
002:      $                       IB, A1, LDA1, A2, LDA2, V, LDV, T, LDT,
003:      $                       WORK, LDWORK, INFO )  
004: 
005: *********************************************************************
006: *     PLASMA core_blas routine (version 2.1.0)                      *
007: *     Author: Hatem Ltaief                                          *
008: *     Release Date: November, 15th 2009                             *
009: *     PLASMA is a software package provided by Univ. of Tennessee,  *
010: *     Univ. of California Berkeley and Univ. of Colorado Denver.    *
011: *********************************************************************
012: *
013: *     .. Scalar Arguments ..
014:       INTEGER            M1, M2, NN, IB
015:       INTEGER            LDA1, LDA2, LDV, LDT, LDWORK
016:       INTEGER            INFO
017: *     ..
018: *     .. Array Arguments ..
019:       REAL            A1(LDA1, * ), A2( LDA2, * )
020:       REAL            V( LDV, * ), T( LDT, * )
021:       REAL            WORK( LDWORK, * )
022: *     ..
023: *     .. Character Arguments ..
024:       CHARACTER          SIDE, TRANS, DIRECT, STOREV
025: *     ..
026: *
027: *  Purpose
028: *  =======
029: *
030: *  CORE_SSSRFB applies a real block reflector H or its transpose H' to a
031: *  real rectangular matrix formed by coupling two tiles A1 and A2,
032: *  from either the left or the right.
033: *
034: *  Arguments
035: *  =========
036: *
037: *  SIDE    (input) CHARACTER*1
038: *          = 'L': apply H or H' from the Left
039: *          = 'R': apply H or H' from the Right
040: *
041: *  TRANS   (input) CHARACTER*1
042: *          = 'N': apply H (No transpose)
043: *          = 'C': apply H' (Conjugate transpose)
044: *
045: *  DIRECT  (input) CHARACTER*1
046: *          Indicates how H is formed from a product of elementary
047: *          reflectors
048: *          = 'F': H = H(1) H(2) . . . H(k) (Forward)
049: *          = 'B': H = H(k) . . . H(2) H(1) (Backward)
050: *
051: *  STOREV  (input) CHARACTER*1
052: *          Indicates how the vectors which define the elementary
053: *          reflectors are stored:
054: *          = 'C': Columnwise
055: *          = 'R': Rowwise
056: *
057: *  M1      (input) INTEGER
058: *          The number of rows of the tile A1. M1 >= 0.
059: *
060: *  M2      (input) INTEGER
061: *          The number of rows of the tile A2. M2 >= 0.
062: *
063: *  NN      (input) INTEGER
064: *          The number of columns of the tiles A1 and A2. NN >= 0.
065: *
066: *  IB      (input) INTEGER
067: *          The inner-blocking size.  IB >= 0.
068: *
069: *  K       (input) INTEGER
070: *          The order of the matrix T (= the number of elementary
071: *          reflectors whose product defines the block reflector).
072: *
073: *  A1      (input/output) REAL array, dimension (LDA1,M1)
074: *          On entry, the M1-by-NN tile A1.
075: *          On exit, A1 is overwritten by the application of Q.
076: *
077: *  LDA1    (input) INTEGER
078: *          The leading dimension of the tile A1. LDA1 >= max(1,M1).
079: *
080: *  A2      (input/output) REAL array, dimension (LDA2,M2)
081: *          On entry, the M2-by-NN tile A2.
082: *          On exit, A2 is overwritten by the application of Q.
083: *
084: *  LDA2    (input) INTEGER
085: *          The leading dimension of the tile A2. LDA2 >= max(1,M2).
086: *
087: *  V       (input) REAL array, dimension
088: *                                (LDV,K) if STOREV = 'C'
089: *                                (LDV,M) if STOREV = 'R' and SIDE = 'L'
090: *                                (LDV,N) if STOREV = 'R' and SIDE = 'R'
091: *          The matrix V. See further details.
092: *
093: *  LDV     (input) INTEGER
094: *          The leading dimension of the array V.
095: *          If STOREV = 'C' and SIDE = 'L', LDV >= max(1,M);
096: *          if STOREV = 'C' and SIDE = 'R', LDV >= max(1,N);
097: *          if STOREV = 'R', LDV >= K.
098: *
099: *  T       (input) REAL array, dimension (LDT,K)
100: *          The triangular K-by-K matrix T in the representation of the
101: *          block reflector.
102: *
103: *  LDT     (input) INTEGER
104: *          The leading dimension of the array T. LDT >= K.
105: *
106: *  WORK    (workspace/output) REAL array, dimension (LDWORK,K)
107: *
108: *  LDWORK  (input) INTEGER
109: *          The dimension of the array WORK.
110: *
111: *  INFO    (output) INTEGER
112: *          = 0:  successful exit
113: *          < 0:  if INFO = -i, the i-th argument had an illegal value
114: *
115: *  =====================================================================
116: *     ..
117: *     .. Local Scalars ..
118:       CHARACTER          TRANST
119:       INTEGER            J
120: *     ..
121: *     .. Parameters ..
122:       REAL            SONE, MSONE
123:       INTEGER            IONE
124:       PARAMETER          ( SONE = 1.0E+0 )
125:       PARAMETER          ( MSONE = -1.0E+0 )
126:       PARAMETER          ( IONE = 1 )
127: *     ..
128: *     .. External Functions ..
129:       LOGICAL            LSAME
130:       EXTERNAL           LSAME
131: *     ..
132: *     Test the input arguments.
133:       INFO = 0
134:       IF( M1.LT.0 ) THEN
135:          INFO = -5
136:       ELSE IF( M2.LT.0 ) THEN
137:          INFO = -6
138:       ELSE IF( NN.LT.0 ) THEN
139:          INFO = -7
140:       ELSE IF( IB.LT.0 ) THEN
141:          INFO = -8
142:       END IF
143:       IF( INFO.NE.0 ) THEN
144:          CALL XERBLA( 'CORE_SSSRFB', -INFO )
145:          RETURN
146:       END IF
147: *
148: *     Quick return if possible.
149: *
150:       IF( M1.LE.0 .OR. M2.LE.0 .OR. NN.LE.0 .OR. IB.LE.0 )
151:      $   RETURN
152: *
153:       IF( LSAME( STOREV, 'C' ) ) THEN
154: *
155:          IF( LSAME( SIDE, 'L' ) ) THEN
156: *
157: *              B = A1+V'*A2
158: *
159:                CALL SLACPY( 'General', IB, NN, A1, LDA1,
160:      $                     WORK, LDWORK )
161:                CALL SGEMM( 'Transpose', 'Notranspose', IB, NN, M2,
162:      $                    SONE, V, LDV, A2, LDA2, SONE,
163:      $                    WORK, LDWORK )
164: *
165: *              A2 = A2 - V*T*B --->  B=T*B, A2=A2-V*B
166: *
167:                CALL STRMM( 'Left', 'Upper', TRANS, 'Nounit', IB, NN, 
168:      $                    SONE, T, LDT, WORK, LDWORK )
169:                CALL SGEMM( 'Notranspose', 'Notranspose', M2, NN, IB,
170:      $                    MSONE, V, LDV, WORK, LDWORK, SONE,
171:      $                    A2, LDA2 )
172: *
173: *              A1 = A1 - T*B
174: *
175:                DO 20 J = 1, NN
176:                   CALL SAXPY( IB, MSONE, WORK( 1, J ), IONE,
177:      $                       A1( 1, J ), IONE )
178:  20            CONTINUE
179: *
180:          ELSE IF( LSAME( SIDE, 'R' ) ) THEN
181: *
182: *              B = A1+V'*A2
183: *
184:                CALL SLACPY( 'General', IB, NN, A1, LDA1,
185:      $                     WORK, LDWORK)
186:                CALL SGEMM( 'Transpose', 'Notranspose', IB, NN, M2,
187:      $                    SONE, V, LDV, A2, LDA2, SONE,
188:      $                    WORK, LDWORK )
189: *
190: *              A2 = A2 - V*T*B --->  B=T*B, A2=A2-V*B
191: *
192:                CALL STRMM( 'Left', 'Upper', TRANS, 'Nounit', IB, NN, 
193:      $                    SONE, T, LDT, WORK, LDWORK )
194:                CALL SGEMM( 'Notranspose', 'Notranspose', M2, NN, IB,
195:      $                    MSONE, V, LDV, WORK, LDWORK, SONE,
196:      $                    A2, LDA2 )
197: *
198: *              A1 = A1 - T*B
199: *
200:                DO 50 J = 1, NN
201:                   CALL SAXPY( IB, MSONE, WORK( 1, J ), IONE,
202:      $                       A1( J, 1 ), IONE )
203:  50            CONTINUE
204:          ENDIF
205: *
206:       ELSE IF( LSAME( STOREV, 'R' ) ) THEN
207: *
208:          IF( LSAME( DIRECT, 'F' ) ) THEN
209: *
210:             IF( LSAME( SIDE, 'L' ) ) THEN
211: *
212: *                 B = A1+V'*A2
213: *
214:                   CALL SLACPY( 'General', IB, NN, A1,
215:      $                        LDA1, WORK, LDWORK)
216:                   CALL SGEMM( 'NoTranspose', 'Notranspose', IB, NN, M2,
217:      $                       SONE, V, LDV, A2, LDA2,
218:      $                       SONE, WORK, LDWORK )
219: *
220: *                 A2 = A2 - V*T*B --->  B=T*B, A2=A2-V*B
221: *
222:                   CALL STRMM( 'Left', 'Up', TRANS, 'Nounit', IB,
223:      $                       NN, SONE, T, LDT, WORK, LDWORK )
224:                   CALL SGEMM( 'Transpose', 'Notranspose', M2, NN,
225:      $                       IB, MSONE, V, LDV, WORK,
226:      $                       LDWORK, SONE, A2, LDA2 )
227: *
228: *                 A1 = A1 - T*B
229: *
230:                   DO 80 J = 1, NN
231:                      CALL SAXPY( IB, MSONE, WORK( 1, J ), IONE,
232:      $                          A1( 1, J ), IONE )
233:  80               CONTINUE
234: *
235:             ELSE IF( LSAME( SIDE, 'R' ) ) THEN
236: *
237: *                      B = A1+A2*V
238: *
239:                        CALL SLACPY( 'General', NN, IB, A1,
240:      $                             LDA1, WORK, LDWORK )
241:                        CALL SGEMM( 'Notranspose', 'Transpose', NN,
242:      $                            IB, M2, SONE, A2, LDA2, V,
243:      $                            LDV, SONE, WORK, LDWORK )
244: *
245: *                      A2 = A2 - B*T*V' --->  B=B*T, A2=A2-B*V'
246: *
247:                        CALL STRMM( 'Right', 'Up', TRANS, 
248:      $                            'Nounit', NN, IB, SONE,
249:      $                            T, LDT, WORK, LDWORK )
250:                        CALL SGEMM( 'Notranspose', 'Notranspose', NN, M2,
251:      $                            IB, MSONE, WORK, LDWORK, V,
252:      $                            LDV, SONE, A2, LDA2 )
253: *
254: *                      A1 = A1 - B*T
255: *
256:                        DO 110 J = 1, NN
257:                           CALL SAXPY( IB, MSONE, WORK( J, 1 ), LDWORK,
258:      $                               A1( J, 1 ), LDA1 )
259:  110                   CONTINUE
260: *
261:             ENDIF
262: *
263:          ELSE
264: *
265:             IF( LSAME( SIDE, 'L' ) ) THEN
266: *
267: *                 B = A1+V'*A2
268: *
269:                   CALL SLACPY( 'General', IB, NN, A1,
270:      $                        LDA1, WORK, LDWORK)
271:                   CALL SGEMM( 'Transpose', 'Notranspose', IB, NN,
272:      $                       M2, SONE, V, LDV, A2, LDA2,
273:      $                       SONE, WORK, LDWORK )
274: *
275: *                 A2 = A2 - V*T*B --->  B=T*B, A2=A2-V*B
276: *
277:                   CALL STRMM( 'Left', 'Up', TRANS, 'Nounit', IB,
278:      $                       NN, SONE, T, LDT, WORK, IB )
279:                   CALL SGEMM( 'Notranspose', 'Notranspose', M2, NN,
280:      $                       IB, MSONE, V, LDV, WORK,
281:      $                       LDWORK, SONE, A2, LDA2 )
282: *
283: *                 A1 = A1 - T*B
284: *
285:                   DO 140 J = 1, NN
286:                      CALL SAXPY( IB, MSONE, WORK( 1, J ), IONE,
287:      $                          A1( 1, J ), IONE )
288:  140              CONTINUE
289: *
290:             ELSE IF( LSAME( SIDE, 'R' ) ) THEN
291: *
292: *                 B = A1+A2*V
293: *
294:                   CALL SLACPY( 'General', NN, IB, A1,
295:      $                        LDA1, WORK, LDWORK )
296:                   CALL SGEMM( 'Notranspose', 'Transpose', NN, IB,
297:      $                       M2, SONE, A2, LDA2, V,
298:      $                       LDV, SONE, WORK, LDWORK )
299: *
300: *                 A2 = A2 - B*T*V' --->  B=B*T, A2=A2-B*V'
301: *
302:                   CALL STRMM( 'Right', 'Up', TRANS, 
303:      $                       'Nounit', NN, IB, SONE,
304:      $                       T, LDT, WORK, LDWORK )
305:                   CALL SGEMM( 'Notranspose', 'Notranspose', NN, M2,
306:      $                       IB, MSONE, WORK, LDWORK, V,
307:      $                       LDV, SONE, A2, LDA2 )
308: *
309: *                 A1 = A1 - B*T
310: *
311:                   DO 170 J = 1, NN
312:                      CALL SAXPY( IB, MSONE, WORK( J, 1 ), LDWORK,
313:      $                          A1( J, 1 ), LDA1 )
314:  170              CONTINUE
315:             ENDIF
316: *
317:          ENDIF
318: *
319:       ENDIF
320: *
321: *     End of CORE_SSSRFB.
322: *
323:       END
324: