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