ergo
template_blas_symm.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_SYMM_HEADER
38 #define TEMPLATE_BLAS_SYMM_HEADER
39 
40 
41 template<class Treal>
42 int 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 integer
Definition: template_blas_common.h:40
#define maxMACRO(a, b)
Definition: template_blas_common.h:45
#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
int template_blas_erbla(const char *srname, integer *info)
Definition: template_blas_common.cc:146
#define c___ref(a_1, a_2)
#define a_ref(a_1, a_2)
bool logical
Definition: template_blas_common.h:41
side
Definition: Matrix.h:75
logical template_blas_lsame(const char *ca, const char *cb)
Definition: template_blas_common.cc:46