ergo
template_blas_symv.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_SYMV_HEADER
38 #define TEMPLATE_BLAS_SYMV_HEADER
39 
40 
41 template<class Treal>
42 int template_blas_symv(const char *uplo, const integer *n, const Treal *alpha,
43  const Treal *a, const integer *lda, const Treal *x, const integer *incx, const Treal
44  *beta, Treal *y, const integer *incy)
45 {
46  /* System generated locals */
47  integer a_dim1, a_offset, i__1, i__2;
48  /* Local variables */
49  integer info;
50  Treal temp1, temp2;
51  integer i__, j;
52  integer ix, iy, jx, jy, kx, ky;
53 #define a_ref(a_1,a_2) a[(a_2)*a_dim1 + a_1]
54 /* Purpose
55  =======
56  DSYMV performs the matrix-vector operation
57  y := alpha*A*x + beta*y,
58  where alpha and beta are scalars, x and y are n element vectors and
59  A is an n by n symmetric matrix.
60  Parameters
61  ==========
62  UPLO - CHARACTER*1.
63  On entry, UPLO specifies whether the upper or lower
64  triangular part of the array A is to be referenced as
65  follows:
66  UPLO = 'U' or 'u' Only the upper triangular part of A
67  is to be referenced.
68  UPLO = 'L' or 'l' Only the lower triangular part of A
69  is to be referenced.
70  Unchanged on exit.
71  N - INTEGER.
72  On entry, N specifies the order of the matrix A.
73  N must be at least zero.
74  Unchanged on exit.
75  ALPHA - DOUBLE PRECISION.
76  On entry, ALPHA specifies the scalar alpha.
77  Unchanged on exit.
78  A - DOUBLE PRECISION array of DIMENSION ( LDA, n ).
79  Before entry with UPLO = 'U' or 'u', the leading n by n
80  upper triangular part of the array A must contain the upper
81  triangular part of the symmetric matrix and the strictly
82  lower triangular part of A is not referenced.
83  Before entry with UPLO = 'L' or 'l', the leading n by n
84  lower triangular part of the array A must contain the lower
85  triangular part of the symmetric matrix and the strictly
86  upper triangular part of A is not referenced.
87  Unchanged on exit.
88  LDA - INTEGER.
89  On entry, LDA specifies the first dimension of A as declared
90  in the calling (sub) program. LDA must be at least
91  max( 1, n ).
92  Unchanged on exit.
93  X - DOUBLE PRECISION array of dimension at least
94  ( 1 + ( n - 1 )*abs( INCX ) ).
95  Before entry, the incremented array X must contain the n
96  element vector x.
97  Unchanged on exit.
98  INCX - INTEGER.
99  On entry, INCX specifies the increment for the elements of
100  X. INCX must not be zero.
101  Unchanged on exit.
102  BETA - DOUBLE PRECISION.
103  On entry, BETA specifies the scalar beta. When BETA is
104  supplied as zero then Y need not be set on input.
105  Unchanged on exit.
106  Y - DOUBLE PRECISION array of dimension at least
107  ( 1 + ( n - 1 )*abs( INCY ) ).
108  Before entry, the incremented array Y must contain the n
109  element vector y. On exit, Y is overwritten by the updated
110  vector y.
111  INCY - INTEGER.
112  On entry, INCY specifies the increment for the elements of
113  Y. INCY must not be zero.
114  Unchanged on exit.
115  Level 2 Blas routine.
116  -- Written on 22-October-1986.
117  Jack Dongarra, Argonne National Lab.
118  Jeremy Du Croz, Nag Central Office.
119  Sven Hammarling, Nag Central Office.
120  Richard Hanson, Sandia National Labs.
121  Test the input parameters.
122  Parameter adjustments */
123  a_dim1 = *lda;
124  a_offset = 1 + a_dim1 * 1;
125  a -= a_offset;
126  --x;
127  --y;
128  /* Function Body */
129  info = 0;
130  if (! template_blas_lsame(uplo, "U") && ! template_blas_lsame(uplo, "L")) {
131  info = 1;
132  } else if (*n < 0) {
133  info = 2;
134  } else if (*lda < maxMACRO(1,*n)) {
135  info = 5;
136  } else if (*incx == 0) {
137  info = 7;
138  } else if (*incy == 0) {
139  info = 10;
140  }
141  if (info != 0) {
142  template_blas_erbla("SYMV ", &info);
143  return 0;
144  }
145 /* Quick return if possible. */
146  if (*n == 0 || (*alpha == 0. && *beta == 1.) ) {
147  return 0;
148  }
149 /* Set up the start points in X and Y. */
150  if (*incx > 0) {
151  kx = 1;
152  } else {
153  kx = 1 - (*n - 1) * *incx;
154  }
155  if (*incy > 0) {
156  ky = 1;
157  } else {
158  ky = 1 - (*n - 1) * *incy;
159  }
160 /* Start the operations. In this version the elements of A are
161  accessed sequentially with one pass through the triangular part
162  of A.
163  First form y := beta*y. */
164  if (*beta != 1.) {
165  if (*incy == 1) {
166  if (*beta == 0.) {
167  i__1 = *n;
168  for (i__ = 1; i__ <= i__1; ++i__) {
169  y[i__] = 0.;
170 /* L10: */
171  }
172  } else {
173  i__1 = *n;
174  for (i__ = 1; i__ <= i__1; ++i__) {
175  y[i__] = *beta * y[i__];
176 /* L20: */
177  }
178  }
179  } else {
180  iy = ky;
181  if (*beta == 0.) {
182  i__1 = *n;
183  for (i__ = 1; i__ <= i__1; ++i__) {
184  y[iy] = 0.;
185  iy += *incy;
186 /* L30: */
187  }
188  } else {
189  i__1 = *n;
190  for (i__ = 1; i__ <= i__1; ++i__) {
191  y[iy] = *beta * y[iy];
192  iy += *incy;
193 /* L40: */
194  }
195  }
196  }
197  }
198  if (*alpha == 0.) {
199  return 0;
200  }
201  if (template_blas_lsame(uplo, "U")) {
202 /* Form y when A is stored in upper triangle. */
203  if (*incx == 1 && *incy == 1) {
204  i__1 = *n;
205  for (j = 1; j <= i__1; ++j) {
206  temp1 = *alpha * x[j];
207  temp2 = 0.;
208  i__2 = j - 1;
209  for (i__ = 1; i__ <= i__2; ++i__) {
210  y[i__] += temp1 * a_ref(i__, j);
211  temp2 += a_ref(i__, j) * x[i__];
212 /* L50: */
213  }
214  y[j] = y[j] + temp1 * a_ref(j, j) + *alpha * temp2;
215 /* L60: */
216  }
217  } else {
218  jx = kx;
219  jy = ky;
220  i__1 = *n;
221  for (j = 1; j <= i__1; ++j) {
222  temp1 = *alpha * x[jx];
223  temp2 = 0.;
224  ix = kx;
225  iy = ky;
226  i__2 = j - 1;
227  for (i__ = 1; i__ <= i__2; ++i__) {
228  y[iy] += temp1 * a_ref(i__, j);
229  temp2 += a_ref(i__, j) * x[ix];
230  ix += *incx;
231  iy += *incy;
232 /* L70: */
233  }
234  y[jy] = y[jy] + temp1 * a_ref(j, j) + *alpha * temp2;
235  jx += *incx;
236  jy += *incy;
237 /* L80: */
238  }
239  }
240  } else {
241 /* Form y when A is stored in lower triangle. */
242  if (*incx == 1 && *incy == 1) {
243  i__1 = *n;
244  for (j = 1; j <= i__1; ++j) {
245  temp1 = *alpha * x[j];
246  temp2 = 0.;
247  y[j] += temp1 * a_ref(j, j);
248  i__2 = *n;
249  for (i__ = j + 1; i__ <= i__2; ++i__) {
250  y[i__] += temp1 * a_ref(i__, j);
251  temp2 += a_ref(i__, j) * x[i__];
252 /* L90: */
253  }
254  y[j] += *alpha * temp2;
255 /* L100: */
256  }
257  } else {
258  jx = kx;
259  jy = ky;
260  i__1 = *n;
261  for (j = 1; j <= i__1; ++j) {
262  temp1 = *alpha * x[jx];
263  temp2 = 0.;
264  y[jy] += temp1 * a_ref(j, j);
265  ix = jx;
266  iy = jy;
267  i__2 = *n;
268  for (i__ = j + 1; i__ <= i__2; ++i__) {
269  ix += *incx;
270  iy += *incy;
271  y[iy] += temp1 * a_ref(i__, j);
272  temp2 += a_ref(i__, j) * x[ix];
273 /* L110: */
274  }
275  y[jy] += *alpha * temp2;
276  jx += *incx;
277  jy += *incy;
278 /* L120: */
279  }
280  }
281  }
282  return 0;
283 /* End of DSYMV . */
284 } /* dsymv_ */
285 #undef a_ref
286 
287 #endif
int integer
Definition: template_blas_common.h:40
#define maxMACRO(a, b)
Definition: template_blas_common.h:45
#define a_ref(a_1, a_2)
int template_blas_erbla(const char *srname, integer *info)
Definition: template_blas_common.cc:146
int template_blas_symv(const char *uplo, const integer *n, const Treal *alpha, const Treal *a, const integer *lda, const Treal *x, const integer *incx, const Treal *beta, Treal *y, const integer *incy)
Definition: template_blas_symv.h:42
logical template_blas_lsame(const char *ca, const char *cb)
Definition: template_blas_common.cc:46