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