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
40template<class Treal>
41int 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;
157L10:
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;
205L30:
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
static const real gu
Definition: fun-pz81.c:68
@ right
Definition: Matrix.h:75
@ left
Definition: Matrix.h:75
Treal template_blas_log(Treal x)
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
Treal template_lapack_lamch(const char *cmach, Treal dummyReal)
Definition: template_lapack_lamch.h:202
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