ergo
template_blas_syrk.h
Go to the documentation of this file.
1/* Ergo, version 3.8, a program for linear scaling electronic structure
2 * calculations.
3 * Copyright (C) 2019 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
41template<class Treal>
42int 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 template_blas_erbla(const char *srname, integer *info)
Definition: template_blas_common.cc:146
logical template_blas_lsame(const char *ca, const char *cb)
Definition: template_blas_common.cc:46
int integer
Definition: template_blas_common.h:40
#define maxMACRO(a, b)
Definition: template_blas_common.h:45
bool logical
Definition: template_blas_common.h:41
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
#define c___ref(a_1, a_2)
#define a_ref(a_1, a_2)