KTH framework for Nek5000 toolboxes; testing version  0.0.1
dsyr2k.f
Go to the documentation of this file.
1  SUBROUTINE dsyr2k( UPLO, TRANS, N, K, ALPHA, A, LDA, B, LDB,
2  $ BETA, C, LDC )
3 * .. Scalar Arguments ..
4  CHARACTER*1 UPLO, TRANS
5  INTEGER N, K, LDA, LDB, LDC
6  DOUBLE PRECISION ALPHA, BETA
7 * .. Array Arguments ..
8  DOUBLE PRECISION A( LDA, * ), B( LDB, * ), C( LDC, * )
9 * ..
10 *
11 * Purpose
12 * =======
13 *
14 * DSYR2K performs one of the symmetric rank 2k operations
15 *
16 * C := alpha*A*B' + alpha*B*A' + beta*C,
17 *
18 * or
19 *
20 * C := alpha*A'*B + alpha*B'*A + beta*C,
21 *
22 * where alpha and beta are scalars, C is an n by n symmetric matrix
23 * and A and B are n by k matrices in the first case and k by n
24 * matrices 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*B' + alpha*B*A' +
47 * beta*C.
48 *
49 * TRANS = 'T' or 't' C := alpha*A'*B + alpha*B'*A +
50 * beta*C.
51 *
52 * TRANS = 'C' or 'c' C := alpha*A'*B + alpha*B'*A +
53 * beta*C.
54 *
55 * Unchanged on exit.
56 *
57 * N - INTEGER.
58 * On entry, N specifies the order of the matrix C. N must be
59 * at least zero.
60 * Unchanged on exit.
61 *
62 * K - INTEGER.
63 * On entry with TRANS = 'N' or 'n', K specifies the number
64 * of columns of the matrices A and B, and on entry with
65 * TRANS = 'T' or 't' or 'C' or 'c', K specifies the number
66 * of rows of the matrices A and B. K must be at least zero.
67 * Unchanged on exit.
68 *
69 * ALPHA - DOUBLE PRECISION.
70 * On entry, ALPHA specifies the scalar alpha.
71 * Unchanged on exit.
72 *
73 * A - DOUBLE PRECISION array of DIMENSION ( LDA, ka ), where ka is
74 * k when TRANS = 'N' or 'n', and is n otherwise.
75 * Before entry with TRANS = 'N' or 'n', the leading n by k
76 * part of the array A must contain the matrix A, otherwise
77 * the leading k by n part of the array A must contain the
78 * matrix A.
79 * Unchanged on exit.
80 *
81 * LDA - INTEGER.
82 * On entry, LDA specifies the first dimension of A as declared
83 * in the calling (sub) program. When TRANS = 'N' or 'n'
84 * then LDA must be at least max( 1, n ), otherwise LDA must
85 * be at least max( 1, k ).
86 * Unchanged on exit.
87 *
88 * B - DOUBLE PRECISION array of DIMENSION ( LDB, kb ), where kb is
89 * k when TRANS = 'N' or 'n', and is n otherwise.
90 * Before entry with TRANS = 'N' or 'n', the leading n by k
91 * part of the array B must contain the matrix B, otherwise
92 * the leading k by n part of the array B must contain the
93 * matrix B.
94 * Unchanged on exit.
95 *
96 * LDB - INTEGER.
97 * On entry, LDB specifies the first dimension of B as declared
98 * in the calling (sub) program. When TRANS = 'N' or 'n'
99 * then LDB must be at least max( 1, n ), otherwise LDB must
100 * be at least max( 1, k ).
101 * Unchanged on exit.
102 *
103 * BETA - DOUBLE PRECISION.
104 * On entry, BETA specifies the scalar beta.
105 * Unchanged on exit.
106 *
107 * C - DOUBLE PRECISION array of DIMENSION ( LDC, n ).
108 * Before entry with UPLO = 'U' or 'u', the leading n by n
109 * upper triangular part of the array C must contain the upper
110 * triangular part of the symmetric matrix and the strictly
111 * lower triangular part of C is not referenced. On exit, the
112 * upper triangular part of the array C is overwritten by the
113 * upper triangular part of the updated matrix.
114 * Before entry with UPLO = 'L' or 'l', the leading n by n
115 * lower triangular part of the array C must contain the lower
116 * triangular part of the symmetric matrix and the strictly
117 * upper triangular part of C is not referenced. On exit, the
118 * lower triangular part of the array C is overwritten by the
119 * lower triangular part of the updated matrix.
120 *
121 * LDC - INTEGER.
122 * On entry, LDC specifies the first dimension of C as declared
123 * in the calling (sub) program. LDC must be at least
124 * max( 1, n ).
125 * Unchanged on exit.
126 *
127 *
128 * Level 3 Blas routine.
129 *
130 *
131 * -- Written on 8-February-1989.
132 * Jack Dongarra, Argonne National Laboratory.
133 * Iain Duff, AERE Harwell.
134 * Jeremy Du Croz, Numerical Algorithms Group Ltd.
135 * Sven Hammarling, Numerical Algorithms Group Ltd.
136 *
137 *
138 * .. External Functions ..
139  LOGICAL LSAME
140  EXTERNAL lsame
141 * .. External Subroutines ..
142  EXTERNAL xerbla
143 * .. Intrinsic Functions ..
144  INTRINSIC max
145 * .. Local Scalars ..
146  LOGICAL UPPER
147  INTEGER I, INFO, J, L, NROWA
148  DOUBLE PRECISION TEMP1, TEMP2
149 * .. Parameters ..
150  DOUBLE PRECISION ONE , ZERO
151  parameter( one = 1.0d+0, zero = 0.0d+0 )
152 * ..
153 * .. Executable Statements ..
154 *
155 * Test the input parameters.
156 *
157  IF( lsame( trans, 'N' ) )THEN
158  nrowa = n
159  ELSE
160  nrowa = k
161  END IF
162  upper = lsame( uplo, 'U' )
163 *
164  info = 0
165  IF( ( .NOT.upper ).AND.
166  $ ( .NOT.lsame( uplo , 'L' ) ) )THEN
167  info = 1
168  ELSE IF( ( .NOT.lsame( trans, 'N' ) ).AND.
169  $ ( .NOT.lsame( trans, 'T' ) ).AND.
170  $ ( .NOT.lsame( trans, 'C' ) ) )THEN
171  info = 2
172  ELSE IF( n .LT.0 )THEN
173  info = 3
174  ELSE IF( k .LT.0 )THEN
175  info = 4
176  ELSE IF( lda.LT.max( 1, nrowa ) )THEN
177  info = 7
178  ELSE IF( ldb.LT.max( 1, nrowa ) )THEN
179  info = 9
180  ELSE IF( ldc.LT.max( 1, n ) )THEN
181  info = 12
182  END IF
183  IF( info.NE.0 )THEN
184  CALL xerbla( 'DSYR2K', info )
185  RETURN
186  END IF
187 *
188 * Quick return if possible.
189 *
190  IF( ( n.EQ.0 ).OR.
191  $ ( ( ( alpha.EQ.zero ).OR.( k.EQ.0 ) ).AND.( beta.EQ.one ) ) )
192  $ RETURN
193 *
194 * And when alpha.eq.zero.
195 *
196  IF( alpha.EQ.zero )THEN
197  IF( upper )THEN
198  IF( beta.EQ.zero )THEN
199  DO 20, j = 1, n
200  DO 10, i = 1, j
201  c( i, j ) = zero
202  10 CONTINUE
203  20 CONTINUE
204  ELSE
205  DO 40, j = 1, n
206  DO 30, i = 1, j
207  c( i, j ) = beta*c( i, j )
208  30 CONTINUE
209  40 CONTINUE
210  END IF
211  ELSE
212  IF( beta.EQ.zero )THEN
213  DO 60, j = 1, n
214  DO 50, i = j, n
215  c( i, j ) = zero
216  50 CONTINUE
217  60 CONTINUE
218  ELSE
219  DO 80, j = 1, n
220  DO 70, i = j, n
221  c( i, j ) = beta*c( i, j )
222  70 CONTINUE
223  80 CONTINUE
224  END IF
225  END IF
226  RETURN
227  END IF
228 *
229 * Start the operations.
230 *
231  IF( lsame( trans, 'N' ) )THEN
232 *
233 * Form C := alpha*A*B' + alpha*B*A' + C.
234 *
235  IF( upper )THEN
236  DO 130, j = 1, n
237  IF( beta.EQ.zero )THEN
238  DO 90, i = 1, j
239  c( i, j ) = zero
240  90 CONTINUE
241  ELSE IF( beta.NE.one )THEN
242  DO 100, i = 1, j
243  c( i, j ) = beta*c( i, j )
244  100 CONTINUE
245  END IF
246  DO 120, l = 1, k
247  IF( ( a( j, l ).NE.zero ).OR.
248  $ ( b( j, l ).NE.zero ) )THEN
249  temp1 = alpha*b( j, l )
250  temp2 = alpha*a( j, l )
251  DO 110, i = 1, j
252  c( i, j ) = c( i, j ) +
253  $ a( i, l )*temp1 + b( i, l )*temp2
254  110 CONTINUE
255  END IF
256  120 CONTINUE
257  130 CONTINUE
258  ELSE
259  DO 180, j = 1, n
260  IF( beta.EQ.zero )THEN
261  DO 140, i = j, n
262  c( i, j ) = zero
263  140 CONTINUE
264  ELSE IF( beta.NE.one )THEN
265  DO 150, i = j, n
266  c( i, j ) = beta*c( i, j )
267  150 CONTINUE
268  END IF
269  DO 170, l = 1, k
270  IF( ( a( j, l ).NE.zero ).OR.
271  $ ( b( j, l ).NE.zero ) )THEN
272  temp1 = alpha*b( j, l )
273  temp2 = alpha*a( j, l )
274  DO 160, i = j, n
275  c( i, j ) = c( i, j ) +
276  $ a( i, l )*temp1 + b( i, l )*temp2
277  160 CONTINUE
278  END IF
279  170 CONTINUE
280  180 CONTINUE
281  END IF
282  ELSE
283 *
284 * Form C := alpha*A'*B + alpha*B'*A + C.
285 *
286  IF( upper )THEN
287  DO 210, j = 1, n
288  DO 200, i = 1, j
289  temp1 = zero
290  temp2 = zero
291  DO 190, l = 1, k
292  temp1 = temp1 + a( l, i )*b( l, j )
293  temp2 = temp2 + b( l, i )*a( l, j )
294  190 CONTINUE
295  IF( beta.EQ.zero )THEN
296  c( i, j ) = alpha*temp1 + alpha*temp2
297  ELSE
298  c( i, j ) = beta *c( i, j ) +
299  $ alpha*temp1 + alpha*temp2
300  END IF
301  200 CONTINUE
302  210 CONTINUE
303  ELSE
304  DO 240, j = 1, n
305  DO 230, i = j, n
306  temp1 = zero
307  temp2 = zero
308  DO 220, l = 1, k
309  temp1 = temp1 + a( l, i )*b( l, j )
310  temp2 = temp2 + b( l, i )*a( l, j )
311  220 CONTINUE
312  IF( beta.EQ.zero )THEN
313  c( i, j ) = alpha*temp1 + alpha*temp2
314  ELSE
315  c( i, j ) = beta *c( i, j ) +
316  $ alpha*temp1 + alpha*temp2
317  END IF
318  230 CONTINUE
319  240 CONTINUE
320  END IF
321  END IF
322 *
323  RETURN
324 *
325 * End of DSYR2K.
326 *
327  END
dsyr2k
subroutine dsyr2k(UPLO, TRANS, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
Definition: dsyr2k.f:3
xerbla
subroutine xerbla(SRNAME, INFO)
Definition: xerbla.f:2