ergo
template_blas_syrk.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_SYRK_HEADER
38 #define TEMPLATE_BLAS_SYRK_HEADER
39 
40 
41 template<class Treal>
42 int template_blas_syrk(const char *uplo, const char *trans, const integer *n, const integer *k,
43  const Treal *alpha, const Treal *a, const integer *lda, const Treal *beta,
44  Treal *c__, const integer *ldc)
45 {
46  /* System generated locals */
47  integer a_dim1, a_offset, c_dim1, c_offset, i__1, i__2, i__3;
48  /* Local variables */
49  integer info;
50  Treal temp;
51  integer i__, j, l;
52  integer nrowa;
53  logical upper;
54 #define a_ref(a_1,a_2) a[(a_2)*a_dim1 + a_1]
55 #define c___ref(a_1,a_2) c__[(a_2)*c_dim1 + a_1]
56 /* Purpose
57  =======
58  DSYRK performs one of the symmetric rank k operations
59  C := alpha*A*A' + beta*C,
60  or
61  C := alpha*A'*A + beta*C,
62  where alpha and beta are scalars, C is an n by n symmetric matrix
63  and A is an n by k matrix in the first case and a k by n matrix
64  in the second case.
65  Parameters
66  ==========
67  UPLO - CHARACTER*1.
68  On entry, UPLO specifies whether the upper or lower
69  triangular part of the array C is to be referenced as
70  follows:
71  UPLO = 'U' or 'u' Only the upper triangular part of C
72  is to be referenced.
73  UPLO = 'L' or 'l' Only the lower triangular part of C
74  is to be referenced.
75  Unchanged on exit.
76  TRANS - CHARACTER*1.
77  On entry, TRANS specifies the operation to be performed as
78  follows:
79  TRANS = 'N' or 'n' C := alpha*A*A' + beta*C.
80  TRANS = 'T' or 't' C := alpha*A'*A + beta*C.
81  TRANS = 'C' or 'c' C := alpha*A'*A + beta*C.
82  Unchanged on exit.
83  N - INTEGER.
84  On entry, N specifies the order of the matrix C. N must be
85  at least zero.
86  Unchanged on exit.
87  K - INTEGER.
88  On entry with TRANS = 'N' or 'n', K specifies the number
89  of columns of the matrix A, and on entry with
90  TRANS = 'T' or 't' or 'C' or 'c', K specifies the number
91  of rows of the matrix A. K must be at least zero.
92  Unchanged on exit.
93  ALPHA - DOUBLE PRECISION.
94  On entry, ALPHA specifies the scalar alpha.
95  Unchanged on exit.
96  A - DOUBLE PRECISION array of DIMENSION ( LDA, ka ), where ka is
97  k when TRANS = 'N' or 'n', and is n otherwise.
98  Before entry with TRANS = 'N' or 'n', the leading n by k
99  part of the array A must contain the matrix A, otherwise
100  the leading k by n part of the array A must contain the
101  matrix A.
102  Unchanged on exit.
103  LDA - INTEGER.
104  On entry, LDA specifies the first dimension of A as declared
105  in the calling (sub) program. When TRANS = 'N' or 'n'
106  then LDA must be at least max( 1, n ), otherwise LDA must
107  be at least max( 1, k ).
108  Unchanged on exit.
109  BETA - DOUBLE PRECISION.
110  On entry, BETA specifies the scalar beta.
111  Unchanged on exit.
112  C - DOUBLE PRECISION array of DIMENSION ( LDC, n ).
113  Before entry with UPLO = 'U' or 'u', the leading n by n
114  upper triangular part of the array C must contain the upper
115  triangular part of the symmetric matrix and the strictly
116  lower triangular part of C is not referenced. On exit, the
117  upper triangular part of the array C is overwritten by the
118  upper triangular part of the updated matrix.
119  Before entry with UPLO = 'L' or 'l', the leading n by n
120  lower triangular part of the array C must contain the lower
121  triangular part of the symmetric matrix and the strictly
122  upper triangular part of C is not referenced. On exit, the
123  lower triangular part of the array C is overwritten by the
124  lower triangular part of the updated matrix.
125  LDC - INTEGER.
126  On entry, LDC specifies the first dimension of C as declared
127  in the calling (sub) program. LDC must be at least
128  max( 1, n ).
129  Unchanged on exit.
130  Level 3 Blas routine.
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  Test the input parameters.
137  Parameter adjustments */
138  a_dim1 = *lda;
139  a_offset = 1 + a_dim1 * 1;
140  a -= a_offset;
141  c_dim1 = *ldc;
142  c_offset = 1 + c_dim1 * 1;
143  c__ -= c_offset;
144  /* Function Body */
145  if (template_blas_lsame(trans, "N")) {
146  nrowa = *n;
147  } else {
148  nrowa = *k;
149  }
150  upper = template_blas_lsame(uplo, "U");
151  info = 0;
152  if (! upper && ! template_blas_lsame(uplo, "L")) {
153  info = 1;
154  } else if (! template_blas_lsame(trans, "N") && ! template_blas_lsame(trans,
155  "T") && ! template_blas_lsame(trans, "C")) {
156  info = 2;
157  } else if (*n < 0) {
158  info = 3;
159  } else if (*k < 0) {
160  info = 4;
161  } else if (*lda < maxMACRO(1,nrowa)) {
162  info = 7;
163  } else if (*ldc < maxMACRO(1,*n)) {
164  info = 10;
165  }
166  if (info != 0) {
167  template_blas_erbla("DSYRK ", &info);
168  return 0;
169  }
170 /* Quick return if possible. */
171  if (*n == 0 || ( (*alpha == 0. || *k == 0) && *beta == 1. ) ) {
172  return 0;
173  }
174 /* And when alpha.eq.zero. */
175  if (*alpha == 0.) {
176  if (upper) {
177  if (*beta == 0.) {
178  i__1 = *n;
179  for (j = 1; j <= i__1; ++j) {
180  i__2 = j;
181  for (i__ = 1; i__ <= i__2; ++i__) {
182  c___ref(i__, j) = 0.;
183 /* L10: */
184  }
185 /* L20: */
186  }
187  } else {
188  i__1 = *n;
189  for (j = 1; j <= i__1; ++j) {
190  i__2 = j;
191  for (i__ = 1; i__ <= i__2; ++i__) {
192  c___ref(i__, j) = *beta * c___ref(i__, j);
193 /* L30: */
194  }
195 /* L40: */
196  }
197  }
198  } else {
199  if (*beta == 0.) {
200  i__1 = *n;
201  for (j = 1; j <= i__1; ++j) {
202  i__2 = *n;
203  for (i__ = j; i__ <= i__2; ++i__) {
204  c___ref(i__, j) = 0.;
205 /* L50: */
206  }
207 /* L60: */
208  }
209  } else {
210  i__1 = *n;
211  for (j = 1; j <= i__1; ++j) {
212  i__2 = *n;
213  for (i__ = j; i__ <= i__2; ++i__) {
214  c___ref(i__, j) = *beta * c___ref(i__, j);
215 /* L70: */
216  }
217 /* L80: */
218  }
219  }
220  }
221  return 0;
222  }
223 /* Start the operations. */
224  if (template_blas_lsame(trans, "N")) {
225 /* Form C := alpha*A*A' + beta*C. */
226  if (upper) {
227  i__1 = *n;
228  for (j = 1; j <= i__1; ++j) {
229  if (*beta == 0.) {
230  i__2 = j;
231  for (i__ = 1; i__ <= i__2; ++i__) {
232  c___ref(i__, j) = 0.;
233 /* L90: */
234  }
235  } else if (*beta != 1.) {
236  i__2 = j;
237  for (i__ = 1; i__ <= i__2; ++i__) {
238  c___ref(i__, j) = *beta * c___ref(i__, j);
239 /* L100: */
240  }
241  }
242  i__2 = *k;
243  for (l = 1; l <= i__2; ++l) {
244  if (a_ref(j, l) != 0.) {
245  temp = *alpha * a_ref(j, l);
246  i__3 = j;
247  for (i__ = 1; i__ <= i__3; ++i__) {
248  c___ref(i__, j) = c___ref(i__, j) + temp * a_ref(
249  i__, l);
250 /* L110: */
251  }
252  }
253 /* L120: */
254  }
255 /* L130: */
256  }
257  } else {
258  i__1 = *n;
259  for (j = 1; j <= i__1; ++j) {
260  if (*beta == 0.) {
261  i__2 = *n;
262  for (i__ = j; i__ <= i__2; ++i__) {
263  c___ref(i__, j) = 0.;
264 /* L140: */
265  }
266  } else if (*beta != 1.) {
267  i__2 = *n;
268  for (i__ = j; i__ <= i__2; ++i__) {
269  c___ref(i__, j) = *beta * c___ref(i__, j);
270 /* L150: */
271  }
272  }
273  i__2 = *k;
274  for (l = 1; l <= i__2; ++l) {
275  if (a_ref(j, l) != 0.) {
276  temp = *alpha * a_ref(j, l);
277  i__3 = *n;
278  for (i__ = j; i__ <= i__3; ++i__) {
279  c___ref(i__, j) = c___ref(i__, j) + temp * a_ref(
280  i__, l);
281 /* L160: */
282  }
283  }
284 /* L170: */
285  }
286 /* L180: */
287  }
288  }
289  } else {
290 /* Form C := alpha*A'*A + beta*C. */
291  if (upper) {
292  i__1 = *n;
293  for (j = 1; j <= i__1; ++j) {
294  i__2 = j;
295  for (i__ = 1; i__ <= i__2; ++i__) {
296  temp = 0.;
297  i__3 = *k;
298  for (l = 1; l <= i__3; ++l) {
299  temp += a_ref(l, i__) * a_ref(l, j);
300 /* L190: */
301  }
302  if (*beta == 0.) {
303  c___ref(i__, j) = *alpha * temp;
304  } else {
305  c___ref(i__, j) = *alpha * temp + *beta * c___ref(i__,
306  j);
307  }
308 /* L200: */
309  }
310 /* L210: */
311  }
312  } else {
313  i__1 = *n;
314  for (j = 1; j <= i__1; ++j) {
315  i__2 = *n;
316  for (i__ = j; i__ <= i__2; ++i__) {
317  temp = 0.;
318  i__3 = *k;
319  for (l = 1; l <= i__3; ++l) {
320  temp += a_ref(l, i__) * a_ref(l, j);
321 /* L220: */
322  }
323  if (*beta == 0.) {
324  c___ref(i__, j) = *alpha * temp;
325  } else {
326  c___ref(i__, j) = *alpha * temp + *beta * c___ref(i__,
327  j);
328  }
329 /* L230: */
330  }
331 /* L240: */
332  }
333  }
334  }
335  return 0;
336 /* End of DSYRK . */
337 } /* dsyrk_ */
338 #undef c___ref
339 #undef a_ref
340 
341 #endif
int integer
Definition: template_blas_common.h:40
#define maxMACRO(a, b)
Definition: template_blas_common.h:45
int template_blas_syrk(const char *uplo, const char *trans, const integer *n, const integer *k, const Treal *alpha, const Treal *a, const integer *lda, const Treal *beta, Treal *c__, const integer *ldc)
Definition: template_blas_syrk.h:42
int template_blas_erbla(const char *srname, integer *info)
Definition: template_blas_common.cc:146
#define a_ref(a_1, a_2)
bool logical
Definition: template_blas_common.h:41
#define c___ref(a_1, a_2)
logical template_blas_lsame(const char *ca, const char *cb)
Definition: template_blas_common.cc:46