1 SUBROUTINE DSYMM(SIDE,UPLO,M,N,ALPHA,A,LDA,B,LDB,BETA,C,LDC)
2 * .. Scalar Arguments ..
3 DOUBLE PRECISION ALPHA,BETA
4 INTEGER LDA,LDB,LDC,M,N
5 CHARACTER SIDE,UPLO
6 * ..
7 * .. Array Arguments ..
8 DOUBLE PRECISION A(LDA,*),B(LDB,*),C(LDC,*)
9 * ..
10 *
11 * Purpose
12 * =======
13 *
14 * DSYMM 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 - DOUBLE PRECISION.
62 * On entry, ALPHA specifies the scalar alpha.
63 * Unchanged on exit.
64 *
65 * A - DOUBLE PRECISION 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 - DOUBLE PRECISION 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 - DOUBLE PRECISION.
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 - DOUBLE PRECISION 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 DOUBLE PRECISION TEMP1,TEMP2
150 INTEGER I,INFO,J,K,NROWA
151 LOGICAL UPPER
152 * ..
153 * .. Parameters ..
154 DOUBLE PRECISION ONE,ZERO
155 PARAMETER (ONE=1.0D+0,ZERO=0.0D+0)
156 * ..
157 *
158 * Set NROWA as the number of rows of A.
159 *
160 IF (LSAME(SIDE,'L')) THEN
161 NROWA = M
162 ELSE
163 NROWA = N
164 END IF
165 UPPER = LSAME(UPLO,'U')
166 *
167 * Test the input parameters.
168 *
169 INFO = 0
170 IF ((.NOT.LSAME(SIDE,'L')) .AND. (.NOT.LSAME(SIDE,'R'))) THEN
171 INFO = 1
172 ELSE IF ((.NOT.UPPER) .AND. (.NOT.LSAME(UPLO,'L'))) THEN
173 INFO = 2
174 ELSE IF (M.LT.0) THEN
175 INFO = 3
176 ELSE IF (N.LT.0) THEN
177 INFO = 4
178 ELSE IF (LDA.LT.MAX(1,NROWA)) THEN
179 INFO = 7
180 ELSE IF (LDB.LT.MAX(1,M)) THEN
181 INFO = 9
182 ELSE IF (LDC.LT.MAX(1,M)) THEN
183 INFO = 12
184 END IF
185 IF (INFO.NE.0) THEN
186 CALL XERBLA('DSYMM ',INFO)
187 RETURN
188 END IF
189 *
190 * Quick return if possible.
191 *
192 IF ((M.EQ.0) .OR. (N.EQ.0) .OR.
193 + ((ALPHA.EQ.ZERO).AND. (BETA.EQ.ONE))) RETURN
194 *
195 * And when alpha.eq.zero.
196 *
197 IF (ALPHA.EQ.ZERO) THEN
198 IF (BETA.EQ.ZERO) THEN
199 DO 20 J = 1,N
200 DO 10 I = 1,M
201 C(I,J) = ZERO
202 10 CONTINUE
203 20 CONTINUE
204 ELSE
205 DO 40 J = 1,N
206 DO 30 I = 1,M
207 C(I,J) = BETA*C(I,J)
208 30 CONTINUE
209 40 CONTINUE
210 END IF
211 RETURN
212 END IF
213 *
214 * Start the operations.
215 *
216 IF (LSAME(SIDE,'L')) THEN
217 *
218 * Form C := alpha*A*B + beta*C.
219 *
220 IF (UPPER) THEN
221 DO 70 J = 1,N
222 DO 60 I = 1,M
223 TEMP1 = ALPHA*B(I,J)
224 TEMP2 = ZERO
225 DO 50 K = 1,I - 1
226 C(K,J) = C(K,J) + TEMP1*A(K,I)
227 TEMP2 = TEMP2 + B(K,J)*A(K,I)
228 50 CONTINUE
229 IF (BETA.EQ.ZERO) THEN
230 C(I,J) = TEMP1*A(I,I) + ALPHA*TEMP2
231 ELSE
232 C(I,J) = BETA*C(I,J) + TEMP1*A(I,I) +
233 + ALPHA*TEMP2
234 END IF
235 60 CONTINUE
236 70 CONTINUE
237 ELSE
238 DO 100 J = 1,N
239 DO 90 I = M,1,-1
240 TEMP1 = ALPHA*B(I,J)
241 TEMP2 = ZERO
242 DO 80 K = I + 1,M
243 C(K,J) = C(K,J) + TEMP1*A(K,I)
244 TEMP2 = TEMP2 + B(K,J)*A(K,I)
245 80 CONTINUE
246 IF (BETA.EQ.ZERO) THEN
247 C(I,J) = TEMP1*A(I,I) + ALPHA*TEMP2
248 ELSE
249 C(I,J) = BETA*C(I,J) + TEMP1*A(I,I) +
250 + ALPHA*TEMP2
251 END IF
252 90 CONTINUE
253 100 CONTINUE
254 END IF
255 ELSE
256 *
257 * Form C := alpha*B*A + beta*C.
258 *
259 DO 170 J = 1,N
260 TEMP1 = ALPHA*A(J,J)
261 IF (BETA.EQ.ZERO) THEN
262 DO 110 I = 1,M
263 C(I,J) = TEMP1*B(I,J)
264 110 CONTINUE
265 ELSE
266 DO 120 I = 1,M
267 C(I,J) = BETA*C(I,J) + TEMP1*B(I,J)
268 120 CONTINUE
269 END IF
270 DO 140 K = 1,J - 1
271 IF (UPPER) THEN
272 TEMP1 = ALPHA*A(K,J)
273 ELSE
274 TEMP1 = ALPHA*A(J,K)
275 END IF
276 DO 130 I = 1,M
277 C(I,J) = C(I,J) + TEMP1*B(I,K)
278 130 CONTINUE
279 140 CONTINUE
280 DO 160 K = J + 1,N
281 IF (UPPER) THEN
282 TEMP1 = ALPHA*A(J,K)
283 ELSE
284 TEMP1 = ALPHA*A(K,J)
285 END IF
286 DO 150 I = 1,M
287 C(I,J) = C(I,J) + TEMP1*B(I,K)
288 150 CONTINUE
289 160 CONTINUE
290 170 CONTINUE
291 END IF
292 *
293 RETURN
294 *
295 * End of DSYMM .
296 *
297 END
2 * .. Scalar Arguments ..
3 DOUBLE PRECISION ALPHA,BETA
4 INTEGER LDA,LDB,LDC,M,N
5 CHARACTER SIDE,UPLO
6 * ..
7 * .. Array Arguments ..
8 DOUBLE PRECISION A(LDA,*),B(LDB,*),C(LDC,*)
9 * ..
10 *
11 * Purpose
12 * =======
13 *
14 * DSYMM 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 - DOUBLE PRECISION.
62 * On entry, ALPHA specifies the scalar alpha.
63 * Unchanged on exit.
64 *
65 * A - DOUBLE PRECISION 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 - DOUBLE PRECISION 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 - DOUBLE PRECISION.
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 - DOUBLE PRECISION 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 DOUBLE PRECISION TEMP1,TEMP2
150 INTEGER I,INFO,J,K,NROWA
151 LOGICAL UPPER
152 * ..
153 * .. Parameters ..
154 DOUBLE PRECISION ONE,ZERO
155 PARAMETER (ONE=1.0D+0,ZERO=0.0D+0)
156 * ..
157 *
158 * Set NROWA as the number of rows of A.
159 *
160 IF (LSAME(SIDE,'L')) THEN
161 NROWA = M
162 ELSE
163 NROWA = N
164 END IF
165 UPPER = LSAME(UPLO,'U')
166 *
167 * Test the input parameters.
168 *
169 INFO = 0
170 IF ((.NOT.LSAME(SIDE,'L')) .AND. (.NOT.LSAME(SIDE,'R'))) THEN
171 INFO = 1
172 ELSE IF ((.NOT.UPPER) .AND. (.NOT.LSAME(UPLO,'L'))) THEN
173 INFO = 2
174 ELSE IF (M.LT.0) THEN
175 INFO = 3
176 ELSE IF (N.LT.0) THEN
177 INFO = 4
178 ELSE IF (LDA.LT.MAX(1,NROWA)) THEN
179 INFO = 7
180 ELSE IF (LDB.LT.MAX(1,M)) THEN
181 INFO = 9
182 ELSE IF (LDC.LT.MAX(1,M)) THEN
183 INFO = 12
184 END IF
185 IF (INFO.NE.0) THEN
186 CALL XERBLA('DSYMM ',INFO)
187 RETURN
188 END IF
189 *
190 * Quick return if possible.
191 *
192 IF ((M.EQ.0) .OR. (N.EQ.0) .OR.
193 + ((ALPHA.EQ.ZERO).AND. (BETA.EQ.ONE))) RETURN
194 *
195 * And when alpha.eq.zero.
196 *
197 IF (ALPHA.EQ.ZERO) THEN
198 IF (BETA.EQ.ZERO) THEN
199 DO 20 J = 1,N
200 DO 10 I = 1,M
201 C(I,J) = ZERO
202 10 CONTINUE
203 20 CONTINUE
204 ELSE
205 DO 40 J = 1,N
206 DO 30 I = 1,M
207 C(I,J) = BETA*C(I,J)
208 30 CONTINUE
209 40 CONTINUE
210 END IF
211 RETURN
212 END IF
213 *
214 * Start the operations.
215 *
216 IF (LSAME(SIDE,'L')) THEN
217 *
218 * Form C := alpha*A*B + beta*C.
219 *
220 IF (UPPER) THEN
221 DO 70 J = 1,N
222 DO 60 I = 1,M
223 TEMP1 = ALPHA*B(I,J)
224 TEMP2 = ZERO
225 DO 50 K = 1,I - 1
226 C(K,J) = C(K,J) + TEMP1*A(K,I)
227 TEMP2 = TEMP2 + B(K,J)*A(K,I)
228 50 CONTINUE
229 IF (BETA.EQ.ZERO) THEN
230 C(I,J) = TEMP1*A(I,I) + ALPHA*TEMP2
231 ELSE
232 C(I,J) = BETA*C(I,J) + TEMP1*A(I,I) +
233 + ALPHA*TEMP2
234 END IF
235 60 CONTINUE
236 70 CONTINUE
237 ELSE
238 DO 100 J = 1,N
239 DO 90 I = M,1,-1
240 TEMP1 = ALPHA*B(I,J)
241 TEMP2 = ZERO
242 DO 80 K = I + 1,M
243 C(K,J) = C(K,J) + TEMP1*A(K,I)
244 TEMP2 = TEMP2 + B(K,J)*A(K,I)
245 80 CONTINUE
246 IF (BETA.EQ.ZERO) THEN
247 C(I,J) = TEMP1*A(I,I) + ALPHA*TEMP2
248 ELSE
249 C(I,J) = BETA*C(I,J) + TEMP1*A(I,I) +
250 + ALPHA*TEMP2
251 END IF
252 90 CONTINUE
253 100 CONTINUE
254 END IF
255 ELSE
256 *
257 * Form C := alpha*B*A + beta*C.
258 *
259 DO 170 J = 1,N
260 TEMP1 = ALPHA*A(J,J)
261 IF (BETA.EQ.ZERO) THEN
262 DO 110 I = 1,M
263 C(I,J) = TEMP1*B(I,J)
264 110 CONTINUE
265 ELSE
266 DO 120 I = 1,M
267 C(I,J) = BETA*C(I,J) + TEMP1*B(I,J)
268 120 CONTINUE
269 END IF
270 DO 140 K = 1,J - 1
271 IF (UPPER) THEN
272 TEMP1 = ALPHA*A(K,J)
273 ELSE
274 TEMP1 = ALPHA*A(J,K)
275 END IF
276 DO 130 I = 1,M
277 C(I,J) = C(I,J) + TEMP1*B(I,K)
278 130 CONTINUE
279 140 CONTINUE
280 DO 160 K = J + 1,N
281 IF (UPPER) THEN
282 TEMP1 = ALPHA*A(J,K)
283 ELSE
284 TEMP1 = ALPHA*A(K,J)
285 END IF
286 DO 150 I = 1,M
287 C(I,J) = C(I,J) + TEMP1*B(I,K)
288 150 CONTINUE
289 160 CONTINUE
290 170 CONTINUE
291 END IF
292 *
293 RETURN
294 *
295 * End of DSYMM .
296 *
297 END