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