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