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