ergo
template_blas_spr2.h
Go to the documentation of this file.
1 /* Ergo, version 3.3, a program for linear scaling electronic structure
2  * calculations.
3  * Copyright (C) 2013 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_SPR2_HEADER
36 #define TEMPLATE_BLAS_SPR2_HEADER
37 
38 #include "template_blas_common.h"
39 
40 template<class Treal>
41 int template_blas_spr2(const char *uplo, const integer *n, const Treal *alpha,
42  const Treal *x, const integer *incx, const Treal *y, const integer *incy,
43  Treal *ap)
44 {
45  /* System generated locals */
46  integer i__1, i__2;
47  /* Local variables */
48  integer info;
49  Treal temp1, temp2;
50  integer i__, j, k;
51  integer kk, ix, iy, jx, jy, kx, ky;
52 /* Purpose
53  =======
54  DSPR2 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
57  n by n symmetric matrix, supplied in packed form.
58  Parameters
59  ==========
60  UPLO - CHARACTER*1.
61  On entry, UPLO specifies whether the upper or lower
62  triangular part of the matrix A is supplied in the packed
63  array AP as follows:
64  UPLO = 'U' or 'u' The upper triangular part of A is
65  supplied in AP.
66  UPLO = 'L' or 'l' The lower triangular part of A is
67  supplied in AP.
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  AP - DOUBLE PRECISION array of DIMENSION at least
95  ( ( n*( n + 1 ) )/2 ).
96  Before entry with UPLO = 'U' or 'u', the array AP must
97  contain the upper triangular part of the symmetric matrix
98  packed sequentially, column by column, so that AP( 1 )
99  contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 1, 2 )
100  and a( 2, 2 ) respectively, and so on. On exit, the array
101  AP is overwritten by the upper triangular part of the
102  updated matrix.
103  Before entry with UPLO = 'L' or 'l', the array AP must
104  contain the lower triangular part of the symmetric matrix
105  packed sequentially, column by column, so that AP( 1 )
106  contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 2, 1 )
107  and a( 3, 1 ) respectively, and so on. On exit, the array
108  AP is overwritten by the lower triangular part of the
109  updated matrix.
110  Level 2 Blas routine.
111  -- Written on 22-October-1986.
112  Jack Dongarra, Argonne National Lab.
113  Jeremy Du Croz, Nag Central Office.
114  Sven Hammarling, Nag Central Office.
115  Richard Hanson, Sandia National Labs.
116  Test the input parameters.
117  Parameter adjustments */
118  --ap;
119  --y;
120  --x;
121  /* Initialization added by Elias to get rid of compiler warnings. */
122  jx = jy = kx = ky = 0;
123  /* Function Body */
124  info = 0;
125  if (! template_blas_lsame(uplo, "U") && ! template_blas_lsame(uplo, "L")) {
126  info = 1;
127  } else if (*n < 0) {
128  info = 2;
129  } else if (*incx == 0) {
130  info = 5;
131  } else if (*incy == 0) {
132  info = 7;
133  }
134  if (info != 0) {
135  template_blas_erbla("SPR2 ", &info);
136  return 0;
137  }
138 /* Quick return if possible. */
139  if (*n == 0 || *alpha == 0.) {
140  return 0;
141  }
142 /* Set up the start points in X and Y if the increments are not both
143  unity. */
144  if (*incx != 1 || *incy != 1) {
145  if (*incx > 0) {
146  kx = 1;
147  } else {
148  kx = 1 - (*n - 1) * *incx;
149  }
150  if (*incy > 0) {
151  ky = 1;
152  } else {
153  ky = 1 - (*n - 1) * *incy;
154  }
155  jx = kx;
156  jy = ky;
157  }
158 /* Start the operations. In this version the elements of the array AP
159  are accessed sequentially with one pass through AP. */
160  kk = 1;
161  if (template_blas_lsame(uplo, "U")) {
162 /* Form A when upper triangle is stored in AP. */
163  if (*incx == 1 && *incy == 1) {
164  i__1 = *n;
165  for (j = 1; j <= i__1; ++j) {
166  if (x[j] != 0. || y[j] != 0.) {
167  temp1 = *alpha * y[j];
168  temp2 = *alpha * x[j];
169  k = kk;
170  i__2 = j;
171  for (i__ = 1; i__ <= i__2; ++i__) {
172  ap[k] = ap[k] + x[i__] * temp1 + y[i__] * temp2;
173  ++k;
174 /* L10: */
175  }
176  }
177  kk += j;
178 /* L20: */
179  }
180  } else {
181  i__1 = *n;
182  for (j = 1; j <= i__1; ++j) {
183  if (x[jx] != 0. || y[jy] != 0.) {
184  temp1 = *alpha * y[jy];
185  temp2 = *alpha * x[jx];
186  ix = kx;
187  iy = ky;
188  i__2 = kk + j - 1;
189  for (k = kk; k <= i__2; ++k) {
190  ap[k] = ap[k] + x[ix] * temp1 + y[iy] * temp2;
191  ix += *incx;
192  iy += *incy;
193 /* L30: */
194  }
195  }
196  jx += *incx;
197  jy += *incy;
198  kk += j;
199 /* L40: */
200  }
201  }
202  } else {
203 /* Form A when lower triangle is stored in AP. */
204  if (*incx == 1 && *incy == 1) {
205  i__1 = *n;
206  for (j = 1; j <= i__1; ++j) {
207  if (x[j] != 0. || y[j] != 0.) {
208  temp1 = *alpha * y[j];
209  temp2 = *alpha * x[j];
210  k = kk;
211  i__2 = *n;
212  for (i__ = j; i__ <= i__2; ++i__) {
213  ap[k] = ap[k] + x[i__] * temp1 + y[i__] * temp2;
214  ++k;
215 /* L50: */
216  }
217  }
218  kk = kk + *n - j + 1;
219 /* L60: */
220  }
221  } else {
222  i__1 = *n;
223  for (j = 1; j <= i__1; ++j) {
224  if (x[jx] != 0. || y[jy] != 0.) {
225  temp1 = *alpha * y[jy];
226  temp2 = *alpha * x[jx];
227  ix = jx;
228  iy = jy;
229  i__2 = kk + *n - j;
230  for (k = kk; k <= i__2; ++k) {
231  ap[k] = ap[k] + x[ix] * temp1 + y[iy] * temp2;
232  ix += *incx;
233  iy += *incy;
234 /* L70: */
235  }
236  }
237  jx += *incx;
238  jy += *incy;
239  kk = kk + *n - j + 1;
240 /* L80: */
241  }
242  }
243  }
244  return 0;
245 /* End of DSPR2 . */
246 } /* dspr2_ */
247 
248 #endif
int integer
Definition: template_blas_common.h:38
int template_blas_erbla(const char *srname, integer *info)
Definition: template_blas_common.cc:144
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:41
logical template_blas_lsame(const char *ca, const char *cb)
Definition: template_blas_common.cc:44