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