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