ergo
template_blas_symv.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_SYMV_HEADER
38#define TEMPLATE_BLAS_SYMV_HEADER
39
40
41template<class Treal>
42int 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 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
#define a_ref(a_1, a_2)
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