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