1
SUBROUTINE DSYRK ( UPLO, TRANS, N, K, ALPHA, A, LDA,
3
* .. Scalar Arguments ..
4
CHARACTER*1 UPLO, TRANS
6
DOUBLE PRECISION ALPHA, BETA
7
* .. Array Arguments ..
8
DOUBLE PRECISION A( LDA, * ), C( LDC, * )
14
* DSYRK performs one of the symmetric rank k operations
16
* C := alpha*A*A' + beta*C,
20
* C := alpha*A'*A + beta*C,
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
30
* On entry, UPLO specifies whether the upper or lower
31
* triangular part of the array C is to be referenced as
34
* UPLO = 'U' or 'u' Only the upper triangular part of C
35
* is to be referenced.
37
* UPLO = 'L' or 'l' Only the lower triangular part of C
38
* is to be referenced.
42
* TRANS - CHARACTER*1.
43
* On entry, TRANS specifies the operation to be performed as
46
* TRANS = 'N' or 'n' C := alpha*A*A' + beta*C.
48
* TRANS = 'T' or 't' C := alpha*A'*A + beta*C.
50
* TRANS = 'C' or 'c' C := alpha*A'*A + beta*C.
55
* On entry, N specifies the order of the matrix C. N must be
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.
66
* ALPHA - DOUBLE PRECISION.
67
* On entry, ALPHA specifies the scalar alpha.
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
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 ).
85
* BETA - DOUBLE PRECISION.
86
* On entry, BETA specifies the scalar beta.
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.
104
* On entry, LDC specifies the first dimension of C as declared
105
* in the calling (sub) program. LDC must be at least
110
* Level 3 Blas routine.
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.
119
* .. External Functions ..
122
* .. External Subroutines ..
124
* .. Intrinsic Functions ..
126
* .. Local Scalars ..
128
INTEGER I, INFO, J, L, NROWA
129
DOUBLE PRECISION TEMP
131
DOUBLE PRECISION ONE , ZERO
132
PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
134
* .. Executable Statements ..
136
* Test the input parameters.
138
IF( LSAME( TRANS, 'N' ) )THEN
143
UPPER = LSAME( UPLO, 'U' )
146
IF( ( .NOT.UPPER ).AND.
147
$ ( .NOT.LSAME( UPLO , 'L' ) ) )THEN
149
ELSE IF( ( .NOT.LSAME( TRANS, 'N' ) ).AND.
150
$ ( .NOT.LSAME( TRANS, 'T' ) ).AND.
151
$ ( .NOT.LSAME( TRANS, 'C' ) ) )THEN
153
ELSE IF( N .LT.0 )THEN
155
ELSE IF( K .LT.0 )THEN
157
ELSE IF( LDA.LT.MAX( 1, NROWA ) )THEN
159
ELSE IF( LDC.LT.MAX( 1, N ) )THEN
163
CALL XERBLA( 'DSYRK ', INFO )
167
* Quick return if possible.
170
$ ( ( ( ALPHA.EQ.ZERO ).OR.( K.EQ.0 ) ).AND.( BETA.EQ.ONE ) ) )
173
* And when alpha.eq.zero.
175
IF( ALPHA.EQ.ZERO )THEN
177
IF( BETA.EQ.ZERO )THEN
186
C( I, J ) = BETA*C( I, J )
191
IF( BETA.EQ.ZERO )THEN
200
C( I, J ) = BETA*C( I, J )
208
* Start the operations.
210
IF( LSAME( TRANS, 'N' ) )THEN
212
* Form C := alpha*A*A' + beta*C.
216
IF( BETA.EQ.ZERO )THEN
220
ELSE IF( BETA.NE.ONE )THEN
222
C( I, J ) = BETA*C( I, J )
226
IF( A( J, L ).NE.ZERO )THEN
227
TEMP = ALPHA*A( J, L )
229
C( I, J ) = C( I, J ) + TEMP*A( I, L )
236
IF( BETA.EQ.ZERO )THEN
240
ELSE IF( BETA.NE.ONE )THEN
242
C( I, J ) = BETA*C( I, J )
246
IF( A( J, L ).NE.ZERO )THEN
247
TEMP = ALPHA*A( J, L )
249
C( I, J ) = C( I, J ) + TEMP*A( I, L )
257
* Form C := alpha*A'*A + beta*C.
264
TEMP = TEMP + A( L, I )*A( L, J )
266
IF( BETA.EQ.ZERO )THEN
267
C( I, J ) = ALPHA*TEMP
269
C( I, J ) = ALPHA*TEMP + BETA*C( I, J )
278
TEMP = TEMP + A( L, I )*A( L, J )
280
IF( BETA.EQ.ZERO )THEN
281
C( I, J ) = ALPHA*TEMP
283
C( I, J ) = ALPHA*TEMP + BETA*C( I, J )