ergo
template_lapack_larrk.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_LARRK_HEADER
38
#define TEMPLATE_LAPACK_LARRK_HEADER
39
40
template
<
class
Treal>
41
int
template_lapack_larrk
(
integer
*n,
integer
*iw, Treal *gl,
42
Treal *
gu
, Treal *d__, Treal *e2, Treal *pivmin,
43
Treal *reltol, Treal *w, Treal *werr,
integer
*info)
44
{
45
/* System generated locals */
46
integer
i__1;
47
Treal d__1, d__2;
48
49
50
/* Local variables */
51
integer
i__, it;
52
Treal mid, eps, tmp1, tmp2,
left
, atoli,
right
;
53
integer
itmax;
54
Treal rtoli, tnorm;
55
integer
negcnt;
56
57
58
/* -- LAPACK auxiliary routine (version 3.2) -- */
59
/* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
60
/* November 2006 */
61
62
/* .. Scalar Arguments .. */
63
/* .. */
64
/* .. Array Arguments .. */
65
/* .. */
66
67
/* Purpose */
68
/* ======= */
69
70
/* DLARRK computes one eigenvalue of a symmetric tridiagonal */
71
/* matrix T to suitable accuracy. This is an auxiliary code to be */
72
/* called from DSTEMR. */
73
74
/* To avoid overflow, the matrix must be scaled so that its */
75
/* largest element is no greater than overflow**(1/2) * */
76
/* underflow**(1/4) in absolute value, and for greatest */
77
/* accuracy, it should not be much smaller than that. */
78
79
/* See W. Kahan "Accurate Eigenvalues of a Symmetric Tridiagonal */
80
/* Matrix", Report CS41, Computer Science Dept., Stanford */
81
/* University, July 21, 1966. */
82
83
/* Arguments */
84
/* ========= */
85
86
/* N (input) INTEGER */
87
/* The order of the tridiagonal matrix T. N >= 0. */
88
89
/* IW (input) INTEGER */
90
/* The index of the eigenvalues to be returned. */
91
92
/* GL (input) DOUBLE PRECISION */
93
/* GU (input) DOUBLE PRECISION */
94
/* An upper and a lower bound on the eigenvalue. */
95
96
/* D (input) DOUBLE PRECISION array, dimension (N) */
97
/* The n diagonal elements of the tridiagonal matrix T. */
98
99
/* E2 (input) DOUBLE PRECISION array, dimension (N-1) */
100
/* The (n-1) squared off-diagonal elements of the tridiagonal matrix T. */
101
102
/* PIVMIN (input) DOUBLE PRECISION */
103
/* The minimum pivot allowed in the Sturm sequence for T. */
104
105
/* RELTOL (input) DOUBLE PRECISION */
106
/* The minimum relative width of an interval. When an interval */
107
/* is narrower than RELTOL times the larger (in */
108
/* magnitude) endpoint, then it is considered to be */
109
/* sufficiently small, i.e., converged. Note: this should */
110
/* always be at least radix*machine epsilon. */
111
112
/* W (output) DOUBLE PRECISION */
113
114
/* WERR (output) DOUBLE PRECISION */
115
/* The error bound on the corresponding eigenvalue approximation */
116
/* in W. */
117
118
/* INFO (output) INTEGER */
119
/* = 0: Eigenvalue converged */
120
/* = -1: Eigenvalue did NOT converge */
121
122
/* Internal Parameters */
123
/* =================== */
124
125
/* FUDGE DOUBLE PRECISION, default = 2 */
126
/* A "fudge factor" to widen the Gershgorin intervals. */
127
128
/* ===================================================================== */
129
130
/* .. Parameters .. */
131
/* .. */
132
/* .. Local Scalars .. */
133
/* .. */
134
/* .. External Functions .. */
135
/* .. */
136
/* .. Intrinsic Functions .. */
137
/* .. */
138
/* .. Executable Statements .. */
139
140
/* Get machine constants */
141
/* Parameter adjustments */
142
--e2;
143
--d__;
144
145
/* Function Body */
146
eps =
template_lapack_lamch
(
"P"
, (Treal)0);
147
/* Computing MAX */
148
d__1 =
absMACRO
(*gl), d__2 =
absMACRO
(*
gu
);
149
tnorm =
maxMACRO
(d__1,d__2);
150
rtoli = *reltol;
151
atoli = *pivmin * 4.;
152
itmax = (
integer
) ((
template_blas_log
(tnorm + *pivmin) -
template_blas_log
(*pivmin)) /
template_blas_log
(2.)) + 2;
153
*info = -1;
154
left
= *gl - tnorm * 2. * eps * *n - *pivmin * 4.;
155
right
= *
gu
+ tnorm * 2. * eps * *n + *pivmin * 4.;
156
it = 0;
157
L10:
158
159
/* Check if interval converged or maximum number of iterations reached */
160
161
tmp1 = (d__1 =
right
-
left
,
absMACRO
(d__1));
162
/* Computing MAX */
163
d__1 =
absMACRO
(
right
), d__2 =
absMACRO
(
left
);
164
tmp2 =
maxMACRO
(d__1,d__2);
165
/* Computing MAX */
166
d__1 =
maxMACRO
(atoli,*pivmin), d__2 = rtoli * tmp2;
167
if
(tmp1 <
maxMACRO
(d__1,d__2)) {
168
*info = 0;
169
goto
L30;
170
}
171
if
(it > itmax) {
172
goto
L30;
173
}
174
175
/* Count number of negative pivots for mid-point */
176
177
++it;
178
mid = (
left
+
right
) * .5;
179
negcnt = 0;
180
tmp1 = d__[1] - mid;
181
if
(
absMACRO
(tmp1) < *pivmin) {
182
tmp1 = -(*pivmin);
183
}
184
if
(tmp1 <= 0.) {
185
++negcnt;
186
}
187
188
i__1 = *n;
189
for
(i__ = 2; i__ <= i__1; ++i__) {
190
tmp1 = d__[i__] - e2[i__ - 1] / tmp1 - mid;
191
if
(
absMACRO
(tmp1) < *pivmin) {
192
tmp1 = -(*pivmin);
193
}
194
if
(tmp1 <= 0.) {
195
++negcnt;
196
}
197
/* L20: */
198
}
199
if
(negcnt >= *iw) {
200
right
= mid;
201
}
else
{
202
left
= mid;
203
}
204
goto
L10;
205
L30:
206
207
/* Converged or maximum number of iterations reached */
208
209
*w = (
left
+
right
) * .5;
210
*werr = (d__1 =
right
-
left
,
absMACRO
(d__1)) * .5;
211
return
0;
212
213
/* End of DLARRK */
214
215
}
/* dlarrk_ */
216
217
#endif
gu
static const real gu
Definition:
fun-pz81.c:68
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
mat::left
@ left
Definition:
Matrix.h:75
template_lapack_larrk
int template_lapack_larrk(integer *n, integer *iw, Treal *gl, Treal *gu, Treal *d__, Treal *e2, Treal *pivmin, Treal *reltol, Treal *w, Treal *werr, integer *info)
Definition:
template_lapack_larrk.h:41
mat::right
@ right
Definition:
Matrix.h:75
template_blas_log
Treal template_blas_log(Treal x)
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_larrk.h
Generated on Sun Aug 16 2020 00:00:00 for ergo by
1.8.18