ergo
template_blas_syr2k.h
Go to the documentation of this file.
1/* Ergo, version 3.8.2, a program for linear scaling electronic structure
2 * calculations.
3 * Copyright (C) 2023 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
41template<class Treal>
42int 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 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_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
#define c___ref(a_1, a_2)
#define a_ref(a_1, a_2)
#define b_ref(a_1, a_2)