ergo
template_blas_syr2.h
Go to the documentation of this file.
1 /* Ergo, version 3.8, a program for linear scaling electronic structure
2  * calculations.
3  * Copyright (C) 2019 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_SYR2_HEADER
38 #define TEMPLATE_BLAS_SYR2_HEADER
39 
40 
41 template<class Treal>
42 int template_blas_syr2(const char *uplo, const integer *n, const Treal *alpha,
43  const Treal *x, const integer *incx, const Treal *y, const integer *incy,
44  Treal *a, const integer *lda)
45 {
46  /* System generated locals */
47  integer a_dim1, a_offset, i__1, i__2;
48  /* Local variables */
49  integer info;
50  Treal temp1, temp2;
51  integer i__, j;
52  integer ix, iy, jx, jy, kx, ky;
53 #define a_ref(a_1,a_2) a[(a_2)*a_dim1 + a_1]
54 /* Purpose
55  =======
56  DSYR2 performs the symmetric rank 2 operation
57  A := alpha*x*y' + alpha*y*x' + A,
58  where alpha is a scalar, x and y are n element vectors and A is an n
59  by n symmetric matrix.
60  Parameters
61  ==========
62  UPLO - CHARACTER*1.
63  On entry, UPLO specifies whether the upper or lower
64  triangular part of the array A is to be referenced as
65  follows:
66  UPLO = 'U' or 'u' Only the upper triangular part of A
67  is to be referenced.
68  UPLO = 'L' or 'l' Only the lower triangular part of A
69  is to be referenced.
70  Unchanged on exit.
71  N - INTEGER.
72  On entry, N specifies the order of the matrix A.
73  N must be at least zero.
74  Unchanged on exit.
75  ALPHA - DOUBLE PRECISION.
76  On entry, ALPHA specifies the scalar alpha.
77  Unchanged on exit.
78  X - DOUBLE PRECISION array of dimension at least
79  ( 1 + ( n - 1 )*abs( INCX ) ).
80  Before entry, the incremented array X must contain the n
81  element vector x.
82  Unchanged on exit.
83  INCX - INTEGER.
84  On entry, INCX specifies the increment for the elements of
85  X. INCX must not be zero.
86  Unchanged on exit.
87  Y - DOUBLE PRECISION array of dimension at least
88  ( 1 + ( n - 1 )*abs( INCY ) ).
89  Before entry, the incremented array Y must contain the n
90  element vector y.
91  Unchanged on exit.
92  INCY - INTEGER.
93  On entry, INCY specifies the increment for the elements of
94  Y. INCY must not be zero.
95  Unchanged on exit.
96  A - DOUBLE PRECISION array of DIMENSION ( LDA, n ).
97  Before entry with UPLO = 'U' or 'u', the leading n by n
98  upper triangular part of the array A must contain the upper
99  triangular part of the symmetric matrix and the strictly
100  lower triangular part of A is not referenced. On exit, the
101  upper triangular part of the array A is overwritten by the
102  upper triangular part of the updated matrix.
103  Before entry with UPLO = 'L' or 'l', the leading n by n
104  lower triangular part of the array A must contain the lower
105  triangular part of the symmetric matrix and the strictly
106  upper triangular part of A is not referenced. On exit, the
107  lower triangular part of the array A is overwritten by the
108  lower triangular part of the updated matrix.
109  LDA - INTEGER.
110  On entry, LDA specifies the first dimension of A as declared
111  in the calling (sub) program. LDA must be at least
112  max( 1, n ).
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  --x;
123  --y;
124  a_dim1 = *lda;
125  a_offset = 1 + a_dim1 * 1;
126  a -= a_offset;
127  /* Initialization added by Elias to get rid of compiler warnings. */
128  jx = jy = kx = ky = 0;
129  /* Function Body */
130  info = 0;
131  if (! template_blas_lsame(uplo, "U") && ! template_blas_lsame(uplo, "L")) {
132  info = 1;
133  } else if (*n < 0) {
134  info = 2;
135  } else if (*incx == 0) {
136  info = 5;
137  } else if (*incy == 0) {
138  info = 7;
139  } else if (*lda < maxMACRO(1,*n)) {
140  info = 9;
141  }
142  if (info != 0) {
143  template_blas_erbla("SYR2 ", &info);
144  return 0;
145  }
146 /* Quick return if possible. */
147  if (*n == 0 || *alpha == 0.) {
148  return 0;
149  }
150 /* Set up the start points in X and Y if the increments are not both
151  unity. */
152  if (*incx != 1 || *incy != 1) {
153  if (*incx > 0) {
154  kx = 1;
155  } else {
156  kx = 1 - (*n - 1) * *incx;
157  }
158  if (*incy > 0) {
159  ky = 1;
160  } else {
161  ky = 1 - (*n - 1) * *incy;
162  }
163  jx = kx;
164  jy = ky;
165  }
166 /* Start the operations. In this version the elements of A are
167  accessed sequentially with one pass through the triangular part
168  of A. */
169  if (template_blas_lsame(uplo, "U")) {
170 /* Form A when A is stored in the upper triangle. */
171  if (*incx == 1 && *incy == 1) {
172  i__1 = *n;
173  for (j = 1; j <= i__1; ++j) {
174  if (x[j] != 0. || y[j] != 0.) {
175  temp1 = *alpha * y[j];
176  temp2 = *alpha * x[j];
177  i__2 = j;
178  for (i__ = 1; i__ <= i__2; ++i__) {
179  a_ref(i__, j) = a_ref(i__, j) + x[i__] * temp1 + y[
180  i__] * temp2;
181 /* L10: */
182  }
183  }
184 /* L20: */
185  }
186  } else {
187  i__1 = *n;
188  for (j = 1; j <= i__1; ++j) {
189  if (x[jx] != 0. || y[jy] != 0.) {
190  temp1 = *alpha * y[jy];
191  temp2 = *alpha * x[jx];
192  ix = kx;
193  iy = ky;
194  i__2 = j;
195  for (i__ = 1; i__ <= i__2; ++i__) {
196  a_ref(i__, j) = a_ref(i__, j) + x[ix] * temp1 + y[iy]
197  * temp2;
198  ix += *incx;
199  iy += *incy;
200 /* L30: */
201  }
202  }
203  jx += *incx;
204  jy += *incy;
205 /* L40: */
206  }
207  }
208  } else {
209 /* Form A when A is stored in the lower triangle. */
210  if (*incx == 1 && *incy == 1) {
211  i__1 = *n;
212  for (j = 1; j <= i__1; ++j) {
213  if (x[j] != 0. || y[j] != 0.) {
214  temp1 = *alpha * y[j];
215  temp2 = *alpha * x[j];
216  i__2 = *n;
217  for (i__ = j; i__ <= i__2; ++i__) {
218  a_ref(i__, j) = a_ref(i__, j) + x[i__] * temp1 + y[
219  i__] * temp2;
220 /* L50: */
221  }
222  }
223 /* L60: */
224  }
225  } else {
226  i__1 = *n;
227  for (j = 1; j <= i__1; ++j) {
228  if (x[jx] != 0. || y[jy] != 0.) {
229  temp1 = *alpha * y[jy];
230  temp2 = *alpha * x[jx];
231  ix = jx;
232  iy = jy;
233  i__2 = *n;
234  for (i__ = j; i__ <= i__2; ++i__) {
235  a_ref(i__, j) = a_ref(i__, j) + x[ix] * temp1 + y[iy]
236  * temp2;
237  ix += *incx;
238  iy += *incy;
239 /* L70: */
240  }
241  }
242  jx += *incx;
243  jy += *incy;
244 /* L80: */
245  }
246  }
247  }
248  return 0;
249 /* End of DSYR2 . */
250 } /* dsyr2_ */
251 #undef a_ref
252 
253 #endif
template_blas_erbla
int template_blas_erbla(const char *srname, integer *info)
Definition: template_blas_common.cc:146
template_blas_lsame
logical template_blas_lsame(const char *ca, const char *cb)
Definition: template_blas_common.cc:46
a_ref
#define a_ref(a_1, a_2)
integer
int integer
Definition: template_blas_common.h:40
template_blas_syr2
int template_blas_syr2(const char *uplo, const integer *n, const Treal *alpha, const Treal *x, const integer *incx, const Treal *y, const integer *incy, Treal *a, const integer *lda)
Definition: template_blas_syr2.h:42
maxMACRO
#define maxMACRO(a, b)
Definition: template_blas_common.h:45