ergo
template_lapack_lagts.h
Go to the documentation of this file.
1/* Ergo, version 3.8.2, a program for linear scaling electronic structure
2 * calculations.
3 * Copyright (C) 2023 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
41template<class Treal>
42int 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);
252L40:
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);
319L70:
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
int template_blas_erbla(const char *srname, integer *info)
Definition template_blas_common.cc:146
int integer
Definition template_blas_common.h:40
#define absMACRO(x)
Definition template_blas_common.h:47
#define maxMACRO(a, b)
Definition template_blas_common.h:45
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
Treal template_lapack_d_sign(const Treal *a, const Treal *b)
Definition template_lapack_lamch.h:48