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