ergo
template_lapack_larrk.h
Go to the documentation of this file.
1 /* Ergo, version 3.3, a program for linear scaling electronic structure
2  * calculations.
3  * Copyright (C) 2013 Elias Rudberg, Emanuel H. Rubensson, and Pawel Salek.
4  *
5  * This program is free software: you can redistribute it and/or modify
6  * it under the terms of the GNU General Public License as published by
7  * the Free Software Foundation, either version 3 of the License, or
8  * (at your option) any later version.
9  *
10  * This program is distributed in the hope that it will be useful,
11  * but WITHOUT ANY WARRANTY; without even the implied warranty of
12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13  * GNU General Public License for more details.
14  *
15  * You should have received a copy of the GNU General Public License
16  * along with this program. If not, see <http://www.gnu.org/licenses/>.
17  *
18  * Primary academic reference:
19  * Kohn−Sham Density Functional Theory Electronic Structure Calculations
20  * with Linearly Scaling Computational Time and Memory Usage,
21  * Elias Rudberg, Emanuel H. Rubensson, and Pawel Salek,
22  * J. Chem. Theory Comput. 7, 340 (2011),
23  * <http://dx.doi.org/10.1021/ct100611z>
24  *
25  * For further information about Ergo, see <http://www.ergoscf.org>.
26  */
27 
28  /* This file belongs to the template_lapack part of the Ergo source
29  * code. The source files in the template_lapack directory are modified
30  * versions of files originally distributed as CLAPACK, see the
31  * Copyright/license notice in the file template_lapack/COPYING.
32  */
33 
34 
35 #ifndef TEMPLATE_LAPACK_LARRK_HEADER
36 #define TEMPLATE_LAPACK_LARRK_HEADER
37 
38 template<class Treal>
39 int template_lapack_larrk(integer *n, integer *iw, Treal *gl,
40  Treal *gu, Treal *d__, Treal *e2, Treal *pivmin,
41  Treal *reltol, Treal *w, Treal *werr, integer *info)
42 {
43  /* System generated locals */
44  integer i__1;
45  Treal d__1, d__2;
46 
47 
48  /* Local variables */
49  integer i__, it;
50  Treal mid, eps, tmp1, tmp2, left, atoli, right;
51  integer itmax;
52  Treal rtoli, tnorm;
53  integer negcnt;
54 
55 
56 /* -- LAPACK auxiliary routine (version 3.2) -- */
57 /* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
58 /* November 2006 */
59 
60 /* .. Scalar Arguments .. */
61 /* .. */
62 /* .. Array Arguments .. */
63 /* .. */
64 
65 /* Purpose */
66 /* ======= */
67 
68 /* DLARRK computes one eigenvalue of a symmetric tridiagonal */
69 /* matrix T to suitable accuracy. This is an auxiliary code to be */
70 /* called from DSTEMR. */
71 
72 /* To avoid overflow, the matrix must be scaled so that its */
73 /* largest element is no greater than overflow**(1/2) * */
74 /* underflow**(1/4) in absolute value, and for greatest */
75 /* accuracy, it should not be much smaller than that. */
76 
77 /* See W. Kahan "Accurate Eigenvalues of a Symmetric Tridiagonal */
78 /* Matrix", Report CS41, Computer Science Dept., Stanford */
79 /* University, July 21, 1966. */
80 
81 /* Arguments */
82 /* ========= */
83 
84 /* N (input) INTEGER */
85 /* The order of the tridiagonal matrix T. N >= 0. */
86 
87 /* IW (input) INTEGER */
88 /* The index of the eigenvalues to be returned. */
89 
90 /* GL (input) DOUBLE PRECISION */
91 /* GU (input) DOUBLE PRECISION */
92 /* An upper and a lower bound on the eigenvalue. */
93 
94 /* D (input) DOUBLE PRECISION array, dimension (N) */
95 /* The n diagonal elements of the tridiagonal matrix T. */
96 
97 /* E2 (input) DOUBLE PRECISION array, dimension (N-1) */
98 /* The (n-1) squared off-diagonal elements of the tridiagonal matrix T. */
99 
100 /* PIVMIN (input) DOUBLE PRECISION */
101 /* The minimum pivot allowed in the Sturm sequence for T. */
102 
103 /* RELTOL (input) DOUBLE PRECISION */
104 /* The minimum relative width of an interval. When an interval */
105 /* is narrower than RELTOL times the larger (in */
106 /* magnitude) endpoint, then it is considered to be */
107 /* sufficiently small, i.e., converged. Note: this should */
108 /* always be at least radix*machine epsilon. */
109 
110 /* W (output) DOUBLE PRECISION */
111 
112 /* WERR (output) DOUBLE PRECISION */
113 /* The error bound on the corresponding eigenvalue approximation */
114 /* in W. */
115 
116 /* INFO (output) INTEGER */
117 /* = 0: Eigenvalue converged */
118 /* = -1: Eigenvalue did NOT converge */
119 
120 /* Internal Parameters */
121 /* =================== */
122 
123 /* FUDGE DOUBLE PRECISION, default = 2 */
124 /* A "fudge factor" to widen the Gershgorin intervals. */
125 
126 /* ===================================================================== */
127 
128 /* .. Parameters .. */
129 /* .. */
130 /* .. Local Scalars .. */
131 /* .. */
132 /* .. External Functions .. */
133 /* .. */
134 /* .. Intrinsic Functions .. */
135 /* .. */
136 /* .. Executable Statements .. */
137 
138 /* Get machine constants */
139  /* Parameter adjustments */
140  --e2;
141  --d__;
142 
143  /* Function Body */
144  eps = template_lapack_lamch("P", (Treal)0);
145 /* Computing MAX */
146  d__1 = absMACRO(*gl), d__2 = absMACRO(*gu);
147  tnorm = maxMACRO(d__1,d__2);
148  rtoli = *reltol;
149  atoli = *pivmin * 4.;
150  itmax = (integer) ((template_blas_log(tnorm + *pivmin) - template_blas_log(*pivmin)) / template_blas_log(2.)) + 2;
151  *info = -1;
152  left = *gl - tnorm * 2. * eps * *n - *pivmin * 4.;
153  right = *gu + tnorm * 2. * eps * *n + *pivmin * 4.;
154  it = 0;
155 L10:
156 
157 /* Check if interval converged or maximum number of iterations reached */
158 
159  tmp1 = (d__1 = right - left, absMACRO(d__1));
160 /* Computing MAX */
161  d__1 = absMACRO(right), d__2 = absMACRO(left);
162  tmp2 = maxMACRO(d__1,d__2);
163 /* Computing MAX */
164  d__1 = maxMACRO(atoli,*pivmin), d__2 = rtoli * tmp2;
165  if (tmp1 < maxMACRO(d__1,d__2)) {
166  *info = 0;
167  goto L30;
168  }
169  if (it > itmax) {
170  goto L30;
171  }
172 
173 /* Count number of negative pivots for mid-point */
174 
175  ++it;
176  mid = (left + right) * .5;
177  negcnt = 0;
178  tmp1 = d__[1] - mid;
179  if (absMACRO(tmp1) < *pivmin) {
180  tmp1 = -(*pivmin);
181  }
182  if (tmp1 <= 0.) {
183  ++negcnt;
184  }
185 
186  i__1 = *n;
187  for (i__ = 2; i__ <= i__1; ++i__) {
188  tmp1 = d__[i__] - e2[i__ - 1] / tmp1 - mid;
189  if (absMACRO(tmp1) < *pivmin) {
190  tmp1 = -(*pivmin);
191  }
192  if (tmp1 <= 0.) {
193  ++negcnt;
194  }
195 /* L20: */
196  }
197  if (negcnt >= *iw) {
198  right = mid;
199  } else {
200  left = mid;
201  }
202  goto L10;
203 L30:
204 
205 /* Converged or maximum number of iterations reached */
206 
207  *w = (left + right) * .5;
208  *werr = (d__1 = right - left, absMACRO(d__1)) * .5;
209  return 0;
210 
211 /* End of DLARRK */
212 
213 } /* dlarrk_ */
214 
215 #endif
static const real gu
Definition: fun-pz81.c:71
#define absMACRO(x)
Definition: template_blas_common.h:45
Definition: Matrix.h:73
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:39
int integer
Definition: template_blas_common.h:38
#define maxMACRO(a, b)
Definition: template_blas_common.h:43
Treal template_blas_log(Treal x)
Treal template_lapack_lamch(const char *cmach, Treal dummyReal)
Definition: template_lapack_lamch.h:199
Definition: Matrix.h:73