ergo
template_lapack_lagts.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_LAPACK_LAGTS_HEADER
36 #define TEMPLATE_LAPACK_LAGTS_HEADER
37 
38 
39 template<class Treal>
40 int template_lapack_lagts(const integer *job, const integer *n, const Treal *a,
41  const Treal *b, const Treal *c__, const Treal *d__, const integer *in,
42  Treal *y, Treal *tol, integer *info)
43 {
44 /* -- LAPACK auxiliary routine (version 3.0) --
45  Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
46  Courant Institute, Argonne National Lab, and Rice University
47  October 31, 1992
48 
49 
50  Purpose
51  =======
52 
53  DLAGTS may be used to solve one of the systems of equations
54 
55  (T - lambda*I)*x = y or (T - lambda*I)'*x = y,
56 
57  where T is an n by n tridiagonal matrix, for x, following the
58  factorization of (T - lambda*I) as
59 
60  (T - lambda*I) = P*L*U ,
61 
62  by routine DLAGTF. The choice of equation to be solved is
63  controlled by the argument JOB, and in each case there is an option
64  to perturb zero or very small diagonal elements of U, this option
65  being intended for use in applications such as inverse iteration.
66 
67  Arguments
68  =========
69 
70  JOB (input) INTEGER
71  Specifies the job to be performed by DLAGTS as follows:
72  = 1: The equations (T - lambda*I)x = y are to be solved,
73  but diagonal elements of U are not to be perturbed.
74  = -1: The equations (T - lambda*I)x = y are to be solved
75  and, if overflow would otherwise occur, the diagonal
76  elements of U are to be perturbed. See argument TOL
77  below.
78  = 2: The equations (T - lambda*I)'x = y are to be solved,
79  but diagonal elements of U are not to be perturbed.
80  = -2: The equations (T - lambda*I)'x = y are to be solved
81  and, if overflow would otherwise occur, the diagonal
82  elements of U are to be perturbed. See argument TOL
83  below.
84 
85  N (input) INTEGER
86  The order of the matrix T.
87 
88  A (input) DOUBLE PRECISION array, dimension (N)
89  On entry, A must contain the diagonal elements of U as
90  returned from DLAGTF.
91 
92  B (input) DOUBLE PRECISION array, dimension (N-1)
93  On entry, B must contain the first super-diagonal elements of
94  U as returned from DLAGTF.
95 
96  C (input) DOUBLE PRECISION array, dimension (N-1)
97  On entry, C must contain the sub-diagonal elements of L as
98  returned from DLAGTF.
99 
100  D (input) DOUBLE PRECISION array, dimension (N-2)
101  On entry, D must contain the second super-diagonal elements
102  of U as returned from DLAGTF.
103 
104  IN (input) INTEGER array, dimension (N)
105  On entry, IN must contain details of the matrix P as returned
106  from DLAGTF.
107 
108  Y (input/output) DOUBLE PRECISION array, dimension (N)
109  On entry, the right hand side vector y.
110  On exit, Y is overwritten by the solution vector x.
111 
112  TOL (input/output) DOUBLE PRECISION
113  On entry, with JOB .lt. 0, TOL should be the minimum
114  perturbation to be made to very small diagonal elements of U.
115  TOL should normally be chosen as about eps*norm(U), where eps
116  is the relative machine precision, but if TOL is supplied as
117  non-positive, then it is reset to eps*max( abs( u(i,j) ) ).
118  If JOB .gt. 0 then TOL is not referenced.
119 
120  On exit, TOL is changed as described above, only if TOL is
121  non-positive on entry. Otherwise TOL is unchanged.
122 
123  INFO (output) INTEGER
124  = 0 : successful exit
125  .lt. 0: if INFO = -i, the i-th argument had an illegal value
126  .gt. 0: overflow would occur when computing the INFO(th)
127  element of the solution vector x. This can only occur
128  when JOB is supplied as positive and either means
129  that a diagonal element of U is very small, or that
130  the elements of the right-hand side vector y are very
131  large.
132 
133  =====================================================================
134 
135 
136  Parameter adjustments */
137  /* System generated locals */
138  integer i__1;
139  Treal d__1, d__2, d__3, d__4, d__5;
140  /* Local variables */
141  Treal temp, pert;
142  integer k;
143  Treal absak, sfmin, ak;
144  Treal bignum, eps;
145 
146  --y;
147  --in;
148  --d__;
149  --c__;
150  --b;
151  --a;
152 
153  /* Function Body */
154  *info = 0;
155  if (absMACRO(*job) > 2 || *job == 0) {
156  *info = -1;
157  } else if (*n < 0) {
158  *info = -2;
159  }
160  if (*info != 0) {
161  i__1 = -(*info);
162  template_blas_erbla("LAGTS ", &i__1);
163  return 0;
164  }
165 
166  if (*n == 0) {
167  return 0;
168  }
169 
170  eps = template_lapack_lamch("Epsilon", (Treal)0);
171  sfmin = template_lapack_lamch("Safe minimum", (Treal)0);
172  bignum = 1. / sfmin;
173 
174  if (*job < 0) {
175  if (*tol <= 0.) {
176  *tol = absMACRO(a[1]);
177  if (*n > 1) {
178 /* Computing MAX */
179  d__1 = *tol, d__2 = absMACRO(a[2]), d__1 = maxMACRO(d__1,d__2), d__2 =
180  absMACRO(b[1]);
181  *tol = maxMACRO(d__1,d__2);
182  }
183  i__1 = *n;
184  for (k = 3; k <= i__1; ++k) {
185 /* Computing MAX */
186  d__4 = *tol, d__5 = (d__1 = a[k], absMACRO(d__1)), d__4 = maxMACRO(d__4,
187  d__5), d__5 = (d__2 = b[k - 1], absMACRO(d__2)), d__4 =
188  maxMACRO(d__4,d__5), d__5 = (d__3 = d__[k - 2], absMACRO(d__3));
189  *tol = maxMACRO(d__4,d__5);
190 /* L10: */
191  }
192  *tol *= eps;
193  if (*tol == 0.) {
194  *tol = eps;
195  }
196  }
197  }
198 
199  if (absMACRO(*job) == 1) {
200  i__1 = *n;
201  for (k = 2; k <= i__1; ++k) {
202  if (in[k - 1] == 0) {
203  y[k] -= c__[k - 1] * y[k - 1];
204  } else {
205  temp = y[k - 1];
206  y[k - 1] = y[k];
207  y[k] = temp - c__[k - 1] * y[k];
208  }
209 /* L20: */
210  }
211  if (*job == 1) {
212  for (k = *n; k >= 1; --k) {
213  if (k <= *n - 2) {
214  temp = y[k] - b[k] * y[k + 1] - d__[k] * y[k + 2];
215  } else if (k == *n - 1) {
216  temp = y[k] - b[k] * y[k + 1];
217  } else {
218  temp = y[k];
219  }
220  ak = a[k];
221  absak = absMACRO(ak);
222  if (absak < 1.) {
223  if (absak < sfmin) {
224  if (absak == 0. || absMACRO(temp) * sfmin > absak) {
225  *info = k;
226  return 0;
227  } else {
228  temp *= bignum;
229  ak *= bignum;
230  }
231  } else if (absMACRO(temp) > absak * bignum) {
232  *info = k;
233  return 0;
234  }
235  }
236  y[k] = temp / ak;
237 /* L30: */
238  }
239  } else {
240  for (k = *n; k >= 1; --k) {
241  if (k <= *n - 2) {
242  temp = y[k] - b[k] * y[k + 1] - d__[k] * y[k + 2];
243  } else if (k == *n - 1) {
244  temp = y[k] - b[k] * y[k + 1];
245  } else {
246  temp = y[k];
247  }
248  ak = a[k];
249  pert = template_lapack_d_sign(tol, &ak);
250 L40:
251  absak = absMACRO(ak);
252  if (absak < 1.) {
253  if (absak < sfmin) {
254  if (absak == 0. || absMACRO(temp) * sfmin > absak) {
255  ak += pert;
256  pert *= 2;
257  goto L40;
258  } else {
259  temp *= bignum;
260  ak *= bignum;
261  }
262  } else if (absMACRO(temp) > absak * bignum) {
263  ak += pert;
264  pert *= 2;
265  goto L40;
266  }
267  }
268  y[k] = temp / ak;
269 /* L50: */
270  }
271  }
272  } else {
273 
274 /* Come to here if JOB = 2 or -2 */
275 
276  if (*job == 2) {
277  i__1 = *n;
278  for (k = 1; k <= i__1; ++k) {
279  if (k >= 3) {
280  temp = y[k] - b[k - 1] * y[k - 1] - d__[k - 2] * y[k - 2];
281  } else if (k == 2) {
282  temp = y[k] - b[k - 1] * y[k - 1];
283  } else {
284  temp = y[k];
285  }
286  ak = a[k];
287  absak = absMACRO(ak);
288  if (absak < 1.) {
289  if (absak < sfmin) {
290  if (absak == 0. || absMACRO(temp) * sfmin > absak) {
291  *info = k;
292  return 0;
293  } else {
294  temp *= bignum;
295  ak *= bignum;
296  }
297  } else if (absMACRO(temp) > absak * bignum) {
298  *info = k;
299  return 0;
300  }
301  }
302  y[k] = temp / ak;
303 /* L60: */
304  }
305  } else {
306  i__1 = *n;
307  for (k = 1; k <= i__1; ++k) {
308  if (k >= 3) {
309  temp = y[k] - b[k - 1] * y[k - 1] - d__[k - 2] * y[k - 2];
310  } else if (k == 2) {
311  temp = y[k] - b[k - 1] * y[k - 1];
312  } else {
313  temp = y[k];
314  }
315  ak = a[k];
316  pert = template_lapack_d_sign(tol, &ak);
317 L70:
318  absak = absMACRO(ak);
319  if (absak < 1.) {
320  if (absak < sfmin) {
321  if (absak == 0. || absMACRO(temp) * sfmin > absak) {
322  ak += pert;
323  pert *= 2;
324  goto L70;
325  } else {
326  temp *= bignum;
327  ak *= bignum;
328  }
329  } else if (absMACRO(temp) > absak * bignum) {
330  ak += pert;
331  pert *= 2;
332  goto L70;
333  }
334  }
335  y[k] = temp / ak;
336 /* L80: */
337  }
338  }
339 
340  for (k = *n; k >= 2; --k) {
341  if (in[k - 1] == 0) {
342  y[k - 1] -= c__[k - 1] * y[k];
343  } else {
344  temp = y[k - 1];
345  y[k - 1] = y[k];
346  y[k] = temp - c__[k - 1] * y[k];
347  }
348 /* L90: */
349  }
350  }
351 
352 /* End of DLAGTS */
353 
354  return 0;
355 } /* dlagts_ */
356 
357 #endif
#define absMACRO(x)
Definition: template_blas_common.h:45
int integer
Definition: template_blas_common.h:38
#define maxMACRO(a, b)
Definition: template_blas_common.h:43
Treal template_lapack_d_sign(const Treal *a, const Treal *b)
Definition: template_lapack_lamch.h:46
int template_blas_erbla(const char *srname, integer *info)
Definition: template_blas_common.cc:144
int template_lapack_lagts(const integer *job, const integer *n, const Treal *a, const Treal *b, const Treal *c__, const Treal *d__, const integer *in, Treal *y, Treal *tol, integer *info)
Definition: template_lapack_lagts.h:40
Treal template_lapack_lamch(const char *cmach, Treal dummyReal)
Definition: template_lapack_lamch.h:199