ergo
template_blas_syr2.h
Go to the documentation of this file.
1/* Ergo, version 3.8.2, a program for linear scaling electronic structure
2 * calculations.
3 * Copyright (C) 2023 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
41template<class Treal>
42int 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
int template_blas_erbla(const char *srname, integer *info)
Definition: template_blas_common.cc:146
logical template_blas_lsame(const char *ca, const char *cb)
Definition: template_blas_common.cc:46
int integer
Definition: template_blas_common.h:40
#define maxMACRO(a, b)
Definition: template_blas_common.h:45
#define a_ref(a_1, a_2)
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