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