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