ergo
template_lapack_lagtf.h
Go to the documentation of this file.
1
/* Ergo, version 3.8, a program for linear scaling electronic structure
2
* calculations.
3
* Copyright (C) 2019 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
41
template
<
class
Treal>
42
int
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
template_lapack_lamch
Treal template_lapack_lamch(const char *cmach, Treal dummyReal)
Definition:
template_lapack_lamch.h:202
absMACRO
#define absMACRO(x)
Definition:
template_blas_common.h:47
template_lapack_lagtf
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
template_blas_erbla
int template_blas_erbla(const char *srname, integer *info)
Definition:
template_blas_common.cc:146
integer
int integer
Definition:
template_blas_common.h:40
maxMACRO
#define maxMACRO(a, b)
Definition:
template_blas_common.h:45
source
matrix
template_lapack
lapack
template_lapack_lagtf.h
Generated on Sun Aug 16 2020 00:00:00 for ergo by
1.8.18