Book a Demo!
CoCalc Logo Icon
StoreFeaturesDocsShareSupportNewsAboutPoliciesSign UpSign In
ElmerCSC
GitHub Repository: ElmerCSC/elmerfem
Path: blob/devel/mathlibs/src/blas/csyrk.f
5227 views
1
SUBROUTINE CSYRK ( 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
COMPLEX ALPHA, BETA
7
* .. Array Arguments ..
8
COMPLEX A( LDA, * ), C( LDC, * )
9
* ..
10
*
11
* Purpose
12
* =======
13
*
14
* CSYRK 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
* 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 .
65
* On entry, ALPHA specifies the scalar alpha.
66
* Unchanged on exit.
67
*
68
* A - COMPLEX 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 .
84
* On entry, BETA specifies the scalar beta.
85
* Unchanged on exit.
86
*
87
* C - COMPLEX 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
*
108
* Level 3 Blas routine.
109
*
110
* -- Written on 8-February-1989.
111
* Jack Dongarra, Argonne National Laboratory.
112
* Iain Duff, AERE Harwell.
113
* Jeremy Du Croz, Numerical Algorithms Group Ltd.
114
* Sven Hammarling, Numerical Algorithms Group Ltd.
115
*
116
*
117
* .. External Functions ..
118
LOGICAL LSAME
119
EXTERNAL LSAME
120
* .. External Subroutines ..
121
EXTERNAL XERBLA
122
* .. Intrinsic Functions ..
123
INTRINSIC MAX
124
* .. Local Scalars ..
125
LOGICAL UPPER
126
INTEGER I, INFO, J, L, NROWA
127
COMPLEX TEMP
128
* .. Parameters ..
129
COMPLEX ONE
130
PARAMETER ( ONE = ( 1.0E+0, 0.0E+0 ) )
131
COMPLEX ZERO
132
PARAMETER ( ZERO = ( 0.0E+0, 0.0E+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' ) ) )THEN
151
INFO = 2
152
ELSE IF( N .LT.0 )THEN
153
INFO = 3
154
ELSE IF( K .LT.0 )THEN
155
INFO = 4
156
ELSE IF( LDA.LT.MAX( 1, NROWA ) )THEN
157
INFO = 7
158
ELSE IF( LDC.LT.MAX( 1, N ) )THEN
159
INFO = 10
160
END IF
161
IF( INFO.NE.0 )THEN
162
CALL XERBLA( 'CSYRK ', INFO )
163
RETURN
164
END IF
165
*
166
* Quick return if possible.
167
*
168
IF( ( N.EQ.0 ).OR.
169
$ ( ( ( ALPHA.EQ.ZERO ).OR.( K.EQ.0 ) ).AND.( BETA.EQ.ONE ) ) )
170
$ RETURN
171
*
172
* And when alpha.eq.zero.
173
*
174
IF( ALPHA.EQ.ZERO )THEN
175
IF( UPPER )THEN
176
IF( BETA.EQ.ZERO )THEN
177
DO 20, J = 1, N
178
DO 10, I = 1, J
179
C( I, J ) = ZERO
180
10 CONTINUE
181
20 CONTINUE
182
ELSE
183
DO 40, J = 1, N
184
DO 30, I = 1, J
185
C( I, J ) = BETA*C( I, J )
186
30 CONTINUE
187
40 CONTINUE
188
END IF
189
ELSE
190
IF( BETA.EQ.ZERO )THEN
191
DO 60, J = 1, N
192
DO 50, I = J, N
193
C( I, J ) = ZERO
194
50 CONTINUE
195
60 CONTINUE
196
ELSE
197
DO 80, J = 1, N
198
DO 70, I = J, N
199
C( I, J ) = BETA*C( I, J )
200
70 CONTINUE
201
80 CONTINUE
202
END IF
203
END IF
204
RETURN
205
END IF
206
*
207
* Start the operations.
208
*
209
IF( LSAME( TRANS, 'N' ) )THEN
210
*
211
* Form C := alpha*A*A' + beta*C.
212
*
213
IF( UPPER )THEN
214
DO 130, J = 1, N
215
IF( BETA.EQ.ZERO )THEN
216
DO 90, I = 1, J
217
C( I, J ) = ZERO
218
90 CONTINUE
219
ELSE IF( BETA.NE.ONE )THEN
220
DO 100, I = 1, J
221
C( I, J ) = BETA*C( I, J )
222
100 CONTINUE
223
END IF
224
DO 120, L = 1, K
225
IF( A( J, L ).NE.ZERO )THEN
226
TEMP = ALPHA*A( J, L )
227
DO 110, I = 1, J
228
C( I, J ) = C( I, J ) + TEMP*A( I, L )
229
110 CONTINUE
230
END IF
231
120 CONTINUE
232
130 CONTINUE
233
ELSE
234
DO 180, J = 1, N
235
IF( BETA.EQ.ZERO )THEN
236
DO 140, I = J, N
237
C( I, J ) = ZERO
238
140 CONTINUE
239
ELSE IF( BETA.NE.ONE )THEN
240
DO 150, I = J, N
241
C( I, J ) = BETA*C( I, J )
242
150 CONTINUE
243
END IF
244
DO 170, L = 1, K
245
IF( A( J, L ).NE.ZERO )THEN
246
TEMP = ALPHA*A( J, L )
247
DO 160, I = J, N
248
C( I, J ) = C( I, J ) + TEMP*A( I, L )
249
160 CONTINUE
250
END IF
251
170 CONTINUE
252
180 CONTINUE
253
END IF
254
ELSE
255
*
256
* Form C := alpha*A'*A + beta*C.
257
*
258
IF( UPPER )THEN
259
DO 210, J = 1, N
260
DO 200, I = 1, J
261
TEMP = ZERO
262
DO 190, L = 1, K
263
TEMP = TEMP + A( L, I )*A( L, J )
264
190 CONTINUE
265
IF( BETA.EQ.ZERO )THEN
266
C( I, J ) = ALPHA*TEMP
267
ELSE
268
C( I, J ) = ALPHA*TEMP + BETA*C( I, J )
269
END IF
270
200 CONTINUE
271
210 CONTINUE
272
ELSE
273
DO 240, J = 1, N
274
DO 230, I = J, N
275
TEMP = ZERO
276
DO 220, L = 1, K
277
TEMP = TEMP + A( L, I )*A( L, J )
278
220 CONTINUE
279
IF( BETA.EQ.ZERO )THEN
280
C( I, J ) = ALPHA*TEMP
281
ELSE
282
C( I, J ) = ALPHA*TEMP + BETA*C( I, J )
283
END IF
284
230 CONTINUE
285
240 CONTINUE
286
END IF
287
END IF
288
*
289
RETURN
290
*
291
* End of CSYRK .
292
*
293
END
294
295