1      SUBROUTINE MB01RX( SIDE, UPLO, TRANS, M, N, ALPHA, BETA, R, LDR,
2     $                   A, LDA, B, LDB, INFO )
3C
4C     RELEASE 4.0, WGS COPYRIGHT 1999.
5C
6C     PURPOSE
7C
8C     To compute either the upper or lower triangular part of one of the
9C     matrix formulas
10C        _
11C        R = alpha*R + beta*op( A )*B,                               (1)
12C        _
13C        R = alpha*R + beta*B*op( A ),                               (2)
14C                                             _
15C     where alpha and beta are scalars, R and R are m-by-m matrices,
16C     op( A ) and B are m-by-n and n-by-m matrices for (1), or n-by-m
17C     and m-by-n matrices for (2), respectively, and op( A ) is one of
18C
19C        op( A ) = A   or   op( A ) = A',  the transpose of A.
20C
21C     The result is overwritten on R.
22C
23C     ARGUMENTS
24C
25C     Mode Parameters
26C
27C     SIDE    CHARACTER*1
28C             Specifies whether the matrix A appears on the left or
29C             right in the matrix product as follows:
30C                     _
31C             = 'L':  R = alpha*R + beta*op( A )*B;
32C                     _
33C             = 'R':  R = alpha*R + beta*B*op( A ).
34C
35C     UPLO    CHARACTER*1                               _
36C             Specifies which triangles of the matrices R and R are
37C             computed and given, respectively, as follows:
38C             = 'U':  the upper triangular part;
39C             = 'L':  the lower triangular part.
40C
41C     TRANS   CHARACTER*1
42C             Specifies the form of op( A ) to be used in the matrix
43C             multiplication as follows:
44C             = 'N':  op( A ) = A;
45C             = 'T':  op( A ) = A';
46C             = 'C':  op( A ) = A'.
47C
48C     Input/Output Parameters
49C
50C     M       (input) INTEGER           _
51C             The order of the matrices R and R, the number of rows of
52C             the matrix op( A ) and the number of columns of the
53C             matrix B, for SIDE = 'L', or the number of rows of the
54C             matrix B and the number of columns of the matrix op( A ),
55C             for SIDE = 'R'.  M >= 0.
56C
57C     N       (input) INTEGER
58C             The number of rows of the matrix B and the number of
59C             columns of the matrix op( A ), for SIDE = 'L', or the
60C             number of rows of the matrix op( A ) and the number of
61C             columns of the matrix B, for SIDE = 'R'.  N >= 0.
62C
63C     ALPHA   (input) DOUBLE PRECISION
64C             The scalar alpha. When alpha is zero then R need not be
65C             set before entry.
66C
67C     BETA    (input) DOUBLE PRECISION
68C             The scalar beta. When beta is zero then A and B are not
69C             referenced.
70C
71C     R       (input/output) DOUBLE PRECISION array, dimension (LDR,M)
72C             On entry with UPLO = 'U', the leading M-by-M upper
73C             triangular part of this array must contain the upper
74C             triangular part of the matrix R; the strictly lower
75C             triangular part of the array is not referenced.
76C             On entry with UPLO = 'L', the leading M-by-M lower
77C             triangular part of this array must contain the lower
78C             triangular part of the matrix R; the strictly upper
79C             triangular part of the array is not referenced.
80C             On exit, the leading M-by-M upper triangular part (if
81C             UPLO = 'U'), or lower triangular part (if UPLO = 'L') of
82C             this array contains the corresponding triangular part of
83C                                 _
84C             the computed matrix R.
85C
86C     LDR     INTEGER
87C             The leading dimension of array R.  LDR >= MAX(1,M).
88C
89C     A       (input) DOUBLE PRECISION array, dimension (LDA,k), where
90C             k = N  when  SIDE = 'L', and TRANS = 'N', or
91C                          SIDE = 'R', and TRANS = 'T';
92C             k = M  when  SIDE = 'R', and TRANS = 'N', or
93C                          SIDE = 'L', and TRANS = 'T'.
94C             On entry, if SIDE = 'L', and TRANS = 'N', or
95C                          SIDE = 'R', and TRANS = 'T',
96C             the leading M-by-N part of this array must contain the
97C             matrix A.
98C             On entry, if SIDE = 'R', and TRANS = 'N', or
99C                          SIDE = 'L', and TRANS = 'T',
100C             the leading N-by-M part of this array must contain the
101C             matrix A.
102C
103C     LDA     INTEGER
104C             The leading dimension of array A.  LDA >= MAX(1,l), where
105C             l = M  when  SIDE = 'L', and TRANS = 'N', or
106C                          SIDE = 'R', and TRANS = 'T';
107C             l = N  when  SIDE = 'R', and TRANS = 'N', or
108C                          SIDE = 'L', and TRANS = 'T'.
109C
110C     B       (input) DOUBLE PRECISION array, dimension (LDB,p), where
111C             p = M  when  SIDE = 'L';
112C             p = N  when  SIDE = 'R'.
113C             On entry, the leading N-by-M part, if SIDE = 'L', or
114C             M-by-N part, if SIDE = 'R', of this array must contain the
115C             matrix B.
116C
117C     LDB     INTEGER
118C             The leading dimension of array B.
119C             LDB >= MAX(1,N), if SIDE = 'L';
120C             LDB >= MAX(1,M), if SIDE = 'R'.
121C
122C     Error Indicator
123C
124C     INFO    INTEGER
125C             = 0:  successful exit;
126C             < 0:  if INFO = -i, the i-th argument had an illegal
127C                   value.
128C
129C     METHOD
130C
131C     The matrix expression is evaluated taking the triangular
132C     structure into account. BLAS 2 operations are used. A block
133C     algorithm can be easily constructed; it can use BLAS 3 GEMM
134C     operations for most computations, and calls of this BLAS 2
135C     algorithm for computing the triangles.
136C
137C     FURTHER COMMENTS
138C
139C     The main application of this routine is when the result should
140C     be a symmetric matrix, e.g., when B = X*op( A )', for (1), or
141C     B = op( A )'*X, for (2), where B is already available and X = X'.
142C
143C     CONTRIBUTORS
144C
145C     V. Sima, Katholieke Univ. Leuven, Belgium, Feb. 1999.
146C
147C     REVISIONS
148C
149C     -
150C
151C     KEYWORDS
152C
153C     Elementary matrix operations, matrix algebra, matrix operations.
154C
155C     ******************************************************************
156C
157C     .. Parameters ..
158      DOUBLE PRECISION  ZERO, ONE
159      PARAMETER         ( ZERO = 0.0D0, ONE = 1.0D0 )
160C     .. Scalar Arguments ..
161      CHARACTER         SIDE, TRANS, UPLO
162      INTEGER           INFO, LDA, LDB, LDR, M, N
163      DOUBLE PRECISION  ALPHA, BETA
164C     .. Array Arguments ..
165      DOUBLE PRECISION  A(LDA,*), B(LDB,*), R(LDR,*)
166C     .. Local Scalars ..
167      LOGICAL           LSIDE, LTRANS, LUPLO
168      INTEGER           J
169C     .. External Functions ..
170      LOGICAL           LSAME
171      EXTERNAL          LSAME
172C     .. External Subroutines ..
173      EXTERNAL          DGEMV, DLASCL, DLASET, XERBLA
174C     .. Intrinsic Functions ..
175      INTRINSIC         MAX
176C     .. Executable Statements ..
177C
178C     Test the input scalar arguments.
179C
180      INFO   = 0
181      LSIDE  = LSAME( SIDE,  'L' )
182      LUPLO  = LSAME( UPLO,  'U' )
183      LTRANS = LSAME( TRANS, 'T' ) .OR. LSAME( TRANS, 'C' )
184C
185      IF(      ( .NOT.LSIDE  ).AND.( .NOT.LSAME( SIDE,  'R' ) ) )THEN
186         INFO = -1
187      ELSE IF( ( .NOT.LUPLO  ).AND.( .NOT.LSAME( UPLO,  'L' ) ) )THEN
188         INFO = -2
189      ELSE IF( ( .NOT.LTRANS ).AND.( .NOT.LSAME( TRANS, 'N' ) ) )THEN
190         INFO = -3
191      ELSE IF( M.LT.0 ) THEN
192         INFO = -4
193      ELSE IF( N.LT.0 ) THEN
194         INFO = -5
195      ELSE IF( LDR.LT.MAX( 1, M ) ) THEN
196         INFO = -9
197      ELSE IF( LDA.LT.1 .OR.
198     $   ( ( (      LSIDE .AND. .NOT.LTRANS ) .OR.
199     $       ( .NOT.LSIDE .AND.      LTRANS ) ) .AND. LDA.LT.M ) .OR.
200     $   ( ( (      LSIDE .AND.      LTRANS ) .OR.
201     $       ( .NOT.LSIDE .AND. .NOT.LTRANS ) ) .AND. LDA.LT.N ) ) THEN
202         INFO = -11
203      ELSE IF( LDB.LT.1 .OR.
204     $       (      LSIDE .AND. LDB.LT.N ) .OR.
205     $       ( .NOT.LSIDE .AND. LDB.LT.M ) ) THEN
206         INFO = -13
207      END IF
208C
209      IF ( INFO.NE.0 ) THEN
210C
211C        Error return.
212C
213         CALL XERBLA( 'MB01RX', -INFO )
214         RETURN
215      END IF
216C
217C     Quick return if possible.
218C
219      IF ( M.EQ.0 )
220     $   RETURN
221C
222      IF ( BETA.EQ.ZERO ) THEN
223         IF ( ALPHA.EQ.ZERO ) THEN
224C
225C           Special case when both alpha = 0 and beta = 0.
226C
227            CALL DLASET( UPLO, M, M, ZERO, ZERO, R, LDR )
228         ELSE
229C
230C           Special case beta = 0.
231C
232            IF ( ALPHA.NE.ONE )
233     $         CALL DLASCL( UPLO, 0, 0, ONE, ALPHA, M, M, R, LDR, INFO )
234         END IF
235         RETURN
236      END IF
237C
238      IF ( N.EQ.0 )
239     $   RETURN
240C
241C     General case: beta <> 0.
242C     Compute the required triangle of (1) or (2) using BLAS 2
243C     operations.
244C
245      IF( LSIDE ) THEN
246         IF( LUPLO ) THEN
247            IF ( LTRANS ) THEN
248               DO 10 J = 1, M
249                  CALL DGEMV( TRANS, N, J, BETA, A, LDA, B(1,J), 1,
250     $                        ALPHA, R(1,J), 1 )
251   10          CONTINUE
252            ELSE
253               DO 20 J = 1, M
254                  CALL DGEMV( TRANS, J, N, BETA, A, LDA, B(1,J), 1,
255     $                        ALPHA, R(1,J), 1 )
256   20          CONTINUE
257            END IF
258         ELSE
259            IF ( LTRANS ) THEN
260               DO 30 J = 1, M
261                  CALL DGEMV( TRANS, N, M-J+1, BETA, A(1,J), LDA,
262     $                        B(1,J), 1, ALPHA, R(J,J), 1 )
263   30          CONTINUE
264            ELSE
265               DO 40 J = 1, M
266                  CALL DGEMV( TRANS, M-J+1, N, BETA, A(J,1), LDA,
267     $                        B(1,J), 1, ALPHA, R(J,J), 1 )
268   40          CONTINUE
269            END IF
270         END IF
271C
272      ELSE
273         IF( LUPLO ) THEN
274            IF( LTRANS ) THEN
275               DO 50 J = 1, M
276                  CALL DGEMV( 'NoTranspose', J, N, BETA, B, LDB, A(J,1),
277     $                        LDA, ALPHA, R(1,J), 1 )
278   50          CONTINUE
279            ELSE
280               DO 60 J = 1, M
281                  CALL DGEMV( 'NoTranspose', J, N, BETA, B, LDB, A(1,J),
282     $                        1, ALPHA, R(1,J), 1 )
283   60          CONTINUE
284            END IF
285         ELSE
286            IF( LTRANS ) THEN
287               DO 70 J = 1, M
288                  CALL DGEMV( 'NoTranspose', M-J+1, N, BETA, B(J,1),
289     $                        LDB, A(J,1), LDA, ALPHA, R(J,J), 1 )
290   70           CONTINUE
291            ELSE
292               DO 80 J = 1, M
293                  CALL DGEMV( 'NoTranspose', M-J+1, N, BETA, B(J,1),
294     $                        LDB, A(1,J), 1, ALPHA, R(J,J), 1 )
295   80          CONTINUE
296            END IF
297         END IF
298      END IF
299C
300      RETURN
301C *** Last line of MB01RX ***
302      END
303