ergo
template_blas_trsv.h
Go to the documentation of this file.
1 /* Ergo, version 3.4, a program for linear scaling electronic structure
2  * calculations.
3  * Copyright (C) 2014 Elias Rudberg, Emanuel H. Rubensson, and Pawel Salek.
4  *
5  * This program is free software: you can redistribute it and/or modify
6  * it under the terms of the GNU General Public License as published by
7  * the Free Software Foundation, either version 3 of the License, or
8  * (at your option) any later version.
9  *
10  * This program is distributed in the hope that it will be useful,
11  * but WITHOUT ANY WARRANTY; without even the implied warranty of
12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13  * GNU General Public License for more details.
14  *
15  * You should have received a copy of the GNU General Public License
16  * along with this program. If not, see <http://www.gnu.org/licenses/>.
17  *
18  * Primary academic reference:
19  * Kohn−Sham Density Functional Theory Electronic Structure Calculations
20  * with Linearly Scaling Computational Time and Memory Usage,
21  * Elias Rudberg, Emanuel H. Rubensson, and Pawel Salek,
22  * J. Chem. Theory Comput. 7, 340 (2011),
23  * <http://dx.doi.org/10.1021/ct100611z>
24  *
25  * For further information about Ergo, see <http://www.ergoscf.org>.
26  */
27 
28  /* This file belongs to the template_lapack part of the Ergo source
29  * code. The source files in the template_lapack directory are modified
30  * versions of files originally distributed as CLAPACK, see the
31  * Copyright/license notice in the file template_lapack/COPYING.
32  */
33 
34 
35 #ifndef TEMPLATE_BLAS_TRSV_HEADER
36 #define TEMPLATE_BLAS_TRSV_HEADER
37 
38 
39 template<class Treal>
40 int template_blas_trsv(const char *uplo, const char *trans, const char *diag, const integer *n,
41  const Treal *a, const integer *lda, Treal *x, const integer *incx)
42 {
43  /* System generated locals */
44  integer a_dim1, a_offset, i__1, i__2;
45  /* Local variables */
46  integer info;
47  Treal temp;
48  integer i__, j;
49  integer ix, jx, kx;
50  logical nounit;
51 #define a_ref(a_1,a_2) a[(a_2)*a_dim1 + a_1]
52 /* Purpose
53  =======
54  DTRSV solves one of the systems of equations
55  A*x = b, or A'*x = b,
56  where b and x are n element vectors and A is an n by n unit, or
57  non-unit, upper or lower triangular matrix.
58  No test for singularity or near-singularity is included in this
59  routine. Such tests must be performed before calling this routine.
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 equations to be solved as
70  follows:
71  TRANS = 'N' or 'n' A*x = b.
72  TRANS = 'T' or 't' A'*x = b.
73  TRANS = 'C' or 'c' A'*x = b.
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 right-hand side vector b. On exit, X is overwritten
107  with the solution 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("TRSV ", &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 := inv( A )*x. */
163  if (template_blas_lsame(uplo, "U")) {
164  if (*incx == 1) {
165  for (j = *n; j >= 1; --j) {
166  if (x[j] != 0.) {
167  if (nounit) {
168  x[j] /= a_ref(j, j);
169  }
170  temp = x[j];
171  for (i__ = j - 1; i__ >= 1; --i__) {
172  x[i__] -= temp * a_ref(i__, j);
173 /* L10: */
174  }
175  }
176 /* L20: */
177  }
178  } else {
179  jx = kx + (*n - 1) * *incx;
180  for (j = *n; j >= 1; --j) {
181  if (x[jx] != 0.) {
182  if (nounit) {
183  x[jx] /= a_ref(j, j);
184  }
185  temp = x[jx];
186  ix = jx;
187  for (i__ = j - 1; i__ >= 1; --i__) {
188  ix -= *incx;
189  x[ix] -= temp * a_ref(i__, j);
190 /* L30: */
191  }
192  }
193  jx -= *incx;
194 /* L40: */
195  }
196  }
197  } else {
198  if (*incx == 1) {
199  i__1 = *n;
200  for (j = 1; j <= i__1; ++j) {
201  if (x[j] != 0.) {
202  if (nounit) {
203  x[j] /= a_ref(j, j);
204  }
205  temp = x[j];
206  i__2 = *n;
207  for (i__ = j + 1; i__ <= i__2; ++i__) {
208  x[i__] -= temp * a_ref(i__, j);
209 /* L50: */
210  }
211  }
212 /* L60: */
213  }
214  } else {
215  jx = kx;
216  i__1 = *n;
217  for (j = 1; j <= i__1; ++j) {
218  if (x[jx] != 0.) {
219  if (nounit) {
220  x[jx] /= a_ref(j, j);
221  }
222  temp = x[jx];
223  ix = jx;
224  i__2 = *n;
225  for (i__ = j + 1; i__ <= i__2; ++i__) {
226  ix += *incx;
227  x[ix] -= temp * a_ref(i__, j);
228 /* L70: */
229  }
230  }
231  jx += *incx;
232 /* L80: */
233  }
234  }
235  }
236  } else {
237 /* Form x := inv( A' )*x. */
238  if (template_blas_lsame(uplo, "U")) {
239  if (*incx == 1) {
240  i__1 = *n;
241  for (j = 1; j <= i__1; ++j) {
242  temp = x[j];
243  i__2 = j - 1;
244  for (i__ = 1; i__ <= i__2; ++i__) {
245  temp -= a_ref(i__, j) * x[i__];
246 /* L90: */
247  }
248  if (nounit) {
249  temp /= a_ref(j, j);
250  }
251  x[j] = temp;
252 /* L100: */
253  }
254  } else {
255  jx = kx;
256  i__1 = *n;
257  for (j = 1; j <= i__1; ++j) {
258  temp = x[jx];
259  ix = kx;
260  i__2 = j - 1;
261  for (i__ = 1; i__ <= i__2; ++i__) {
262  temp -= a_ref(i__, j) * x[ix];
263  ix += *incx;
264 /* L110: */
265  }
266  if (nounit) {
267  temp /= a_ref(j, j);
268  }
269  x[jx] = temp;
270  jx += *incx;
271 /* L120: */
272  }
273  }
274  } else {
275  if (*incx == 1) {
276  for (j = *n; j >= 1; --j) {
277  temp = x[j];
278  i__1 = j + 1;
279  for (i__ = *n; i__ >= i__1; --i__) {
280  temp -= a_ref(i__, j) * x[i__];
281 /* L130: */
282  }
283  if (nounit) {
284  temp /= a_ref(j, j);
285  }
286  x[j] = temp;
287 /* L140: */
288  }
289  } else {
290  kx += (*n - 1) * *incx;
291  jx = kx;
292  for (j = *n; j >= 1; --j) {
293  temp = x[jx];
294  ix = kx;
295  i__1 = j + 1;
296  for (i__ = *n; i__ >= i__1; --i__) {
297  temp -= a_ref(i__, j) * x[ix];
298  ix -= *incx;
299 /* L150: */
300  }
301  if (nounit) {
302  temp /= a_ref(j, j);
303  }
304  x[jx] = temp;
305  jx -= *incx;
306 /* L160: */
307  }
308  }
309  }
310  }
311  return 0;
312 /* End of DTRSV . */
313 } /* dtrsv_ */
314 #undef a_ref
315 
316 #endif
#define a_ref(a_1, a_2)
int integer
Definition: template_blas_common.h:38
#define maxMACRO(a, b)
Definition: template_blas_common.h:43
int template_blas_erbla(const char *srname, integer *info)
Definition: template_blas_common.cc:144
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:40
bool logical
Definition: template_blas_common.h:39
logical template_blas_lsame(const char *ca, const char *cb)
Definition: template_blas_common.cc:44