2
************************************************************************
4
SUBROUTINE DSYR2K( UPLO, TRANS, N, K, ALPHA, A, LDA, B, LDB,
6
* .. Scalar Arguments ..
7
CHARACTER*1 UPLO, TRANS
8
INTEGER N, K, LDA, LDB, LDC
9
DOUBLE PRECISION ALPHA, BETA
10
* .. Array Arguments ..
11
DOUBLE PRECISION A( LDA, * ), B( LDB, * ), C( LDC, * )
17
* DSYR2K performs one of the symmetric rank 2k operations
19
* C := alpha*A*B' + alpha*B*A' + beta*C,
23
* C := alpha*A'*B + alpha*B'*A + beta*C,
25
* where alpha and beta are scalars, C is an n by n symmetric matrix
26
* and A and B are n by k matrices in the first case and k by n
27
* matrices in the second case.
33
* On entry, UPLO specifies whether the upper or lower
34
* triangular part of the array C is to be referenced as
37
* UPLO = 'U' or 'u' Only the upper triangular part of C
38
* is to be referenced.
40
* UPLO = 'L' or 'l' Only the lower triangular part of C
41
* is to be referenced.
45
* TRANS - CHARACTER*1.
46
* On entry, TRANS specifies the operation to be performed as
49
* TRANS = 'N' or 'n' C := alpha*A*B' + alpha*B*A' +
52
* TRANS = 'T' or 't' C := alpha*A'*B + alpha*B'*A +
55
* TRANS = 'C' or 'c' C := alpha*A'*B + alpha*B'*A +
61
* On entry, N specifies the order of the matrix C. N must be
66
* On entry with TRANS = 'N' or 'n', K specifies the number
67
* of columns of the matrices A and B, and on entry with
68
* TRANS = 'T' or 't' or 'C' or 'c', K specifies the number
69
* of rows of the matrices A and B. K must be at least zero.
72
* ALPHA - DOUBLE PRECISION.
73
* On entry, ALPHA specifies the scalar alpha.
76
* A - DOUBLE PRECISION array of DIMENSION ( LDA, ka ), where ka is
77
* k when TRANS = 'N' or 'n', and is n otherwise.
78
* Before entry with TRANS = 'N' or 'n', the leading n by k
79
* part of the array A must contain the matrix A, otherwise
80
* the leading k by n part of the array A must contain the
85
* On entry, LDA specifies the first dimension of A as declared
86
* in the calling (sub) program. When TRANS = 'N' or 'n'
87
* then LDA must be at least max( 1, n ), otherwise LDA must
88
* be at least max( 1, k ).
91
* B - DOUBLE PRECISION array of DIMENSION ( LDB, kb ), where kb is
92
* k when TRANS = 'N' or 'n', and is n otherwise.
93
* Before entry with TRANS = 'N' or 'n', the leading n by k
94
* part of the array B must contain the matrix B, otherwise
95
* the leading k by n part of the array B must contain the
100
* On entry, LDB specifies the first dimension of B as declared
101
* in the calling (sub) program. When TRANS = 'N' or 'n'
102
* then LDB must be at least max( 1, n ), otherwise LDB must
103
* be at least max( 1, k ).
106
* BETA - DOUBLE PRECISION.
107
* On entry, BETA specifies the scalar beta.
110
* C - DOUBLE PRECISION array of DIMENSION ( LDC, n ).
111
* Before entry with UPLO = 'U' or 'u', the leading n by n
112
* upper triangular part of the array C must contain the upper
113
* triangular part of the symmetric matrix and the strictly
114
* lower triangular part of C is not referenced. On exit, the
115
* upper triangular part of the array C is overwritten by the
116
* upper triangular part of the updated matrix.
117
* Before entry with UPLO = 'L' or 'l', the leading n by n
118
* lower triangular part of the array C must contain the lower
119
* triangular part of the symmetric matrix and the strictly
120
* upper triangular part of C is not referenced. On exit, the
121
* lower triangular part of the array C is overwritten by the
122
* lower triangular part of the updated matrix.
125
* On entry, LDC specifies the first dimension of C as declared
126
* in the calling (sub) program. LDC must be at least
131
* Level 3 Blas routine.
134
* -- Written on 8-February-1989.
135
* Jack Dongarra, Argonne National Laboratory.
136
* Iain Duff, AERE Harwell.
137
* Jeremy Du Croz, Numerical Algorithms Group Ltd.
138
* Sven Hammarling, Numerical Algorithms Group Ltd.
141
* .. External Functions ..
144
* .. External Subroutines ..
146
* .. Intrinsic Functions ..
148
* .. Local Scalars ..
150
INTEGER I, INFO, J, L, NROWA
151
DOUBLE PRECISION TEMP1, TEMP2
153
DOUBLE PRECISION ONE , ZERO
154
PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
156
* .. Executable Statements ..
158
* Test the input parameters.
160
IF( LSAME( TRANS, 'N' ) )THEN
165
UPPER = LSAME( UPLO, 'U' )
168
IF( ( .NOT.UPPER ).AND.
169
$ ( .NOT.LSAME( UPLO , 'L' ) ) )THEN
171
ELSE IF( ( .NOT.LSAME( TRANS, 'N' ) ).AND.
172
$ ( .NOT.LSAME( TRANS, 'T' ) ).AND.
173
$ ( .NOT.LSAME( TRANS, 'C' ) ) )THEN
175
ELSE IF( N .LT.0 )THEN
177
ELSE IF( K .LT.0 )THEN
179
ELSE IF( LDA.LT.MAX( 1, NROWA ) )THEN
181
ELSE IF( LDB.LT.MAX( 1, NROWA ) )THEN
183
ELSE IF( LDC.LT.MAX( 1, N ) )THEN
187
CALL XERBLA( 'DSYR2K', INFO )
191
* Quick return if possible.
194
$ ( ( ( ALPHA.EQ.ZERO ).OR.( K.EQ.0 ) ).AND.( BETA.EQ.ONE ) ) )
197
* And when alpha.eq.zero.
199
IF( ALPHA.EQ.ZERO )THEN
201
IF( BETA.EQ.ZERO )THEN
210
C( I, J ) = BETA*C( I, J )
215
IF( BETA.EQ.ZERO )THEN
224
C( I, J ) = BETA*C( I, J )
232
* Start the operations.
234
IF( LSAME( TRANS, 'N' ) )THEN
236
* Form C := alpha*A*B' + alpha*B*A' + C.
240
IF( BETA.EQ.ZERO )THEN
244
ELSE IF( BETA.NE.ONE )THEN
246
C( I, J ) = BETA*C( I, J )
250
IF( ( A( J, L ).NE.ZERO ).OR.
251
$ ( B( J, L ).NE.ZERO ) )THEN
252
TEMP1 = ALPHA*B( J, L )
253
TEMP2 = ALPHA*A( J, L )
255
C( I, J ) = C( I, J ) +
256
$ A( I, L )*TEMP1 + B( I, L )*TEMP2
263
IF( BETA.EQ.ZERO )THEN
267
ELSE IF( BETA.NE.ONE )THEN
269
C( I, J ) = BETA*C( I, J )
273
IF( ( A( J, L ).NE.ZERO ).OR.
274
$ ( B( J, L ).NE.ZERO ) )THEN
275
TEMP1 = ALPHA*B( J, L )
276
TEMP2 = ALPHA*A( J, L )
278
C( I, J ) = C( I, J ) +
279
$ A( I, L )*TEMP1 + B( I, L )*TEMP2
287
* Form C := alpha*A'*B + alpha*B'*A + C.
295
TEMP1 = TEMP1 + A( L, I )*B( L, J )
296
TEMP2 = TEMP2 + B( L, I )*A( L, J )
298
IF( BETA.EQ.ZERO )THEN
299
C( I, J ) = ALPHA*TEMP1 + ALPHA*TEMP2
301
C( I, J ) = BETA *C( I, J ) +
302
$ ALPHA*TEMP1 + ALPHA*TEMP2
312
TEMP1 = TEMP1 + A( L, I )*B( L, J )
313
TEMP2 = TEMP2 + B( L, I )*A( L, J )
315
IF( BETA.EQ.ZERO )THEN
316
C( I, J ) = ALPHA*TEMP1 + ALPHA*TEMP2
318
C( I, J ) = BETA *C( I, J ) +
319
$ ALPHA*TEMP1 + ALPHA*TEMP2
added
removed
libcruft/blas/dsyr2k.f
