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