ergo
template_blas_trmv.h
Go to the documentation of this file.
1 /* Ergo, version 3.7, a program for linear scaling electronic structure
2  * calculations.
3  * Copyright (C) 2018 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 
41 template<class Treal>
42 int 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
#define a_ref(a_1, a_2)
int integer
Definition: template_blas_common.h:40
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
#define maxMACRO(a, b)
Definition: template_blas_common.h:45
int template_blas_erbla(const char *srname, integer *info)
Definition: template_blas_common.cc:146
bool logical
Definition: template_blas_common.h:41
logical template_blas_lsame(const char *ca, const char *cb)
Definition: template_blas_common.cc:46