ergo
template_lapack_larrk.h
Go to the documentation of this file.
1 /* Ergo, version 3.7, a program for linear scaling electronic structure
2  * calculations.
3  * Copyright (C) 2018 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
static const real gu
Definition: fun-pz81.c:68
#define absMACRO(x)
Definition: template_blas_common.h:47
Definition: Matrix.h:75
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
int integer
Definition: template_blas_common.h:40
#define maxMACRO(a, b)
Definition: template_blas_common.h:45
Treal template_blas_log(Treal x)
Treal template_lapack_lamch(const char *cmach, Treal dummyReal)
Definition: template_lapack_lamch.h:202
Definition: Matrix.h:75