ergo
template_blas_gemv.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_GEMV_HEADER
38 #define TEMPLATE_BLAS_GEMV_HEADER
39 
40 #include "template_blas_common.h"
41 
42 template<class Treal>
43 int 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 integer
Definition: template_blas_common.h:40
#define maxMACRO(a, b)
Definition: template_blas_common.h:45
int template_blas_erbla(const char *srname, integer *info)
Definition: template_blas_common.cc:146
#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
logical template_blas_lsame(const char *ca, const char *cb)
Definition: template_blas_common.cc:46