ergo
template_blas_tpmv.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_TPMV_HEADER
38 #define TEMPLATE_BLAS_TPMV_HEADER
39 
40 
41 template<class Treal>
42 int template_blas_tpmv(const char *uplo, const char *trans, const char *diag, const integer *n,
43  const Treal *ap, Treal *x, const integer *incx)
44 {
45  /* System generated locals */
46  integer i__1, i__2;
47  /* Local variables */
48  integer info;
49  Treal temp;
50  integer i__, j, k;
51  integer kk, ix, jx, kx;
52  logical nounit;
53 /* Purpose
54  =======
55  DTPMV performs one of the matrix-vector operations
56  x := A*x, or x := A'*x,
57  where x is an n element vector and A is an n by n unit, or non-unit,
58  upper or lower triangular matrix, supplied in packed form.
59  Parameters
60  ==========
61  UPLO - CHARACTER*1.
62  On entry, UPLO specifies whether the matrix is an upper or
63  lower triangular matrix as follows:
64  UPLO = 'U' or 'u' A is an upper triangular matrix.
65  UPLO = 'L' or 'l' A is a lower triangular matrix.
66  Unchanged on exit.
67  TRANS - CHARACTER*1.
68  On entry, TRANS specifies the operation to be performed as
69  follows:
70  TRANS = 'N' or 'n' x := A*x.
71  TRANS = 'T' or 't' x := A'*x.
72  TRANS = 'C' or 'c' x := A'*x.
73  Unchanged on exit.
74  DIAG - CHARACTER*1.
75  On entry, DIAG specifies whether or not A is unit
76  triangular as follows:
77  DIAG = 'U' or 'u' A is assumed to be unit triangular.
78  DIAG = 'N' or 'n' A is not assumed to be unit
79  triangular.
80  Unchanged on exit.
81  N - INTEGER.
82  On entry, N specifies the order of the matrix A.
83  N must be at least zero.
84  Unchanged on exit.
85  AP - DOUBLE PRECISION array of DIMENSION at least
86  ( ( n*( n + 1 ) )/2 ).
87  Before entry with UPLO = 'U' or 'u', the array AP must
88  contain the upper triangular matrix packed sequentially,
89  column by column, so that AP( 1 ) contains a( 1, 1 ),
90  AP( 2 ) and AP( 3 ) contain a( 1, 2 ) and a( 2, 2 )
91  respectively, and so on.
92  Before entry with UPLO = 'L' or 'l', the array AP must
93  contain the lower triangular matrix packed sequentially,
94  column by column, so that AP( 1 ) contains a( 1, 1 ),
95  AP( 2 ) and AP( 3 ) contain a( 2, 1 ) and a( 3, 1 )
96  respectively, and so on.
97  Note that when DIAG = 'U' or 'u', the diagonal elements of
98  A are not referenced, but are assumed to be unity.
99  Unchanged on exit.
100  X - DOUBLE PRECISION array of dimension at least
101  ( 1 + ( n - 1 )*abs( INCX ) ).
102  Before entry, the incremented array X must contain the n
103  element vector x. On exit, X is overwritten with the
104  tranformed vector x.
105  INCX - INTEGER.
106  On entry, INCX specifies the increment for the elements of
107  X. INCX must not be zero.
108  Unchanged on exit.
109  Level 2 Blas routine.
110  -- Written on 22-October-1986.
111  Jack Dongarra, Argonne National Lab.
112  Jeremy Du Croz, Nag Central Office.
113  Sven Hammarling, Nag Central Office.
114  Richard Hanson, Sandia National Labs.
115  Test the input parameters.
116  Parameter adjustments */
117  --x;
118  --ap;
119  /* Initialization added by Elias to get rid of compiler warnings. */
120  kx = 0;
121  /* Function Body */
122  info = 0;
123  if (! template_blas_lsame(uplo, "U") && ! template_blas_lsame(uplo, "L")) {
124  info = 1;
125  } else if (! template_blas_lsame(trans, "N") && ! template_blas_lsame(trans,
126  "T") && ! template_blas_lsame(trans, "C")) {
127  info = 2;
128  } else if (! template_blas_lsame(diag, "U") && ! template_blas_lsame(diag,
129  "N")) {
130  info = 3;
131  } else if (*n < 0) {
132  info = 4;
133  } else if (*incx == 0) {
134  info = 7;
135  }
136  if (info != 0) {
137  template_blas_erbla("TPMV ", &info);
138  return 0;
139  }
140 /* Quick return if possible. */
141  if (*n == 0) {
142  return 0;
143  }
144  nounit = template_blas_lsame(diag, "N");
145 /* Set up the start point in X if the increment is not unity. This
146  will be ( N - 1 )*INCX too small for descending loops. */
147  if (*incx <= 0) {
148  kx = 1 - (*n - 1) * *incx;
149  } else if (*incx != 1) {
150  kx = 1;
151  }
152 /* Start the operations. In this version the elements of AP are
153  accessed sequentially with one pass through AP. */
154  if (template_blas_lsame(trans, "N")) {
155 /* Form x:= A*x. */
156  if (template_blas_lsame(uplo, "U")) {
157  kk = 1;
158  if (*incx == 1) {
159  i__1 = *n;
160  for (j = 1; j <= i__1; ++j) {
161  if (x[j] != 0.) {
162  temp = x[j];
163  k = kk;
164  i__2 = j - 1;
165  for (i__ = 1; i__ <= i__2; ++i__) {
166  x[i__] += temp * ap[k];
167  ++k;
168 /* L10: */
169  }
170  if (nounit) {
171  x[j] *= ap[kk + j - 1];
172  }
173  }
174  kk += j;
175 /* L20: */
176  }
177  } else {
178  jx = kx;
179  i__1 = *n;
180  for (j = 1; j <= i__1; ++j) {
181  if (x[jx] != 0.) {
182  temp = x[jx];
183  ix = kx;
184  i__2 = kk + j - 2;
185  for (k = kk; k <= i__2; ++k) {
186  x[ix] += temp * ap[k];
187  ix += *incx;
188 /* L30: */
189  }
190  if (nounit) {
191  x[jx] *= ap[kk + j - 1];
192  }
193  }
194  jx += *incx;
195  kk += j;
196 /* L40: */
197  }
198  }
199  } else {
200  kk = *n * (*n + 1) / 2;
201  if (*incx == 1) {
202  for (j = *n; j >= 1; --j) {
203  if (x[j] != 0.) {
204  temp = x[j];
205  k = kk;
206  i__1 = j + 1;
207  for (i__ = *n; i__ >= i__1; --i__) {
208  x[i__] += temp * ap[k];
209  --k;
210 /* L50: */
211  }
212  if (nounit) {
213  x[j] *= ap[kk - *n + j];
214  }
215  }
216  kk -= *n - j + 1;
217 /* L60: */
218  }
219  } else {
220  kx += (*n - 1) * *incx;
221  jx = kx;
222  for (j = *n; j >= 1; --j) {
223  if (x[jx] != 0.) {
224  temp = x[jx];
225  ix = kx;
226  i__1 = kk - (*n - (j + 1));
227  for (k = kk; k >= i__1; --k) {
228  x[ix] += temp * ap[k];
229  ix -= *incx;
230 /* L70: */
231  }
232  if (nounit) {
233  x[jx] *= ap[kk - *n + j];
234  }
235  }
236  jx -= *incx;
237  kk -= *n - j + 1;
238 /* L80: */
239  }
240  }
241  }
242  } else {
243 /* Form x := A'*x. */
244  if (template_blas_lsame(uplo, "U")) {
245  kk = *n * (*n + 1) / 2;
246  if (*incx == 1) {
247  for (j = *n; j >= 1; --j) {
248  temp = x[j];
249  if (nounit) {
250  temp *= ap[kk];
251  }
252  k = kk - 1;
253  for (i__ = j - 1; i__ >= 1; --i__) {
254  temp += ap[k] * x[i__];
255  --k;
256 /* L90: */
257  }
258  x[j] = temp;
259  kk -= j;
260 /* L100: */
261  }
262  } else {
263  jx = kx + (*n - 1) * *incx;
264  for (j = *n; j >= 1; --j) {
265  temp = x[jx];
266  ix = jx;
267  if (nounit) {
268  temp *= ap[kk];
269  }
270  i__1 = kk - j + 1;
271  for (k = kk - 1; k >= i__1; --k) {
272  ix -= *incx;
273  temp += ap[k] * x[ix];
274 /* L110: */
275  }
276  x[jx] = temp;
277  jx -= *incx;
278  kk -= j;
279 /* L120: */
280  }
281  }
282  } else {
283  kk = 1;
284  if (*incx == 1) {
285  i__1 = *n;
286  for (j = 1; j <= i__1; ++j) {
287  temp = x[j];
288  if (nounit) {
289  temp *= ap[kk];
290  }
291  k = kk + 1;
292  i__2 = *n;
293  for (i__ = j + 1; i__ <= i__2; ++i__) {
294  temp += ap[k] * x[i__];
295  ++k;
296 /* L130: */
297  }
298  x[j] = temp;
299  kk += *n - j + 1;
300 /* L140: */
301  }
302  } else {
303  jx = kx;
304  i__1 = *n;
305  for (j = 1; j <= i__1; ++j) {
306  temp = x[jx];
307  ix = jx;
308  if (nounit) {
309  temp *= ap[kk];
310  }
311  i__2 = kk + *n - j;
312  for (k = kk + 1; k <= i__2; ++k) {
313  ix += *incx;
314  temp += ap[k] * x[ix];
315 /* L150: */
316  }
317  x[jx] = temp;
318  jx += *incx;
319  kk += *n - j + 1;
320 /* L160: */
321  }
322  }
323  }
324  }
325  return 0;
326 /* End of DTPMV . */
327 } /* dtpmv_ */
328 
329 #endif
int template_blas_tpmv(const char *uplo, const char *trans, const char *diag, const integer *n, const Treal *ap, Treal *x, const integer *incx)
Definition: template_blas_tpmv.h:42
int integer
Definition: template_blas_common.h:40
int template_blas_erbla(const char *srname, integer *info)
Definition: template_blas_common.cc:146
bool logical
Definition: template_blas_common.h:41
logical template_blas_lsame(const char *ca, const char *cb)
Definition: template_blas_common.cc:46