ergo
template_blas_trsv.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_TRSV_HEADER
38 #define TEMPLATE_BLAS_TRSV_HEADER
39 
40 
41 template<class Treal>
42 int 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
a_ref
#define a_ref(a_1, a_2)
logical
bool logical
Definition: template_blas_common.h:41
template_blas_erbla
int template_blas_erbla(const char *srname, integer *info)
Definition: template_blas_common.cc:146
template_blas_lsame
logical template_blas_lsame(const char *ca, const char *cb)
Definition: template_blas_common.cc:46
integer
int integer
Definition: template_blas_common.h:40
template_blas_trsv
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
maxMACRO
#define maxMACRO(a, b)
Definition: template_blas_common.h:45