1 SUBROUTINE CSYMM(SIDE,UPLO,M,N,ALPHA,A,LDA,B,LDB,BETA,C,LDC)
2 * .. Scalar Arguments ..
3 COMPLEX ALPHA,BETA
4 INTEGER LDA,LDB,LDC,M,N
5 CHARACTER SIDE,UPLO
6 * ..
7 * .. Array Arguments ..
8 COMPLEX A(LDA,*),B(LDB,*),C(LDC,*)
9 * ..
10 *
11 * Purpose
12 * =======
13 *
14 * CSYMM performs one of the matrix-matrix operations
15 *
16 * C := alpha*A*B + beta*C,
17 *
18 * or
19 *
20 * C := alpha*B*A + beta*C,
21 *
22 * where alpha and beta are scalars, A is a symmetric matrix and B and
23 * C are m by n matrices.
24 *
25 * Arguments
26 * ==========
27 *
28 * SIDE - CHARACTER*1.
29 * On entry, SIDE specifies whether the symmetric matrix A
30 * appears on the left or right in the operation as follows:
31 *
32 * SIDE = 'L' or 'l' C := alpha*A*B + beta*C,
33 *
34 * SIDE = 'R' or 'r' C := alpha*B*A + beta*C,
35 *
36 * Unchanged on exit.
37 *
38 * UPLO - CHARACTER*1.
39 * On entry, UPLO specifies whether the upper or lower
40 * triangular part of the symmetric matrix A is to be
41 * referenced as follows:
42 *
43 * UPLO = 'U' or 'u' Only the upper triangular part of the
44 * symmetric matrix is to be referenced.
45 *
46 * UPLO = 'L' or 'l' Only the lower triangular part of the
47 * symmetric matrix is to be referenced.
48 *
49 * Unchanged on exit.
50 *
51 * M - INTEGER.
52 * On entry, M specifies the number of rows of the matrix C.
53 * M must be at least zero.
54 * Unchanged on exit.
55 *
56 * N - INTEGER.
57 * On entry, N specifies the number of columns of the matrix C.
58 * N must be at least zero.
59 * Unchanged on exit.
60 *
61 * ALPHA - COMPLEX .
62 * On entry, ALPHA specifies the scalar alpha.
63 * Unchanged on exit.
64 *
65 * A - COMPLEX array of DIMENSION ( LDA, ka ), where ka is
66 * m when SIDE = 'L' or 'l' and is n otherwise.
67 * Before entry with SIDE = 'L' or 'l', the m by m part of
68 * the array A must contain the symmetric matrix, such that
69 * when UPLO = 'U' or 'u', the leading m by m upper triangular
70 * part of the array A must contain the upper triangular part
71 * of the symmetric matrix and the strictly lower triangular
72 * part of A is not referenced, and when UPLO = 'L' or 'l',
73 * the leading m by m lower triangular part of the array A
74 * must contain the lower triangular part of the symmetric
75 * matrix and the strictly upper triangular part of A is not
76 * referenced.
77 * Before entry with SIDE = 'R' or 'r', the n by n part of
78 * the array A must contain the symmetric matrix, such that
79 * when UPLO = 'U' or 'u', the leading n by n upper triangular
80 * part of the array A must contain the upper triangular part
81 * of the symmetric matrix and the strictly lower triangular
82 * part of A is not referenced, and when UPLO = 'L' or 'l',
83 * the leading n by n lower triangular part of the array A
84 * must contain the lower triangular part of the symmetric
85 * matrix and the strictly upper triangular part of A is not
86 * referenced.
87 * Unchanged on exit.
88 *
89 * LDA - INTEGER.
90 * On entry, LDA specifies the first dimension of A as declared
91 * in the calling (sub) program. When SIDE = 'L' or 'l' then
92 * LDA must be at least max( 1, m ), otherwise LDA must be at
93 * least max( 1, n ).
94 * Unchanged on exit.
95 *
96 * B - COMPLEX array of DIMENSION ( LDB, n ).
97 * Before entry, the leading m by n part of the array B must
98 * contain the matrix B.
99 * Unchanged on exit.
100 *
101 * LDB - INTEGER.
102 * On entry, LDB specifies the first dimension of B as declared
103 * in the calling (sub) program. LDB must be at least
104 * max( 1, m ).
105 * Unchanged on exit.
106 *
107 * BETA - COMPLEX .
108 * On entry, BETA specifies the scalar beta. When BETA is
109 * supplied as zero then C need not be set on input.
110 * Unchanged on exit.
111 *
112 * C - COMPLEX array of DIMENSION ( LDC, n ).
113 * Before entry, the leading m by n part of the array C must
114 * contain the matrix C, except when beta is zero, in which
115 * case C need not be set on entry.
116 * On exit, the array C is overwritten by the m by n updated
117 * matrix.
118 *
119 * LDC - INTEGER.
120 * On entry, LDC specifies the first dimension of C as declared
121 * in the calling (sub) program. LDC must be at least
122 * max( 1, m ).
123 * Unchanged on exit.
124 *
125 * Further Details
126 * ===============
127 *
128 * Level 3 Blas routine.
129 *
130 * -- Written on 8-February-1989.
131 * Jack Dongarra, Argonne National Laboratory.
132 * Iain Duff, AERE Harwell.
133 * Jeremy Du Croz, Numerical Algorithms Group Ltd.
134 * Sven Hammarling, Numerical Algorithms Group Ltd.
135 *
136 * =====================================================================
137 *
138 * .. External Functions ..
139 LOGICAL LSAME
140 EXTERNAL LSAME
141 * ..
142 * .. External Subroutines ..
143 EXTERNAL XERBLA
144 * ..
145 * .. Intrinsic Functions ..
146 INTRINSIC MAX
147 * ..
148 * .. Local Scalars ..
149 COMPLEX TEMP1,TEMP2
150 INTEGER I,INFO,J,K,NROWA
151 LOGICAL UPPER
152 * ..
153 * .. Parameters ..
154 COMPLEX ONE
155 PARAMETER (ONE= (1.0E+0,0.0E+0))
156 COMPLEX ZERO
157 PARAMETER (ZERO= (0.0E+0,0.0E+0))
158 * ..
159 *
160 * Set NROWA as the number of rows of A.
161 *
162 IF (LSAME(SIDE,'L')) THEN
163 NROWA = M
164 ELSE
165 NROWA = N
166 END IF
167 UPPER = LSAME(UPLO,'U')
168 *
169 * Test the input parameters.
170 *
171 INFO = 0
172 IF ((.NOT.LSAME(SIDE,'L')) .AND. (.NOT.LSAME(SIDE,'R'))) THEN
173 INFO = 1
174 ELSE IF ((.NOT.UPPER) .AND. (.NOT.LSAME(UPLO,'L'))) THEN
175 INFO = 2
176 ELSE IF (M.LT.0) THEN
177 INFO = 3
178 ELSE IF (N.LT.0) THEN
179 INFO = 4
180 ELSE IF (LDA.LT.MAX(1,NROWA)) THEN
181 INFO = 7
182 ELSE IF (LDB.LT.MAX(1,M)) THEN
183 INFO = 9
184 ELSE IF (LDC.LT.MAX(1,M)) THEN
185 INFO = 12
186 END IF
187 IF (INFO.NE.0) THEN
188 CALL XERBLA('CSYMM ',INFO)
189 RETURN
190 END IF
191 *
192 * Quick return if possible.
193 *
194 IF ((M.EQ.0) .OR. (N.EQ.0) .OR.
195 + ((ALPHA.EQ.ZERO).AND. (BETA.EQ.ONE))) RETURN
196 *
197 * And when alpha.eq.zero.
198 *
199 IF (ALPHA.EQ.ZERO) THEN
200 IF (BETA.EQ.ZERO) THEN
201 DO 20 J = 1,N
202 DO 10 I = 1,M
203 C(I,J) = ZERO
204 10 CONTINUE
205 20 CONTINUE
206 ELSE
207 DO 40 J = 1,N
208 DO 30 I = 1,M
209 C(I,J) = BETA*C(I,J)
210 30 CONTINUE
211 40 CONTINUE
212 END IF
213 RETURN
214 END IF
215 *
216 * Start the operations.
217 *
218 IF (LSAME(SIDE,'L')) THEN
219 *
220 * Form C := alpha*A*B + beta*C.
221 *
222 IF (UPPER) THEN
223 DO 70 J = 1,N
224 DO 60 I = 1,M
225 TEMP1 = ALPHA*B(I,J)
226 TEMP2 = ZERO
227 DO 50 K = 1,I - 1
228 C(K,J) = C(K,J) + TEMP1*A(K,I)
229 TEMP2 = TEMP2 + B(K,J)*A(K,I)
230 50 CONTINUE
231 IF (BETA.EQ.ZERO) THEN
232 C(I,J) = TEMP1*A(I,I) + ALPHA*TEMP2
233 ELSE
234 C(I,J) = BETA*C(I,J) + TEMP1*A(I,I) +
235 + ALPHA*TEMP2
236 END IF
237 60 CONTINUE
238 70 CONTINUE
239 ELSE
240 DO 100 J = 1,N
241 DO 90 I = M,1,-1
242 TEMP1 = ALPHA*B(I,J)
243 TEMP2 = ZERO
244 DO 80 K = I + 1,M
245 C(K,J) = C(K,J) + TEMP1*A(K,I)
246 TEMP2 = TEMP2 + B(K,J)*A(K,I)
247 80 CONTINUE
248 IF (BETA.EQ.ZERO) THEN
249 C(I,J) = TEMP1*A(I,I) + ALPHA*TEMP2
250 ELSE
251 C(I,J) = BETA*C(I,J) + TEMP1*A(I,I) +
252 + ALPHA*TEMP2
253 END IF
254 90 CONTINUE
255 100 CONTINUE
256 END IF
257 ELSE
258 *
259 * Form C := alpha*B*A + beta*C.
260 *
261 DO 170 J = 1,N
262 TEMP1 = ALPHA*A(J,J)
263 IF (BETA.EQ.ZERO) THEN
264 DO 110 I = 1,M
265 C(I,J) = TEMP1*B(I,J)
266 110 CONTINUE
267 ELSE
268 DO 120 I = 1,M
269 C(I,J) = BETA*C(I,J) + TEMP1*B(I,J)
270 120 CONTINUE
271 END IF
272 DO 140 K = 1,J - 1
273 IF (UPPER) THEN
274 TEMP1 = ALPHA*A(K,J)
275 ELSE
276 TEMP1 = ALPHA*A(J,K)
277 END IF
278 DO 130 I = 1,M
279 C(I,J) = C(I,J) + TEMP1*B(I,K)
280 130 CONTINUE
281 140 CONTINUE
282 DO 160 K = J + 1,N
283 IF (UPPER) THEN
284 TEMP1 = ALPHA*A(J,K)
285 ELSE
286 TEMP1 = ALPHA*A(K,J)
287 END IF
288 DO 150 I = 1,M
289 C(I,J) = C(I,J) + TEMP1*B(I,K)
290 150 CONTINUE
291 160 CONTINUE
292 170 CONTINUE
293 END IF
294 *
295 RETURN
296 *
297 * End of CSYMM .
298 *
299 END
2 * .. Scalar Arguments ..
3 COMPLEX ALPHA,BETA
4 INTEGER LDA,LDB,LDC,M,N
5 CHARACTER SIDE,UPLO
6 * ..
7 * .. Array Arguments ..
8 COMPLEX A(LDA,*),B(LDB,*),C(LDC,*)
9 * ..
10 *
11 * Purpose
12 * =======
13 *
14 * CSYMM performs one of the matrix-matrix operations
15 *
16 * C := alpha*A*B + beta*C,
17 *
18 * or
19 *
20 * C := alpha*B*A + beta*C,
21 *
22 * where alpha and beta are scalars, A is a symmetric matrix and B and
23 * C are m by n matrices.
24 *
25 * Arguments
26 * ==========
27 *
28 * SIDE - CHARACTER*1.
29 * On entry, SIDE specifies whether the symmetric matrix A
30 * appears on the left or right in the operation as follows:
31 *
32 * SIDE = 'L' or 'l' C := alpha*A*B + beta*C,
33 *
34 * SIDE = 'R' or 'r' C := alpha*B*A + beta*C,
35 *
36 * Unchanged on exit.
37 *
38 * UPLO - CHARACTER*1.
39 * On entry, UPLO specifies whether the upper or lower
40 * triangular part of the symmetric matrix A is to be
41 * referenced as follows:
42 *
43 * UPLO = 'U' or 'u' Only the upper triangular part of the
44 * symmetric matrix is to be referenced.
45 *
46 * UPLO = 'L' or 'l' Only the lower triangular part of the
47 * symmetric matrix is to be referenced.
48 *
49 * Unchanged on exit.
50 *
51 * M - INTEGER.
52 * On entry, M specifies the number of rows of the matrix C.
53 * M must be at least zero.
54 * Unchanged on exit.
55 *
56 * N - INTEGER.
57 * On entry, N specifies the number of columns of the matrix C.
58 * N must be at least zero.
59 * Unchanged on exit.
60 *
61 * ALPHA - COMPLEX .
62 * On entry, ALPHA specifies the scalar alpha.
63 * Unchanged on exit.
64 *
65 * A - COMPLEX array of DIMENSION ( LDA, ka ), where ka is
66 * m when SIDE = 'L' or 'l' and is n otherwise.
67 * Before entry with SIDE = 'L' or 'l', the m by m part of
68 * the array A must contain the symmetric matrix, such that
69 * when UPLO = 'U' or 'u', the leading m by m upper triangular
70 * part of the array A must contain the upper triangular part
71 * of the symmetric matrix and the strictly lower triangular
72 * part of A is not referenced, and when UPLO = 'L' or 'l',
73 * the leading m by m lower triangular part of the array A
74 * must contain the lower triangular part of the symmetric
75 * matrix and the strictly upper triangular part of A is not
76 * referenced.
77 * Before entry with SIDE = 'R' or 'r', the n by n part of
78 * the array A must contain the symmetric matrix, such that
79 * when UPLO = 'U' or 'u', the leading n by n upper triangular
80 * part of the array A must contain the upper triangular part
81 * of the symmetric matrix and the strictly lower triangular
82 * part of A is not referenced, and when UPLO = 'L' or 'l',
83 * the leading n by n lower triangular part of the array A
84 * must contain the lower triangular part of the symmetric
85 * matrix and the strictly upper triangular part of A is not
86 * referenced.
87 * Unchanged on exit.
88 *
89 * LDA - INTEGER.
90 * On entry, LDA specifies the first dimension of A as declared
91 * in the calling (sub) program. When SIDE = 'L' or 'l' then
92 * LDA must be at least max( 1, m ), otherwise LDA must be at
93 * least max( 1, n ).
94 * Unchanged on exit.
95 *
96 * B - COMPLEX array of DIMENSION ( LDB, n ).
97 * Before entry, the leading m by n part of the array B must
98 * contain the matrix B.
99 * Unchanged on exit.
100 *
101 * LDB - INTEGER.
102 * On entry, LDB specifies the first dimension of B as declared
103 * in the calling (sub) program. LDB must be at least
104 * max( 1, m ).
105 * Unchanged on exit.
106 *
107 * BETA - COMPLEX .
108 * On entry, BETA specifies the scalar beta. When BETA is
109 * supplied as zero then C need not be set on input.
110 * Unchanged on exit.
111 *
112 * C - COMPLEX array of DIMENSION ( LDC, n ).
113 * Before entry, the leading m by n part of the array C must
114 * contain the matrix C, except when beta is zero, in which
115 * case C need not be set on entry.
116 * On exit, the array C is overwritten by the m by n updated
117 * matrix.
118 *
119 * LDC - INTEGER.
120 * On entry, LDC specifies the first dimension of C as declared
121 * in the calling (sub) program. LDC must be at least
122 * max( 1, m ).
123 * Unchanged on exit.
124 *
125 * Further Details
126 * ===============
127 *
128 * Level 3 Blas routine.
129 *
130 * -- Written on 8-February-1989.
131 * Jack Dongarra, Argonne National Laboratory.
132 * Iain Duff, AERE Harwell.
133 * Jeremy Du Croz, Numerical Algorithms Group Ltd.
134 * Sven Hammarling, Numerical Algorithms Group Ltd.
135 *
136 * =====================================================================
137 *
138 * .. External Functions ..
139 LOGICAL LSAME
140 EXTERNAL LSAME
141 * ..
142 * .. External Subroutines ..
143 EXTERNAL XERBLA
144 * ..
145 * .. Intrinsic Functions ..
146 INTRINSIC MAX
147 * ..
148 * .. Local Scalars ..
149 COMPLEX TEMP1,TEMP2
150 INTEGER I,INFO,J,K,NROWA
151 LOGICAL UPPER
152 * ..
153 * .. Parameters ..
154 COMPLEX ONE
155 PARAMETER (ONE= (1.0E+0,0.0E+0))
156 COMPLEX ZERO
157 PARAMETER (ZERO= (0.0E+0,0.0E+0))
158 * ..
159 *
160 * Set NROWA as the number of rows of A.
161 *
162 IF (LSAME(SIDE,'L')) THEN
163 NROWA = M
164 ELSE
165 NROWA = N
166 END IF
167 UPPER = LSAME(UPLO,'U')
168 *
169 * Test the input parameters.
170 *
171 INFO = 0
172 IF ((.NOT.LSAME(SIDE,'L')) .AND. (.NOT.LSAME(SIDE,'R'))) THEN
173 INFO = 1
174 ELSE IF ((.NOT.UPPER) .AND. (.NOT.LSAME(UPLO,'L'))) THEN
175 INFO = 2
176 ELSE IF (M.LT.0) THEN
177 INFO = 3
178 ELSE IF (N.LT.0) THEN
179 INFO = 4
180 ELSE IF (LDA.LT.MAX(1,NROWA)) THEN
181 INFO = 7
182 ELSE IF (LDB.LT.MAX(1,M)) THEN
183 INFO = 9
184 ELSE IF (LDC.LT.MAX(1,M)) THEN
185 INFO = 12
186 END IF
187 IF (INFO.NE.0) THEN
188 CALL XERBLA('CSYMM ',INFO)
189 RETURN
190 END IF
191 *
192 * Quick return if possible.
193 *
194 IF ((M.EQ.0) .OR. (N.EQ.0) .OR.
195 + ((ALPHA.EQ.ZERO).AND. (BETA.EQ.ONE))) RETURN
196 *
197 * And when alpha.eq.zero.
198 *
199 IF (ALPHA.EQ.ZERO) THEN
200 IF (BETA.EQ.ZERO) THEN
201 DO 20 J = 1,N
202 DO 10 I = 1,M
203 C(I,J) = ZERO
204 10 CONTINUE
205 20 CONTINUE
206 ELSE
207 DO 40 J = 1,N
208 DO 30 I = 1,M
209 C(I,J) = BETA*C(I,J)
210 30 CONTINUE
211 40 CONTINUE
212 END IF
213 RETURN
214 END IF
215 *
216 * Start the operations.
217 *
218 IF (LSAME(SIDE,'L')) THEN
219 *
220 * Form C := alpha*A*B + beta*C.
221 *
222 IF (UPPER) THEN
223 DO 70 J = 1,N
224 DO 60 I = 1,M
225 TEMP1 = ALPHA*B(I,J)
226 TEMP2 = ZERO
227 DO 50 K = 1,I - 1
228 C(K,J) = C(K,J) + TEMP1*A(K,I)
229 TEMP2 = TEMP2 + B(K,J)*A(K,I)
230 50 CONTINUE
231 IF (BETA.EQ.ZERO) THEN
232 C(I,J) = TEMP1*A(I,I) + ALPHA*TEMP2
233 ELSE
234 C(I,J) = BETA*C(I,J) + TEMP1*A(I,I) +
235 + ALPHA*TEMP2
236 END IF
237 60 CONTINUE
238 70 CONTINUE
239 ELSE
240 DO 100 J = 1,N
241 DO 90 I = M,1,-1
242 TEMP1 = ALPHA*B(I,J)
243 TEMP2 = ZERO
244 DO 80 K = I + 1,M
245 C(K,J) = C(K,J) + TEMP1*A(K,I)
246 TEMP2 = TEMP2 + B(K,J)*A(K,I)
247 80 CONTINUE
248 IF (BETA.EQ.ZERO) THEN
249 C(I,J) = TEMP1*A(I,I) + ALPHA*TEMP2
250 ELSE
251 C(I,J) = BETA*C(I,J) + TEMP1*A(I,I) +
252 + ALPHA*TEMP2
253 END IF
254 90 CONTINUE
255 100 CONTINUE
256 END IF
257 ELSE
258 *
259 * Form C := alpha*B*A + beta*C.
260 *
261 DO 170 J = 1,N
262 TEMP1 = ALPHA*A(J,J)
263 IF (BETA.EQ.ZERO) THEN
264 DO 110 I = 1,M
265 C(I,J) = TEMP1*B(I,J)
266 110 CONTINUE
267 ELSE
268 DO 120 I = 1,M
269 C(I,J) = BETA*C(I,J) + TEMP1*B(I,J)
270 120 CONTINUE
271 END IF
272 DO 140 K = 1,J - 1
273 IF (UPPER) THEN
274 TEMP1 = ALPHA*A(K,J)
275 ELSE
276 TEMP1 = ALPHA*A(J,K)
277 END IF
278 DO 130 I = 1,M
279 C(I,J) = C(I,J) + TEMP1*B(I,K)
280 130 CONTINUE
281 140 CONTINUE
282 DO 160 K = J + 1,N
283 IF (UPPER) THEN
284 TEMP1 = ALPHA*A(J,K)
285 ELSE
286 TEMP1 = ALPHA*A(K,J)
287 END IF
288 DO 150 I = 1,M
289 C(I,J) = C(I,J) + TEMP1*B(I,K)
290 150 CONTINUE
291 160 CONTINUE
292 170 CONTINUE
293 END IF
294 *
295 RETURN
296 *
297 * End of CSYMM .
298 *
299 END