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