|
2329
|
1 SUBROUTINE DSYRK ( UPLO, TRANS, N, K, ALPHA, A, LDA, |
|
|
2 $ BETA, C, LDC ) |
|
|
3 * .. Scalar Arguments .. |
|
|
4 CHARACTER*1 UPLO, TRANS |
|
|
5 INTEGER N, K, LDA, LDC |
|
|
6 DOUBLE PRECISION ALPHA, BETA |
|
|
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' + beta*C, |
|
|
17 * |
|
|
18 * or |
|
|
19 * |
|
|
20 * C := alpha*A'*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 * Parameters |
|
|
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' + beta*C. |
|
|
47 * |
|
|
48 * TRANS = 'T' or 't' C := alpha*A'*A + beta*C. |
|
|
49 * |
|
|
50 * TRANS = 'C' or 'c' C := alpha*A'*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 * |
|
|
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 * .. External Functions .. |
|
|
120 LOGICAL LSAME |
|
|
121 EXTERNAL LSAME |
|
|
122 * .. External Subroutines .. |
|
|
123 EXTERNAL XERBLA |
|
|
124 * .. Intrinsic Functions .. |
|
|
125 INTRINSIC MAX |
|
|
126 * .. Local Scalars .. |
|
|
127 LOGICAL UPPER |
|
|
128 INTEGER I, INFO, J, L, NROWA |
|
|
129 DOUBLE PRECISION TEMP |
|
|
130 * .. Parameters .. |
|
|
131 DOUBLE PRECISION ONE , ZERO |
|
|
132 PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) |
|
|
133 * .. |
|
|
134 * .. Executable Statements .. |
|
|
135 * |
|
|
136 * Test the input parameters. |
|
|
137 * |
|
|
138 IF( LSAME( TRANS, 'N' ) )THEN |
|
|
139 NROWA = N |
|
|
140 ELSE |
|
|
141 NROWA = K |
|
|
142 END IF |
|
|
143 UPPER = LSAME( UPLO, 'U' ) |
|
|
144 * |
|
|
145 INFO = 0 |
|
|
146 IF( ( .NOT.UPPER ).AND. |
|
|
147 $ ( .NOT.LSAME( UPLO , 'L' ) ) )THEN |
|
|
148 INFO = 1 |
|
|
149 ELSE IF( ( .NOT.LSAME( TRANS, 'N' ) ).AND. |
|
|
150 $ ( .NOT.LSAME( TRANS, 'T' ) ).AND. |
|
|
151 $ ( .NOT.LSAME( TRANS, 'C' ) ) )THEN |
|
|
152 INFO = 2 |
|
|
153 ELSE IF( N .LT.0 )THEN |
|
|
154 INFO = 3 |
|
|
155 ELSE IF( K .LT.0 )THEN |
|
|
156 INFO = 4 |
|
|
157 ELSE IF( LDA.LT.MAX( 1, NROWA ) )THEN |
|
|
158 INFO = 7 |
|
|
159 ELSE IF( LDC.LT.MAX( 1, N ) )THEN |
|
|
160 INFO = 10 |
|
|
161 END IF |
|
|
162 IF( INFO.NE.0 )THEN |
|
|
163 CALL XERBLA( 'DSYRK ', INFO ) |
|
|
164 RETURN |
|
|
165 END IF |
|
|
166 * |
|
|
167 * Quick return if possible. |
|
|
168 * |
|
|
169 IF( ( N.EQ.0 ).OR. |
|
|
170 $ ( ( ( ALPHA.EQ.ZERO ).OR.( K.EQ.0 ) ).AND.( BETA.EQ.ONE ) ) ) |
|
|
171 $ RETURN |
|
|
172 * |
|
|
173 * And when alpha.eq.zero. |
|
|
174 * |
|
|
175 IF( ALPHA.EQ.ZERO )THEN |
|
|
176 IF( UPPER )THEN |
|
|
177 IF( BETA.EQ.ZERO )THEN |
|
|
178 DO 20, J = 1, N |
|
|
179 DO 10, I = 1, J |
|
|
180 C( I, J ) = ZERO |
|
|
181 10 CONTINUE |
|
|
182 20 CONTINUE |
|
|
183 ELSE |
|
|
184 DO 40, J = 1, N |
|
|
185 DO 30, I = 1, J |
|
|
186 C( I, J ) = BETA*C( I, J ) |
|
|
187 30 CONTINUE |
|
|
188 40 CONTINUE |
|
|
189 END IF |
|
|
190 ELSE |
|
|
191 IF( BETA.EQ.ZERO )THEN |
|
|
192 DO 60, J = 1, N |
|
|
193 DO 50, I = J, N |
|
|
194 C( I, J ) = ZERO |
|
|
195 50 CONTINUE |
|
|
196 60 CONTINUE |
|
|
197 ELSE |
|
|
198 DO 80, J = 1, N |
|
|
199 DO 70, I = J, N |
|
|
200 C( I, J ) = BETA*C( I, J ) |
|
|
201 70 CONTINUE |
|
|
202 80 CONTINUE |
|
|
203 END IF |
|
|
204 END IF |
|
|
205 RETURN |
|
|
206 END IF |
|
|
207 * |
|
|
208 * Start the operations. |
|
|
209 * |
|
|
210 IF( LSAME( TRANS, 'N' ) )THEN |
|
|
211 * |
|
|
212 * Form C := alpha*A*A' + beta*C. |
|
|
213 * |
|
|
214 IF( UPPER )THEN |
|
|
215 DO 130, J = 1, N |
|
|
216 IF( BETA.EQ.ZERO )THEN |
|
|
217 DO 90, I = 1, J |
|
|
218 C( I, J ) = ZERO |
|
|
219 90 CONTINUE |
|
|
220 ELSE IF( BETA.NE.ONE )THEN |
|
|
221 DO 100, I = 1, J |
|
|
222 C( I, J ) = BETA*C( I, J ) |
|
|
223 100 CONTINUE |
|
|
224 END IF |
|
|
225 DO 120, L = 1, K |
|
|
226 IF( A( J, L ).NE.ZERO )THEN |
|
|
227 TEMP = ALPHA*A( J, L ) |
|
|
228 DO 110, I = 1, J |
|
|
229 C( I, J ) = C( I, J ) + TEMP*A( I, L ) |
|
|
230 110 CONTINUE |
|
|
231 END IF |
|
|
232 120 CONTINUE |
|
|
233 130 CONTINUE |
|
|
234 ELSE |
|
|
235 DO 180, J = 1, N |
|
|
236 IF( BETA.EQ.ZERO )THEN |
|
|
237 DO 140, I = J, N |
|
|
238 C( I, J ) = ZERO |
|
|
239 140 CONTINUE |
|
|
240 ELSE IF( BETA.NE.ONE )THEN |
|
|
241 DO 150, I = J, N |
|
|
242 C( I, J ) = BETA*C( I, J ) |
|
|
243 150 CONTINUE |
|
|
244 END IF |
|
|
245 DO 170, L = 1, K |
|
|
246 IF( A( J, L ).NE.ZERO )THEN |
|
|
247 TEMP = ALPHA*A( J, L ) |
|
|
248 DO 160, I = J, N |
|
|
249 C( I, J ) = C( I, J ) + TEMP*A( I, L ) |
|
|
250 160 CONTINUE |
|
|
251 END IF |
|
|
252 170 CONTINUE |
|
|
253 180 CONTINUE |
|
|
254 END IF |
|
|
255 ELSE |
|
|
256 * |
|
|
257 * Form C := alpha*A'*A + beta*C. |
|
|
258 * |
|
|
259 IF( UPPER )THEN |
|
|
260 DO 210, J = 1, N |
|
|
261 DO 200, I = 1, J |
|
|
262 TEMP = ZERO |
|
|
263 DO 190, L = 1, K |
|
|
264 TEMP = TEMP + A( L, I )*A( L, J ) |
|
|
265 190 CONTINUE |
|
|
266 IF( BETA.EQ.ZERO )THEN |
|
|
267 C( I, J ) = ALPHA*TEMP |
|
|
268 ELSE |
|
|
269 C( I, J ) = ALPHA*TEMP + BETA*C( I, J ) |
|
|
270 END IF |
|
|
271 200 CONTINUE |
|
|
272 210 CONTINUE |
|
|
273 ELSE |
|
|
274 DO 240, J = 1, N |
|
|
275 DO 230, I = J, N |
|
|
276 TEMP = ZERO |
|
|
277 DO 220, L = 1, K |
|
|
278 TEMP = TEMP + A( L, I )*A( L, J ) |
|
|
279 220 CONTINUE |
|
|
280 IF( BETA.EQ.ZERO )THEN |
|
|
281 C( I, J ) = ALPHA*TEMP |
|
|
282 ELSE |
|
|
283 C( I, J ) = ALPHA*TEMP + BETA*C( I, J ) |
|
|
284 END IF |
|
|
285 230 CONTINUE |
|
|
286 240 CONTINUE |
|
|
287 END IF |
|
|
288 END IF |
|
|
289 * |
|
|
290 RETURN |
|
|
291 * |
|
|
292 * End of DSYRK . |
|
|
293 * |
|
|
294 END |
