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