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