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