ergo
template_blas_syr2.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_SYR2_HEADER
36 #define TEMPLATE_BLAS_SYR2_HEADER
37 
38 
39 template<class Treal>
40 int template_blas_syr2(const char *uplo, const integer *n, const Treal *alpha,
41  const Treal *x, const integer *incx, const Treal *y, const integer *incy,
42  Treal *a, const integer *lda)
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  DSYR2 performs the symmetric rank 2 operation
55  A := alpha*x*y' + alpha*y*x' + A,
56  where alpha is a scalar, x and y are n element vectors and A is an n
57  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  X - DOUBLE PRECISION array of dimension at least
77  ( 1 + ( n - 1 )*abs( INCX ) ).
78  Before entry, the incremented array X must contain the n
79  element vector x.
80  Unchanged on exit.
81  INCX - INTEGER.
82  On entry, INCX specifies the increment for the elements of
83  X. INCX must not be zero.
84  Unchanged on exit.
85  Y - DOUBLE PRECISION array of dimension at least
86  ( 1 + ( n - 1 )*abs( INCY ) ).
87  Before entry, the incremented array Y must contain the n
88  element vector y.
89  Unchanged on exit.
90  INCY - INTEGER.
91  On entry, INCY specifies the increment for the elements of
92  Y. INCY must not be zero.
93  Unchanged on exit.
94  A - DOUBLE PRECISION array of DIMENSION ( LDA, n ).
95  Before entry with UPLO = 'U' or 'u', the leading n by n
96  upper triangular part of the array A must contain the upper
97  triangular part of the symmetric matrix and the strictly
98  lower triangular part of A is not referenced. On exit, the
99  upper triangular part of the array A is overwritten by the
100  upper triangular part of the updated matrix.
101  Before entry with UPLO = 'L' or 'l', the leading n by n
102  lower triangular part of the array A must contain the lower
103  triangular part of the symmetric matrix and the strictly
104  upper triangular part of A is not referenced. On exit, the
105  lower triangular part of the array A is overwritten by the
106  lower triangular part of the updated matrix.
107  LDA - INTEGER.
108  On entry, LDA specifies the first dimension of A as declared
109  in the calling (sub) program. LDA must be at least
110  max( 1, n ).
111  Unchanged on exit.
112  Level 2 Blas routine.
113  -- Written on 22-October-1986.
114  Jack Dongarra, Argonne National Lab.
115  Jeremy Du Croz, Nag Central Office.
116  Sven Hammarling, Nag Central Office.
117  Richard Hanson, Sandia National Labs.
118  Test the input parameters.
119  Parameter adjustments */
120  --x;
121  --y;
122  a_dim1 = *lda;
123  a_offset = 1 + a_dim1 * 1;
124  a -= a_offset;
125  /* Initialization added by Elias to get rid of compiler warnings. */
126  jx = jy = kx = ky = 0;
127  /* Function Body */
128  info = 0;
129  if (! template_blas_lsame(uplo, "U") && ! template_blas_lsame(uplo, "L")) {
130  info = 1;
131  } else if (*n < 0) {
132  info = 2;
133  } else if (*incx == 0) {
134  info = 5;
135  } else if (*incy == 0) {
136  info = 7;
137  } else if (*lda < maxMACRO(1,*n)) {
138  info = 9;
139  }
140  if (info != 0) {
141  template_blas_erbla("SYR2 ", &info);
142  return 0;
143  }
144 /* Quick return if possible. */
145  if (*n == 0 || *alpha == 0.) {
146  return 0;
147  }
148 /* Set up the start points in X and Y if the increments are not both
149  unity. */
150  if (*incx != 1 || *incy != 1) {
151  if (*incx > 0) {
152  kx = 1;
153  } else {
154  kx = 1 - (*n - 1) * *incx;
155  }
156  if (*incy > 0) {
157  ky = 1;
158  } else {
159  ky = 1 - (*n - 1) * *incy;
160  }
161  jx = kx;
162  jy = ky;
163  }
164 /* Start the operations. In this version the elements of A are
165  accessed sequentially with one pass through the triangular part
166  of A. */
167  if (template_blas_lsame(uplo, "U")) {
168 /* Form A when A is stored in the upper triangle. */
169  if (*incx == 1 && *incy == 1) {
170  i__1 = *n;
171  for (j = 1; j <= i__1; ++j) {
172  if (x[j] != 0. || y[j] != 0.) {
173  temp1 = *alpha * y[j];
174  temp2 = *alpha * x[j];
175  i__2 = j;
176  for (i__ = 1; i__ <= i__2; ++i__) {
177  a_ref(i__, j) = a_ref(i__, j) + x[i__] * temp1 + y[
178  i__] * temp2;
179 /* L10: */
180  }
181  }
182 /* L20: */
183  }
184  } else {
185  i__1 = *n;
186  for (j = 1; j <= i__1; ++j) {
187  if (x[jx] != 0. || y[jy] != 0.) {
188  temp1 = *alpha * y[jy];
189  temp2 = *alpha * x[jx];
190  ix = kx;
191  iy = ky;
192  i__2 = j;
193  for (i__ = 1; i__ <= i__2; ++i__) {
194  a_ref(i__, j) = a_ref(i__, j) + x[ix] * temp1 + y[iy]
195  * temp2;
196  ix += *incx;
197  iy += *incy;
198 /* L30: */
199  }
200  }
201  jx += *incx;
202  jy += *incy;
203 /* L40: */
204  }
205  }
206  } else {
207 /* Form A when A is stored in the lower triangle. */
208  if (*incx == 1 && *incy == 1) {
209  i__1 = *n;
210  for (j = 1; j <= i__1; ++j) {
211  if (x[j] != 0. || y[j] != 0.) {
212  temp1 = *alpha * y[j];
213  temp2 = *alpha * x[j];
214  i__2 = *n;
215  for (i__ = j; i__ <= i__2; ++i__) {
216  a_ref(i__, j) = a_ref(i__, j) + x[i__] * temp1 + y[
217  i__] * temp2;
218 /* L50: */
219  }
220  }
221 /* L60: */
222  }
223  } else {
224  i__1 = *n;
225  for (j = 1; j <= i__1; ++j) {
226  if (x[jx] != 0. || y[jy] != 0.) {
227  temp1 = *alpha * y[jy];
228  temp2 = *alpha * x[jx];
229  ix = jx;
230  iy = jy;
231  i__2 = *n;
232  for (i__ = j; i__ <= i__2; ++i__) {
233  a_ref(i__, j) = a_ref(i__, j) + x[ix] * temp1 + y[iy]
234  * temp2;
235  ix += *incx;
236  iy += *incy;
237 /* L70: */
238  }
239  }
240  jx += *incx;
241  jy += *incy;
242 /* L80: */
243  }
244  }
245  }
246  return 0;
247 /* End of DSYR2 . */
248 } /* dsyr2_ */
249 #undef a_ref
250 
251 #endif