ergo
template_blas_spr2.h
Go to the documentation of this file.
1 /* Ergo, version 3.7, a program for linear scaling electronic structure
2  * calculations.
3  * Copyright (C) 2018 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 
40 #include "template_blas_common.h"
41 
42 template<class Treal>
43 int 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 integer
Definition: template_blas_common.h:40
int template_blas_erbla(const char *srname, integer *info)
Definition: template_blas_common.cc:146
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
logical template_blas_lsame(const char *ca, const char *cb)
Definition: template_blas_common.cc:46