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