ergo
template_blas_tpsv.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_TPSV_HEADER
36 #define TEMPLATE_BLAS_TPSV_HEADER
37 
38 
39 template<class Treal>
40 int template_blas_tpsv(const char *uplo, const char *trans, const char *diag, const integer *n,
41  const Treal *ap, Treal *x, const integer *incx)
42 {
43  /* System generated locals */
44  integer i__1, i__2;
45  /* Local variables */
46  integer info;
47  Treal temp;
48  integer i__, j, k;
49  integer kk, ix, jx, kx;
50  logical nounit;
51 /* Purpose
52  =======
53  DTPSV solves one of the systems of equations
54  A*x = b, or A'*x = b,
55  where b and x are n element vectors and A is an n by n unit, or
56  non-unit, upper or lower triangular matrix, supplied in packed form.
57  No test for singularity or near-singularity is included in this
58  routine. Such tests must be performed before calling this routine.
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 equations to be solved as
69  follows:
70  TRANS = 'N' or 'n' A*x = b.
71  TRANS = 'T' or 't' A'*x = b.
72  TRANS = 'C' or 'c' A'*x = b.
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 right-hand side vector b. On exit, X is overwritten
104  with the solution 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("DTPSV ", &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 := inv( A )*x. */
156  if (template_blas_lsame(uplo, "U")) {
157  kk = *n * (*n + 1) / 2;
158  if (*incx == 1) {
159  for (j = *n; j >= 1; --j) {
160  if (x[j] != 0.) {
161  if (nounit) {
162  x[j] /= ap[kk];
163  }
164  temp = x[j];
165  k = kk - 1;
166  for (i__ = j - 1; i__ >= 1; --i__) {
167  x[i__] -= temp * ap[k];
168  --k;
169 /* L10: */
170  }
171  }
172  kk -= j;
173 /* L20: */
174  }
175  } else {
176  jx = kx + (*n - 1) * *incx;
177  for (j = *n; j >= 1; --j) {
178  if (x[jx] != 0.) {
179  if (nounit) {
180  x[jx] /= ap[kk];
181  }
182  temp = x[jx];
183  ix = jx;
184  i__1 = kk - j + 1;
185  for (k = kk - 1; k >= i__1; --k) {
186  ix -= *incx;
187  x[ix] -= temp * ap[k];
188 /* L30: */
189  }
190  }
191  jx -= *incx;
192  kk -= j;
193 /* L40: */
194  }
195  }
196  } else {
197  kk = 1;
198  if (*incx == 1) {
199  i__1 = *n;
200  for (j = 1; j <= i__1; ++j) {
201  if (x[j] != 0.) {
202  if (nounit) {
203  x[j] /= ap[kk];
204  }
205  temp = x[j];
206  k = kk + 1;
207  i__2 = *n;
208  for (i__ = j + 1; i__ <= i__2; ++i__) {
209  x[i__] -= temp * ap[k];
210  ++k;
211 /* L50: */
212  }
213  }
214  kk += *n - j + 1;
215 /* L60: */
216  }
217  } else {
218  jx = kx;
219  i__1 = *n;
220  for (j = 1; j <= i__1; ++j) {
221  if (x[jx] != 0.) {
222  if (nounit) {
223  x[jx] /= ap[kk];
224  }
225  temp = x[jx];
226  ix = jx;
227  i__2 = kk + *n - j;
228  for (k = kk + 1; k <= i__2; ++k) {
229  ix += *incx;
230  x[ix] -= temp * ap[k];
231 /* L70: */
232  }
233  }
234  jx += *incx;
235  kk += *n - j + 1;
236 /* L80: */
237  }
238  }
239  }
240  } else {
241 /* Form x := inv( A' )*x. */
242  if (template_blas_lsame(uplo, "U")) {
243  kk = 1;
244  if (*incx == 1) {
245  i__1 = *n;
246  for (j = 1; j <= i__1; ++j) {
247  temp = x[j];
248  k = kk;
249  i__2 = j - 1;
250  for (i__ = 1; i__ <= i__2; ++i__) {
251  temp -= ap[k] * x[i__];
252  ++k;
253 /* L90: */
254  }
255  if (nounit) {
256  temp /= ap[kk + j - 1];
257  }
258  x[j] = temp;
259  kk += j;
260 /* L100: */
261  }
262  } else {
263  jx = kx;
264  i__1 = *n;
265  for (j = 1; j <= i__1; ++j) {
266  temp = x[jx];
267  ix = kx;
268  i__2 = kk + j - 2;
269  for (k = kk; k <= i__2; ++k) {
270  temp -= ap[k] * x[ix];
271  ix += *incx;
272 /* L110: */
273  }
274  if (nounit) {
275  temp /= ap[kk + j - 1];
276  }
277  x[jx] = temp;
278  jx += *incx;
279  kk += j;
280 /* L120: */
281  }
282  }
283  } else {
284  kk = *n * (*n + 1) / 2;
285  if (*incx == 1) {
286  for (j = *n; j >= 1; --j) {
287  temp = x[j];
288  k = kk;
289  i__1 = j + 1;
290  for (i__ = *n; i__ >= i__1; --i__) {
291  temp -= ap[k] * x[i__];
292  --k;
293 /* L130: */
294  }
295  if (nounit) {
296  temp /= ap[kk - *n + j];
297  }
298  x[j] = temp;
299  kk -= *n - j + 1;
300 /* L140: */
301  }
302  } else {
303  kx += (*n - 1) * *incx;
304  jx = kx;
305  for (j = *n; j >= 1; --j) {
306  temp = x[jx];
307  ix = kx;
308  i__1 = kk - (*n - (j + 1));
309  for (k = kk; k >= i__1; --k) {
310  temp -= ap[k] * x[ix];
311  ix -= *incx;
312 /* L150: */
313  }
314  if (nounit) {
315  temp /= ap[kk - *n + j];
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 DTPSV . */
327 } /* dtpsv_ */
328 
329 #endif
int integer
Definition: template_blas_common.h:38
int template_blas_erbla(const char *srname, integer *info)
Definition: template_blas_common.cc:144
bool logical
Definition: template_blas_common.h:39
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:40
logical template_blas_lsame(const char *ca, const char *cb)
Definition: template_blas_common.cc:44