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