ergo
template_blas_tpsv.h
Go to the documentation of this file.
1/* Ergo, version 3.8, a program for linear scaling electronic structure
2 * calculations.
3 * Copyright (C) 2019 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
41template<class Treal>
42int 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 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
int integer
Definition: template_blas_common.h:40
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