ergo
template_blas_syr2k.h
Go to the documentation of this file.
1 /* Ergo, version 3.7, a program for linear scaling electronic structure
2  * calculations.
3  * Copyright (C) 2018 Elias Rudberg, Emanuel H. Rubensson, Pawel Salek,
4  * and Anastasia Kruchinina.
5  *
6  * This program is free software: you can redistribute it and/or modify
7  * it under the terms of the GNU General Public License as published by
8  * the Free Software Foundation, either version 3 of the License, or
9  * (at your option) any later version.
10  *
11  * This program is distributed in the hope that it will be useful,
12  * but WITHOUT ANY WARRANTY; without even the implied warranty of
13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14  * GNU General Public License for more details.
15  *
16  * You should have received a copy of the GNU General Public License
17  * along with this program. If not, see <http://www.gnu.org/licenses/>.
18  *
19  * Primary academic reference:
20  * Ergo: An open-source program for linear-scaling electronic structure
21  * calculations,
22  * Elias Rudberg, Emanuel H. Rubensson, Pawel Salek, and Anastasia
23  * Kruchinina,
24  * SoftwareX 7, 107 (2018),
25  * <http://dx.doi.org/10.1016/j.softx.2018.03.005>
26  *
27  * For further information about Ergo, see <http://www.ergoscf.org>.
28  */
29 
30  /* This file belongs to the template_lapack part of the Ergo source
31  * code. The source files in the template_lapack directory are modified
32  * versions of files originally distributed as CLAPACK, see the
33  * Copyright/license notice in the file template_lapack/COPYING.
34  */
35 
36 
37 #ifndef TEMPLATE_BLAS_SYR2K_HEADER
38 #define TEMPLATE_BLAS_SYR2K_HEADER
39 
40 
41 template<class Treal>
42 int template_blas_syr2k(const char *uplo, const char *trans, const integer *n,
43  const integer *k, const Treal *alpha, const Treal *a,
44  const integer *lda, const Treal *b, const integer *ldb,
45  const Treal *beta, Treal *c__, const integer *ldc)
46 {
47  /* System generated locals */
48  integer a_dim1, a_offset, b_dim1, b_offset, c_dim1, c_offset, i__1, i__2,
49  i__3;
50  /* Local variables */
51  integer info;
52  Treal temp1, temp2;
53  integer i__, j, l;
54  integer nrowa;
55  logical upper;
56 #define a_ref(a_1,a_2) a[(a_2)*a_dim1 + a_1]
57 #define b_ref(a_1,a_2) b[(a_2)*b_dim1 + a_1]
58 #define c___ref(a_1,a_2) c__[(a_2)*c_dim1 + a_1]
59 /* Purpose
60  =======
61  DSYR2K performs one of the symmetric rank 2k operations
62  C := alpha*A*B' + alpha*B*A' + beta*C,
63  or
64  C := alpha*A'*B + alpha*B'*A + beta*C,
65  where alpha and beta are scalars, C is an n by n symmetric matrix
66  and A and B are n by k matrices in the first case and k by n
67  matrices in the second case.
68  Parameters
69  ==========
70  UPLO - CHARACTER*1.
71  On entry, UPLO specifies whether the upper or lower
72  triangular part of the array C is to be referenced as
73  follows:
74  UPLO = 'U' or 'u' Only the upper triangular part of C
75  is to be referenced.
76  UPLO = 'L' or 'l' Only the lower triangular part of C
77  is to be referenced.
78  Unchanged on exit.
79  TRANS - CHARACTER*1.
80  On entry, TRANS specifies the operation to be performed as
81  follows:
82  TRANS = 'N' or 'n' C := alpha*A*B' + alpha*B*A' +
83  beta*C.
84  TRANS = 'T' or 't' C := alpha*A'*B + alpha*B'*A +
85  beta*C.
86  TRANS = 'C' or 'c' C := alpha*A'*B + alpha*B'*A +
87  beta*C.
88  Unchanged on exit.
89  N - INTEGER.
90  On entry, N specifies the order of the matrix C. N must be
91  at least zero.
92  Unchanged on exit.
93  K - INTEGER.
94  On entry with TRANS = 'N' or 'n', K specifies the number
95  of columns of the matrices A and B, and on entry with
96  TRANS = 'T' or 't' or 'C' or 'c', K specifies the number
97  of rows of the matrices A and B. K must be at least zero.
98  Unchanged on exit.
99  ALPHA - DOUBLE PRECISION.
100  On entry, ALPHA specifies the scalar alpha.
101  Unchanged on exit.
102  A - DOUBLE PRECISION array of DIMENSION ( LDA, ka ), where ka is
103  k when TRANS = 'N' or 'n', and is n otherwise.
104  Before entry with TRANS = 'N' or 'n', the leading n by k
105  part of the array A must contain the matrix A, otherwise
106  the leading k by n part of the array A must contain the
107  matrix A.
108  Unchanged on exit.
109  LDA - INTEGER.
110  On entry, LDA specifies the first dimension of A as declared
111  in the calling (sub) program. When TRANS = 'N' or 'n'
112  then LDA must be at least max( 1, n ), otherwise LDA must
113  be at least max( 1, k ).
114  Unchanged on exit.
115  B - DOUBLE PRECISION array of DIMENSION ( LDB, kb ), where kb is
116  k when TRANS = 'N' or 'n', and is n otherwise.
117  Before entry with TRANS = 'N' or 'n', the leading n by k
118  part of the array B must contain the matrix B, otherwise
119  the leading k by n part of the array B must contain the
120  matrix B.
121  Unchanged on exit.
122  LDB - INTEGER.
123  On entry, LDB specifies the first dimension of B as declared
124  in the calling (sub) program. When TRANS = 'N' or 'n'
125  then LDB must be at least max( 1, n ), otherwise LDB must
126  be at least max( 1, k ).
127  Unchanged on exit.
128  BETA - DOUBLE PRECISION.
129  On entry, BETA specifies the scalar beta.
130  Unchanged on exit.
131  C - DOUBLE PRECISION array of DIMENSION ( LDC, n ).
132  Before entry with UPLO = 'U' or 'u', the leading n by n
133  upper triangular part of the array C must contain the upper
134  triangular part of the symmetric matrix and the strictly
135  lower triangular part of C is not referenced. On exit, the
136  upper triangular part of the array C is overwritten by the
137  upper triangular part of the updated matrix.
138  Before entry with UPLO = 'L' or 'l', the leading n by n
139  lower triangular part of the array C must contain the lower
140  triangular part of the symmetric matrix and the strictly
141  upper triangular part of C is not referenced. On exit, the
142  lower triangular part of the array C is overwritten by the
143  lower triangular part of the updated matrix.
144  LDC - INTEGER.
145  On entry, LDC specifies the first dimension of C as declared
146  in the calling (sub) program. LDC must be at least
147  max( 1, n ).
148  Unchanged on exit.
149  Level 3 Blas routine.
150  -- Written on 8-February-1989.
151  Jack Dongarra, Argonne National Laboratory.
152  Iain Duff, AERE Harwell.
153  Jeremy Du Croz, Numerical Algorithms Group Ltd.
154  Sven Hammarling, Numerical Algorithms Group Ltd.
155  Test the input parameters.
156  Parameter adjustments */
157  a_dim1 = *lda;
158  a_offset = 1 + a_dim1 * 1;
159  a -= a_offset;
160  b_dim1 = *ldb;
161  b_offset = 1 + b_dim1 * 1;
162  b -= b_offset;
163  c_dim1 = *ldc;
164  c_offset = 1 + c_dim1 * 1;
165  c__ -= c_offset;
166  /* Function Body */
167  if (template_blas_lsame(trans, "N")) {
168  nrowa = *n;
169  } else {
170  nrowa = *k;
171  }
172  upper = template_blas_lsame(uplo, "U");
173  info = 0;
174  if (! upper && ! template_blas_lsame(uplo, "L")) {
175  info = 1;
176  } else if (! template_blas_lsame(trans, "N") && ! template_blas_lsame(trans,
177  "T") && ! template_blas_lsame(trans, "C")) {
178  info = 2;
179  } else if (*n < 0) {
180  info = 3;
181  } else if (*k < 0) {
182  info = 4;
183  } else if (*lda < maxMACRO(1,nrowa)) {
184  info = 7;
185  } else if (*ldb < maxMACRO(1,nrowa)) {
186  info = 9;
187  } else if (*ldc < maxMACRO(1,*n)) {
188  info = 12;
189  }
190  if (info != 0) {
191  template_blas_erbla("SYR2K ", &info);
192  return 0;
193  }
194 /* Quick return if possible. */
195  if (*n == 0 || ( (*alpha == 0. || *k == 0) && *beta == 1. ) ) {
196  return 0;
197  }
198 /* And when alpha.eq.zero. */
199  if (*alpha == 0.) {
200  if (upper) {
201  if (*beta == 0.) {
202  i__1 = *n;
203  for (j = 1; j <= i__1; ++j) {
204  i__2 = j;
205  for (i__ = 1; i__ <= i__2; ++i__) {
206  c___ref(i__, j) = 0.;
207 /* L10: */
208  }
209 /* L20: */
210  }
211  } else {
212  i__1 = *n;
213  for (j = 1; j <= i__1; ++j) {
214  i__2 = j;
215  for (i__ = 1; i__ <= i__2; ++i__) {
216  c___ref(i__, j) = *beta * c___ref(i__, j);
217 /* L30: */
218  }
219 /* L40: */
220  }
221  }
222  } else {
223  if (*beta == 0.) {
224  i__1 = *n;
225  for (j = 1; j <= i__1; ++j) {
226  i__2 = *n;
227  for (i__ = j; i__ <= i__2; ++i__) {
228  c___ref(i__, j) = 0.;
229 /* L50: */
230  }
231 /* L60: */
232  }
233  } else {
234  i__1 = *n;
235  for (j = 1; j <= i__1; ++j) {
236  i__2 = *n;
237  for (i__ = j; i__ <= i__2; ++i__) {
238  c___ref(i__, j) = *beta * c___ref(i__, j);
239 /* L70: */
240  }
241 /* L80: */
242  }
243  }
244  }
245  return 0;
246  }
247 /* Start the operations. */
248  if (template_blas_lsame(trans, "N")) {
249 /* Form C := alpha*A*B' + alpha*B*A' + C. */
250  if (upper) {
251  i__1 = *n;
252  for (j = 1; j <= i__1; ++j) {
253  if (*beta == 0.) {
254  i__2 = j;
255  for (i__ = 1; i__ <= i__2; ++i__) {
256  c___ref(i__, j) = 0.;
257 /* L90: */
258  }
259  } else if (*beta != 1.) {
260  i__2 = j;
261  for (i__ = 1; i__ <= i__2; ++i__) {
262  c___ref(i__, j) = *beta * c___ref(i__, j);
263 /* L100: */
264  }
265  }
266  i__2 = *k;
267  for (l = 1; l <= i__2; ++l) {
268  if (a_ref(j, l) != 0. || b_ref(j, l) != 0.) {
269  temp1 = *alpha * b_ref(j, l);
270  temp2 = *alpha * a_ref(j, l);
271  i__3 = j;
272  for (i__ = 1; i__ <= i__3; ++i__) {
273  c___ref(i__, j) = c___ref(i__, j) + a_ref(i__, l)
274  * temp1 + b_ref(i__, l) * temp2;
275 /* L110: */
276  }
277  }
278 /* L120: */
279  }
280 /* L130: */
281  }
282  } else {
283  i__1 = *n;
284  for (j = 1; j <= i__1; ++j) {
285  if (*beta == 0.) {
286  i__2 = *n;
287  for (i__ = j; i__ <= i__2; ++i__) {
288  c___ref(i__, j) = 0.;
289 /* L140: */
290  }
291  } else if (*beta != 1.) {
292  i__2 = *n;
293  for (i__ = j; i__ <= i__2; ++i__) {
294  c___ref(i__, j) = *beta * c___ref(i__, j);
295 /* L150: */
296  }
297  }
298  i__2 = *k;
299  for (l = 1; l <= i__2; ++l) {
300  if (a_ref(j, l) != 0. || b_ref(j, l) != 0.) {
301  temp1 = *alpha * b_ref(j, l);
302  temp2 = *alpha * a_ref(j, l);
303  i__3 = *n;
304  for (i__ = j; i__ <= i__3; ++i__) {
305  c___ref(i__, j) = c___ref(i__, j) + a_ref(i__, l)
306  * temp1 + b_ref(i__, l) * temp2;
307 /* L160: */
308  }
309  }
310 /* L170: */
311  }
312 /* L180: */
313  }
314  }
315  } else {
316 /* Form C := alpha*A'*B + alpha*B'*A + C. */
317  if (upper) {
318  i__1 = *n;
319  for (j = 1; j <= i__1; ++j) {
320  i__2 = j;
321  for (i__ = 1; i__ <= i__2; ++i__) {
322  temp1 = 0.;
323  temp2 = 0.;
324  i__3 = *k;
325  for (l = 1; l <= i__3; ++l) {
326  temp1 += a_ref(l, i__) * b_ref(l, j);
327  temp2 += b_ref(l, i__) * a_ref(l, j);
328 /* L190: */
329  }
330  if (*beta == 0.) {
331  c___ref(i__, j) = *alpha * temp1 + *alpha * temp2;
332  } else {
333  c___ref(i__, j) = *beta * c___ref(i__, j) + *alpha *
334  temp1 + *alpha * temp2;
335  }
336 /* L200: */
337  }
338 /* L210: */
339  }
340  } else {
341  i__1 = *n;
342  for (j = 1; j <= i__1; ++j) {
343  i__2 = *n;
344  for (i__ = j; i__ <= i__2; ++i__) {
345  temp1 = 0.;
346  temp2 = 0.;
347  i__3 = *k;
348  for (l = 1; l <= i__3; ++l) {
349  temp1 += a_ref(l, i__) * b_ref(l, j);
350  temp2 += b_ref(l, i__) * a_ref(l, j);
351 /* L220: */
352  }
353  if (*beta == 0.) {
354  c___ref(i__, j) = *alpha * temp1 + *alpha * temp2;
355  } else {
356  c___ref(i__, j) = *beta * c___ref(i__, j) + *alpha *
357  temp1 + *alpha * temp2;
358  }
359 /* L230: */
360  }
361 /* L240: */
362  }
363  }
364  }
365  return 0;
366 /* End of DSYR2K. */
367 } /* dsyr2k_ */
368 #undef c___ref
369 #undef b_ref
370 #undef a_ref
371 
372 #endif
int integer
Definition: template_blas_common.h:40
#define maxMACRO(a, b)
Definition: template_blas_common.h:45
#define b_ref(a_1, a_2)
#define a_ref(a_1, a_2)
#define c___ref(a_1, a_2)
int template_blas_erbla(const char *srname, integer *info)
Definition: template_blas_common.cc:146
bool logical
Definition: template_blas_common.h:41
int template_blas_syr2k(const char *uplo, const char *trans, const integer *n, const integer *k, const Treal *alpha, const Treal *a, const integer *lda, const Treal *b, const integer *ldb, const Treal *beta, Treal *c__, const integer *ldc)
Definition: template_blas_syr2k.h:42
logical template_blas_lsame(const char *ca, const char *cb)
Definition: template_blas_common.cc:46