ergo
template_blas_spr2.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_SPR2_HEADER
38#define TEMPLATE_BLAS_SPR2_HEADER
39
41
42template<class Treal>
43int template_blas_spr2(const char *uplo, const integer *n, const Treal *alpha,
44 const Treal *x, const integer *incx, const Treal *y, const integer *incy,
45 Treal *ap)
46{
47 /* System generated locals */
48 integer i__1, i__2;
49 /* Local variables */
50 integer info;
51 Treal temp1, temp2;
52 integer i__, j, k;
53 integer kk, ix, iy, jx, jy, kx, ky;
54/* Purpose
55 =======
56 DSPR2 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
59 n by n symmetric matrix, supplied in packed form.
60 Parameters
61 ==========
62 UPLO - CHARACTER*1.
63 On entry, UPLO specifies whether the upper or lower
64 triangular part of the matrix A is supplied in the packed
65 array AP as follows:
66 UPLO = 'U' or 'u' The upper triangular part of A is
67 supplied in AP.
68 UPLO = 'L' or 'l' The lower triangular part of A is
69 supplied in AP.
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 AP - DOUBLE PRECISION array of DIMENSION at least
97 ( ( n*( n + 1 ) )/2 ).
98 Before entry with UPLO = 'U' or 'u', the array AP must
99 contain the upper triangular part of the symmetric matrix
100 packed sequentially, column by column, so that AP( 1 )
101 contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 1, 2 )
102 and a( 2, 2 ) respectively, and so on. On exit, the array
103 AP is overwritten by the upper triangular part of the
104 updated matrix.
105 Before entry with UPLO = 'L' or 'l', the array AP must
106 contain the lower triangular part of the symmetric matrix
107 packed sequentially, column by column, so that AP( 1 )
108 contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 2, 1 )
109 and a( 3, 1 ) respectively, and so on. On exit, the array
110 AP is overwritten by the lower triangular part of the
111 updated matrix.
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 --ap;
121 --y;
122 --x;
123 /* Initialization added by Elias to get rid of compiler warnings. */
124 jx = jy = kx = ky = 0;
125 /* Function Body */
126 info = 0;
127 if (! template_blas_lsame(uplo, "U") && ! template_blas_lsame(uplo, "L")) {
128 info = 1;
129 } else if (*n < 0) {
130 info = 2;
131 } else if (*incx == 0) {
132 info = 5;
133 } else if (*incy == 0) {
134 info = 7;
135 }
136 if (info != 0) {
137 template_blas_erbla("SPR2 ", &info);
138 return 0;
139 }
140/* Quick return if possible. */
141 if (*n == 0 || *alpha == 0.) {
142 return 0;
143 }
144/* Set up the start points in X and Y if the increments are not both
145 unity. */
146 if (*incx != 1 || *incy != 1) {
147 if (*incx > 0) {
148 kx = 1;
149 } else {
150 kx = 1 - (*n - 1) * *incx;
151 }
152 if (*incy > 0) {
153 ky = 1;
154 } else {
155 ky = 1 - (*n - 1) * *incy;
156 }
157 jx = kx;
158 jy = ky;
159 }
160/* Start the operations. In this version the elements of the array AP
161 are accessed sequentially with one pass through AP. */
162 kk = 1;
163 if (template_blas_lsame(uplo, "U")) {
164/* Form A when upper triangle is stored in AP. */
165 if (*incx == 1 && *incy == 1) {
166 i__1 = *n;
167 for (j = 1; j <= i__1; ++j) {
168 if (x[j] != 0. || y[j] != 0.) {
169 temp1 = *alpha * y[j];
170 temp2 = *alpha * x[j];
171 k = kk;
172 i__2 = j;
173 for (i__ = 1; i__ <= i__2; ++i__) {
174 ap[k] = ap[k] + x[i__] * temp1 + y[i__] * temp2;
175 ++k;
176/* L10: */
177 }
178 }
179 kk += j;
180/* L20: */
181 }
182 } else {
183 i__1 = *n;
184 for (j = 1; j <= i__1; ++j) {
185 if (x[jx] != 0. || y[jy] != 0.) {
186 temp1 = *alpha * y[jy];
187 temp2 = *alpha * x[jx];
188 ix = kx;
189 iy = ky;
190 i__2 = kk + j - 1;
191 for (k = kk; k <= i__2; ++k) {
192 ap[k] = ap[k] + x[ix] * temp1 + y[iy] * temp2;
193 ix += *incx;
194 iy += *incy;
195/* L30: */
196 }
197 }
198 jx += *incx;
199 jy += *incy;
200 kk += j;
201/* L40: */
202 }
203 }
204 } else {
205/* Form A when lower triangle is stored in AP. */
206 if (*incx == 1 && *incy == 1) {
207 i__1 = *n;
208 for (j = 1; j <= i__1; ++j) {
209 if (x[j] != 0. || y[j] != 0.) {
210 temp1 = *alpha * y[j];
211 temp2 = *alpha * x[j];
212 k = kk;
213 i__2 = *n;
214 for (i__ = j; i__ <= i__2; ++i__) {
215 ap[k] = ap[k] + x[i__] * temp1 + y[i__] * temp2;
216 ++k;
217/* L50: */
218 }
219 }
220 kk = kk + *n - j + 1;
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 = kk + *n - j;
232 for (k = kk; k <= i__2; ++k) {
233 ap[k] = ap[k] + x[ix] * temp1 + y[iy] * temp2;
234 ix += *incx;
235 iy += *incy;
236/* L70: */
237 }
238 }
239 jx += *incx;
240 jy += *incy;
241 kk = kk + *n - j + 1;
242/* L80: */
243 }
244 }
245 }
246 return 0;
247/* End of DSPR2 . */
248} /* dspr2_ */
249
250#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
int template_blas_spr2(const char *uplo, const integer *n, const Treal *alpha, const Treal *x, const integer *incx, const Treal *y, const integer *incy, Treal *ap)
Definition: template_blas_spr2.h:43