ergo
template_lapack_lagtf.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_LAGTF_HEADER
38#define TEMPLATE_LAPACK_LAGTF_HEADER
39
40
41template<class Treal>
42int template_lapack_lagtf(const integer *n, Treal *a, const Treal *lambda,
43 Treal *b, Treal *c__, const Treal *tol, Treal *d__,
44 integer *in, integer *info)
45{
46/* -- LAPACK routine (version 3.0) --
47 Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
48 Courant Institute, Argonne National Lab, and Rice University
49 June 30, 1999
50
51
52 Purpose
53 =======
54
55 DLAGTF factorizes the matrix (T - lambda*I), where T is an n by n
56 tridiagonal matrix and lambda is a scalar, as
57
58 T - lambda*I = PLU,
59
60 where P is a permutation matrix, L is a unit lower tridiagonal matrix
61 with at most one non-zero sub-diagonal elements per column and U is
62 an upper triangular matrix with at most two non-zero super-diagonal
63 elements per column.
64
65 The factorization is obtained by Gaussian elimination with partial
66 pivoting and implicit row scaling.
67
68 The parameter LAMBDA is included in the routine so that DLAGTF may
69 be used, in conjunction with DLAGTS, to obtain eigenvectors of T by
70 inverse iteration.
71
72 Arguments
73 =========
74
75 N (input) INTEGER
76 The order of the matrix T.
77
78 A (input/output) DOUBLE PRECISION array, dimension (N)
79 On entry, A must contain the diagonal elements of T.
80
81 On exit, A is overwritten by the n diagonal elements of the
82 upper triangular matrix U of the factorization of T.
83
84 LAMBDA (input) DOUBLE PRECISION
85 On entry, the scalar lambda.
86
87 B (input/output) DOUBLE PRECISION array, dimension (N-1)
88 On entry, B must contain the (n-1) super-diagonal elements of
89 T.
90
91 On exit, B is overwritten by the (n-1) super-diagonal
92 elements of the matrix U of the factorization of T.
93
94 C (input/output) DOUBLE PRECISION array, dimension (N-1)
95 On entry, C must contain the (n-1) sub-diagonal elements of
96 T.
97
98 On exit, C is overwritten by the (n-1) sub-diagonal elements
99 of the matrix L of the factorization of T.
100
101 TOL (input) DOUBLE PRECISION
102 On entry, a relative tolerance used to indicate whether or
103 not the matrix (T - lambda*I) is nearly singular. TOL should
104 normally be chose as approximately the largest relative error
105 in the elements of T. For example, if the elements of T are
106 correct to about 4 significant figures, then TOL should be
107 set to about 5*10**(-4). If TOL is supplied as less than eps,
108 where eps is the relative machine precision, then the value
109 eps is used in place of TOL.
110
111 D (output) DOUBLE PRECISION array, dimension (N-2)
112 On exit, D is overwritten by the (n-2) second super-diagonal
113 elements of the matrix U of the factorization of T.
114
115 IN (output) INTEGER array, dimension (N)
116 On exit, IN contains details of the permutation matrix P. If
117 an interchange occurred at the kth step of the elimination,
118 then IN(k) = 1, otherwise IN(k) = 0. The element IN(n)
119 returns the smallest positive integer j such that
120
121 abs( u(j,j) ).le. norm( (T - lambda*I)(j) )*TOL,
122
123 where norm( A(j) ) denotes the sum of the absolute values of
124 the jth row of the matrix A. If no such j exists then IN(n)
125 is returned as zero. If IN(n) is returned as positive, then a
126 diagonal element of U is small, indicating that
127 (T - lambda*I) is singular or nearly singular,
128
129 INFO (output) INTEGER
130 = 0 : successful exit
131 .lt. 0: if INFO = -k, the kth argument had an illegal value
132
133 =====================================================================
134
135
136 Parameter adjustments */
137 /* System generated locals */
138 integer i__1;
139 Treal d__1, d__2;
140 /* Local variables */
141 Treal temp, mult;
142 integer k;
143 Treal scale1, scale2;
144 Treal tl;
145 Treal eps, piv1, piv2;
146
147 --in;
148 --d__;
149 --c__;
150 --b;
151 --a;
152
153 /* Function Body */
154 *info = 0;
155 if (*n < 0) {
156 *info = -1;
157 i__1 = -(*info);
158 template_blas_erbla("LAGTF ", &i__1);
159 return 0;
160 }
161
162 if (*n == 0) {
163 return 0;
164 }
165
166 a[1] -= *lambda;
167 in[*n] = 0;
168 if (*n == 1) {
169 if (a[1] == 0.) {
170 in[1] = 1;
171 }
172 return 0;
173 }
174
175 eps = template_lapack_lamch("Epsilon", (Treal)0);
176
177 tl = maxMACRO(*tol,eps);
178 scale1 = absMACRO(a[1]) + absMACRO(b[1]);
179 i__1 = *n - 1;
180 for (k = 1; k <= i__1; ++k) {
181 a[k + 1] -= *lambda;
182 scale2 = (d__1 = c__[k], absMACRO(d__1)) + (d__2 = a[k + 1], absMACRO(d__2));
183 if (k < *n - 1) {
184 scale2 += (d__1 = b[k + 1], absMACRO(d__1));
185 }
186 if (a[k] == 0.) {
187 piv1 = 0.;
188 } else {
189 piv1 = (d__1 = a[k], absMACRO(d__1)) / scale1;
190 }
191 if (c__[k] == 0.) {
192 in[k] = 0;
193 piv2 = 0.;
194 scale1 = scale2;
195 if (k < *n - 1) {
196 d__[k] = 0.;
197 }
198 } else {
199 piv2 = (d__1 = c__[k], absMACRO(d__1)) / scale2;
200 if (piv2 <= piv1) {
201 in[k] = 0;
202 scale1 = scale2;
203 c__[k] /= a[k];
204 a[k + 1] -= c__[k] * b[k];
205 if (k < *n - 1) {
206 d__[k] = 0.;
207 }
208 } else {
209 in[k] = 1;
210 mult = a[k] / c__[k];
211 a[k] = c__[k];
212 temp = a[k + 1];
213 a[k + 1] = b[k] - mult * temp;
214 if (k < *n - 1) {
215 d__[k] = b[k + 1];
216 b[k + 1] = -mult * d__[k];
217 }
218 b[k] = temp;
219 c__[k] = mult;
220 }
221 }
222 if (maxMACRO(piv1,piv2) <= tl && in[*n] == 0) {
223 in[*n] = k;
224 }
225/* L10: */
226 }
227 if ((d__1 = a[*n], absMACRO(d__1)) <= scale1 * tl && in[*n] == 0) {
228 in[*n] = *n;
229 }
230
231 return 0;
232
233/* End of DLAGTF */
234
235} /* dlagtf_ */
236
237#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_lagtf(const integer *n, Treal *a, const Treal *lambda, Treal *b, Treal *c__, const Treal *tol, Treal *d__, integer *in, integer *info)
Definition: template_lapack_lagtf.h:42
Treal template_lapack_lamch(const char *cmach, Treal dummyReal)
Definition: template_lapack_lamch.h:202