ergo
template_blas_gemv.h
Go to the documentation of this file.
1 /* Ergo, version 3.4, a program for linear scaling electronic structure
2  * calculations.
3  * Copyright (C) 2014 Elias Rudberg, Emanuel H. Rubensson, and Pawel Salek.
4  *
5  * This program is free software: you can redistribute it and/or modify
6  * it under the terms of the GNU General Public License as published by
7  * the Free Software Foundation, either version 3 of the License, or
8  * (at your option) any later version.
9  *
10  * This program is distributed in the hope that it will be useful,
11  * but WITHOUT ANY WARRANTY; without even the implied warranty of
12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13  * GNU General Public License for more details.
14  *
15  * You should have received a copy of the GNU General Public License
16  * along with this program. If not, see <http://www.gnu.org/licenses/>.
17  *
18  * Primary academic reference:
19  * Kohn−Sham Density Functional Theory Electronic Structure Calculations
20  * with Linearly Scaling Computational Time and Memory Usage,
21  * Elias Rudberg, Emanuel H. Rubensson, and Pawel Salek,
22  * J. Chem. Theory Comput. 7, 340 (2011),
23  * <http://dx.doi.org/10.1021/ct100611z>
24  *
25  * For further information about Ergo, see <http://www.ergoscf.org>.
26  */
27 
28  /* This file belongs to the template_lapack part of the Ergo source
29  * code. The source files in the template_lapack directory are modified
30  * versions of files originally distributed as CLAPACK, see the
31  * Copyright/license notice in the file template_lapack/COPYING.
32  */
33 
34 
35 #ifndef TEMPLATE_BLAS_GEMV_HEADER
36 #define TEMPLATE_BLAS_GEMV_HEADER
37 
38 #include "template_blas_common.h"
39 
40 template<class Treal>
41 int template_blas_gemv(const char *trans, const integer *m, const integer *n, const Treal *
42  alpha, const Treal *a, const integer *lda, const Treal *x, const integer *incx,
43  const Treal *beta, Treal *y, const integer *incy)
44 {
45  /* System generated locals */
46  integer a_dim1, a_offset, i__1, i__2;
47  /* Local variables */
48  integer info;
49  Treal temp;
50  integer lenx, leny, i__, j;
51  integer ix, iy, jx, jy, kx, ky;
52 #define a_ref(a_1,a_2) a[(a_2)*a_dim1 + a_1]
53 /* Purpose
54  =======
55  DGEMV performs one of the matrix-vector operations
56  y := alpha*A*x + beta*y, or y := alpha*A'*x + beta*y,
57  where alpha and beta are scalars, x and y are vectors and A is an
58  m by n matrix.
59  Parameters
60  ==========
61  TRANS - CHARACTER*1.
62  On entry, TRANS specifies the operation to be performed as
63  follows:
64  TRANS = 'N' or 'n' y := alpha*A*x + beta*y.
65  TRANS = 'T' or 't' y := alpha*A'*x + beta*y.
66  TRANS = 'C' or 'c' y := alpha*A'*x + beta*y.
67  Unchanged on exit.
68  M - INTEGER.
69  On entry, M specifies the number of rows of the matrix A.
70  M must be at least zero.
71  Unchanged on exit.
72  N - INTEGER.
73  On entry, N specifies the number of columns of the matrix A.
74  N must be at least zero.
75  Unchanged on exit.
76  ALPHA - DOUBLE PRECISION.
77  On entry, ALPHA specifies the scalar alpha.
78  Unchanged on exit.
79  A - DOUBLE PRECISION array of DIMENSION ( LDA, n ).
80  Before entry, the leading m by n part of the array A must
81  contain the matrix of coefficients.
82  Unchanged on exit.
83  LDA - INTEGER.
84  On entry, LDA specifies the first dimension of A as declared
85  in the calling (sub) program. LDA must be at least
86  max( 1, m ).
87  Unchanged on exit.
88  X - DOUBLE PRECISION array of DIMENSION at least
89  ( 1 + ( n - 1 )*abs( INCX ) ) when TRANS = 'N' or 'n'
90  and at least
91  ( 1 + ( m - 1 )*abs( INCX ) ) otherwise.
92  Before entry, the incremented array X must contain the
93  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 + ( m - 1 )*abs( INCY ) ) when TRANS = 'N' or 'n'
105  and at least
106  ( 1 + ( n - 1 )*abs( INCY ) ) otherwise.
107  Before entry with BETA non-zero, the incremented array Y
108  must contain the vector y. On exit, Y is overwritten by the
109  updated vector y.
110  INCY - INTEGER.
111  On entry, INCY specifies the increment for the elements of
112  Y. INCY must not be zero.
113  Unchanged on exit.
114  Level 2 Blas routine.
115  -- Written on 22-October-1986.
116  Jack Dongarra, Argonne National Lab.
117  Jeremy Du Croz, Nag Central Office.
118  Sven Hammarling, Nag Central Office.
119  Richard Hanson, Sandia National Labs.
120  Test the input parameters.
121  Parameter adjustments */
122  a_dim1 = *lda;
123  a_offset = 1 + a_dim1 * 1;
124  a -= a_offset;
125  --x;
126  --y;
127  /* Function Body */
128  info = 0;
129  if (! template_blas_lsame(trans, "N") && ! template_blas_lsame(trans, "T") && ! template_blas_lsame(trans, "C")
130  ) {
131  info = 1;
132  } else if (*m < 0) {
133  info = 2;
134  } else if (*n < 0) {
135  info = 3;
136  } else if (*lda < maxMACRO(1,*m)) {
137  info = 6;
138  } else if (*incx == 0) {
139  info = 8;
140  } else if (*incy == 0) {
141  info = 11;
142  }
143  if (info != 0) {
144  template_blas_erbla("GEMV ", &info);
145  return 0;
146  }
147 /* Quick return if possible. */
148  if (*m == 0 || *n == 0 || (*alpha == 0. && *beta == 1.) ) {
149  return 0;
150  }
151 /* Set LENX and LENY, the lengths of the vectors x and y, and set
152  up the start points in X and Y. */
153  if (template_blas_lsame(trans, "N")) {
154  lenx = *n;
155  leny = *m;
156  } else {
157  lenx = *m;
158  leny = *n;
159  }
160  if (*incx > 0) {
161  kx = 1;
162  } else {
163  kx = 1 - (lenx - 1) * *incx;
164  }
165  if (*incy > 0) {
166  ky = 1;
167  } else {
168  ky = 1 - (leny - 1) * *incy;
169  }
170 /* Start the operations. In this version the elements of A are
171  accessed sequentially with one pass through A.
172  First form y := beta*y. */
173  if (*beta != 1.) {
174  if (*incy == 1) {
175  if (*beta == 0.) {
176  i__1 = leny;
177  for (i__ = 1; i__ <= i__1; ++i__) {
178  y[i__] = 0.;
179 /* L10: */
180  }
181  } else {
182  i__1 = leny;
183  for (i__ = 1; i__ <= i__1; ++i__) {
184  y[i__] = *beta * y[i__];
185 /* L20: */
186  }
187  }
188  } else {
189  iy = ky;
190  if (*beta == 0.) {
191  i__1 = leny;
192  for (i__ = 1; i__ <= i__1; ++i__) {
193  y[iy] = 0.;
194  iy += *incy;
195 /* L30: */
196  }
197  } else {
198  i__1 = leny;
199  for (i__ = 1; i__ <= i__1; ++i__) {
200  y[iy] = *beta * y[iy];
201  iy += *incy;
202 /* L40: */
203  }
204  }
205  }
206  }
207  if (*alpha == 0.) {
208  return 0;
209  }
210  if (template_blas_lsame(trans, "N")) {
211 /* Form y := alpha*A*x + y. */
212  jx = kx;
213  if (*incy == 1) {
214  i__1 = *n;
215  for (j = 1; j <= i__1; ++j) {
216  if (x[jx] != 0.) {
217  temp = *alpha * x[jx];
218  i__2 = *m;
219  for (i__ = 1; i__ <= i__2; ++i__) {
220  y[i__] += temp * a_ref(i__, j);
221 /* L50: */
222  }
223  }
224  jx += *incx;
225 /* L60: */
226  }
227  } else {
228  i__1 = *n;
229  for (j = 1; j <= i__1; ++j) {
230  if (x[jx] != 0.) {
231  temp = *alpha * x[jx];
232  iy = ky;
233  i__2 = *m;
234  for (i__ = 1; i__ <= i__2; ++i__) {
235  y[iy] += temp * a_ref(i__, j);
236  iy += *incy;
237 /* L70: */
238  }
239  }
240  jx += *incx;
241 /* L80: */
242  }
243  }
244  } else {
245 /* Form y := alpha*A'*x + y. */
246  jy = ky;
247  if (*incx == 1) {
248  i__1 = *n;
249  for (j = 1; j <= i__1; ++j) {
250  temp = 0.;
251  i__2 = *m;
252  for (i__ = 1; i__ <= i__2; ++i__) {
253  temp += a_ref(i__, j) * x[i__];
254 /* L90: */
255  }
256  y[jy] += *alpha * temp;
257  jy += *incy;
258 /* L100: */
259  }
260  } else {
261  i__1 = *n;
262  for (j = 1; j <= i__1; ++j) {
263  temp = 0.;
264  ix = kx;
265  i__2 = *m;
266  for (i__ = 1; i__ <= i__2; ++i__) {
267  temp += a_ref(i__, j) * x[ix];
268  ix += *incx;
269 /* L110: */
270  }
271  y[jy] += *alpha * temp;
272  jy += *incy;
273 /* L120: */
274  }
275  }
276  }
277  return 0;
278 /* End of DGEMV . */
279 } /* dgemv_ */
280 #undef a_ref
281 
282 #endif
int integer
Definition: template_blas_common.h:38
#define maxMACRO(a, b)
Definition: template_blas_common.h:43
int template_blas_erbla(const char *srname, integer *info)
Definition: template_blas_common.cc:144
#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:41
logical template_blas_lsame(const char *ca, const char *cb)
Definition: template_blas_common.cc:44