ergo
template_blas_spmv.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_SPMV_HEADER
38#define TEMPLATE_BLAS_SPMV_HEADER
39
40
41template<class Treal>
42int 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
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
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