ergo
template_blas_spr.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_SPR_HEADER
36 #define TEMPLATE_BLAS_SPR_HEADER
37 
38 #include "template_blas_common.h"
39 
40 template<class Treal>
41 int template_blas_spr(const char *uplo, integer *n, Treal *alpha,
42  Treal *x, integer *incx, Treal *ap)
43 {
44  /* System generated locals */
45  integer i__1, i__2;
46  /* Local variables */
47  integer info;
48  Treal temp;
49  integer i__, j, k;
50  integer kk, ix, jx, kx;
51 /* Purpose
52  =======
53  DSPR performs the symmetric rank 1 operation
54  A := alpha*x*x' + A,
55  where alpha is a real scalar, x is an n element vector and A is an
56  n by n symmetric matrix, supplied in packed form.
57  Parameters
58  ==========
59  UPLO - CHARACTER*1.
60  On entry, UPLO specifies whether the upper or lower
61  triangular part of the matrix A is supplied in the packed
62  array AP as follows:
63  UPLO = 'U' or 'u' The upper triangular part of A is
64  supplied in AP.
65  UPLO = 'L' or 'l' The lower triangular part of A is
66  supplied in AP.
67  Unchanged on exit.
68  N - INTEGER.
69  On entry, N specifies the order of the matrix A.
70  N must be at least zero.
71  Unchanged on exit.
72  ALPHA - DOUBLE PRECISION.
73  On entry, ALPHA specifies the scalar alpha.
74  Unchanged on exit.
75  X - DOUBLE PRECISION array of dimension at least
76  ( 1 + ( n - 1 )*abs( INCX ) ).
77  Before entry, the incremented array X must contain the n
78  element vector x.
79  Unchanged on exit.
80  INCX - INTEGER.
81  On entry, INCX specifies the increment for the elements of
82  X. INCX must not be zero.
83  Unchanged on exit.
84  AP - DOUBLE PRECISION array of DIMENSION at least
85  ( ( n*( n + 1 ) )/2 ).
86  Before entry with UPLO = 'U' or 'u', the array AP must
87  contain the upper triangular part of the symmetric matrix
88  packed sequentially, column by column, so that AP( 1 )
89  contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 1, 2 )
90  and a( 2, 2 ) respectively, and so on. On exit, the array
91  AP is overwritten by the upper triangular part of the
92  updated matrix.
93  Before entry with UPLO = 'L' or 'l', the array AP must
94  contain the lower triangular part of the symmetric matrix
95  packed sequentially, column by column, so that AP( 1 )
96  contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 2, 1 )
97  and a( 3, 1 ) respectively, and so on. On exit, the array
98  AP is overwritten by the lower triangular part of the
99  updated matrix.
100  Level 2 Blas routine.
101  -- Written on 22-October-1986.
102  Jack Dongarra, Argonne National Lab.
103  Jeremy Du Croz, Nag Central Office.
104  Sven Hammarling, Nag Central Office.
105  Richard Hanson, Sandia National Labs.
106  Test the input parameters.
107  Parameter adjustments */
108  --ap;
109  --x;
110  /* Initialization added by Elias to get rid of compiler warnings. */
111  kx = 0;
112  /* Function Body */
113  info = 0;
114  if (! template_blas_lsame(uplo, "U") && ! template_blas_lsame(uplo, "L")) {
115  info = 1;
116  } else if (*n < 0) {
117  info = 2;
118  } else if (*incx == 0) {
119  info = 5;
120  }
121  if (info != 0) {
122  template_blas_erbla("SPR ", &info);
123  return 0;
124  }
125 /* Quick return if possible. */
126  if (*n == 0 || *alpha == 0.) {
127  return 0;
128  }
129 /* Set the start point in X if the increment is not unity. */
130  if (*incx <= 0) {
131  kx = 1 - (*n - 1) * *incx;
132  } else if (*incx != 1) {
133  kx = 1;
134  }
135 /* Start the operations. In this version the elements of the array AP
136  are accessed sequentially with one pass through AP. */
137  kk = 1;
138  if (template_blas_lsame(uplo, "U")) {
139 /* Form A when upper triangle is stored in AP. */
140  if (*incx == 1) {
141  i__1 = *n;
142  for (j = 1; j <= i__1; ++j) {
143  if (x[j] != 0.) {
144  temp = *alpha * x[j];
145  k = kk;
146  i__2 = j;
147  for (i__ = 1; i__ <= i__2; ++i__) {
148  ap[k] += x[i__] * temp;
149  ++k;
150 /* L10: */
151  }
152  }
153  kk += j;
154 /* L20: */
155  }
156  } else {
157  jx = kx;
158  i__1 = *n;
159  for (j = 1; j <= i__1; ++j) {
160  if (x[jx] != 0.) {
161  temp = *alpha * x[jx];
162  ix = kx;
163  i__2 = kk + j - 1;
164  for (k = kk; k <= i__2; ++k) {
165  ap[k] += x[ix] * temp;
166  ix += *incx;
167 /* L30: */
168  }
169  }
170  jx += *incx;
171  kk += j;
172 /* L40: */
173  }
174  }
175  } else {
176 /* Form A when lower triangle is stored in AP. */
177  if (*incx == 1) {
178  i__1 = *n;
179  for (j = 1; j <= i__1; ++j) {
180  if (x[j] != 0.) {
181  temp = *alpha * x[j];
182  k = kk;
183  i__2 = *n;
184  for (i__ = j; i__ <= i__2; ++i__) {
185  ap[k] += x[i__] * temp;
186  ++k;
187 /* L50: */
188  }
189  }
190  kk = kk + *n - j + 1;
191 /* L60: */
192  }
193  } else {
194  jx = kx;
195  i__1 = *n;
196  for (j = 1; j <= i__1; ++j) {
197  if (x[jx] != 0.) {
198  temp = *alpha * x[jx];
199  ix = jx;
200  i__2 = kk + *n - j;
201  for (k = kk; k <= i__2; ++k) {
202  ap[k] += x[ix] * temp;
203  ix += *incx;
204 /* L70: */
205  }
206  }
207  jx += *incx;
208  kk = kk + *n - j + 1;
209 /* L80: */
210  }
211  }
212  }
213  return 0;
214 /* End of DSPR . */
215 } /* dspr_ */
216 
217 #endif
int template_blas_spr(const char *uplo, integer *n, Treal *alpha, Treal *x, integer *incx, Treal *ap)
Definition: template_blas_spr.h:41
int integer
Definition: template_blas_common.h:38
int template_blas_erbla(const char *srname, integer *info)
Definition: template_blas_common.cc:144
logical template_blas_lsame(const char *ca, const char *cb)
Definition: template_blas_common.cc:44