ergo
template_lapack_larrb.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_LARRB_HEADER
38 #define TEMPLATE_LAPACK_LARRB_HEADER
39 
40 
41 #include "template_lapack_laneg.h"
42 
43 
44 template<class Treal>
45 int template_lapack_larrb(integer *n, Treal *d__, Treal *lld,
46  integer *ifirst, integer *ilast, Treal *rtol1, Treal *rtol2,
47  integer *offset, Treal *w, Treal *wgap, Treal *werr,
48  Treal *work, integer *iwork, Treal *pivmin, Treal *
49  spdiam, integer *twist, integer *info)
50 {
51  /* System generated locals */
52  integer i__1;
53  Treal d__1, d__2;
54 
55  /* Local variables */
56  integer i__, k, r__, i1, ii, ip;
57  Treal gap, mid, tmp, back, lgap, rgap, left;
58  integer iter, nint, prev, next;
59  Treal cvrgd, right, width;
60  integer negcnt;
61  Treal mnwdth;
62  integer olnint, maxitr;
63 
64 
65 /* -- LAPACK auxiliary routine (version 3.2) -- */
66 /* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
67 /* November 2006 */
68 
69 /* .. Scalar Arguments .. */
70 /* .. */
71 /* .. Array Arguments .. */
72 /* .. */
73 
74 /* Purpose */
75 /* ======= */
76 
77 /* Given the relatively robust representation(RRR) L D L^T, DLARRB */
78 /* does "limited" bisection to refine the eigenvalues of L D L^T, */
79 /* W( IFIRST-OFFSET ) through W( ILAST-OFFSET ), to more accuracy. Initial */
80 /* guesses for these eigenvalues are input in W, the corresponding estimate */
81 /* of the error in these guesses and their gaps are input in WERR */
82 /* and WGAP, respectively. During bisection, intervals */
83 /* [left, right] are maintained by storing their mid-points and */
84 /* semi-widths in the arrays W and WERR respectively. */
85 
86 /* Arguments */
87 /* ========= */
88 
89 /* N (input) INTEGER */
90 /* The order of the matrix. */
91 
92 /* D (input) DOUBLE PRECISION array, dimension (N) */
93 /* The N diagonal elements of the diagonal matrix D. */
94 
95 /* LLD (input) DOUBLE PRECISION array, dimension (N-1) */
96 /* The (N-1) elements L(i)*L(i)*D(i). */
97 
98 /* IFIRST (input) INTEGER */
99 /* The index of the first eigenvalue to be computed. */
100 
101 /* ILAST (input) INTEGER */
102 /* The index of the last eigenvalue to be computed. */
103 
104 /* RTOL1 (input) DOUBLE PRECISION */
105 /* RTOL2 (input) DOUBLE PRECISION */
106 /* Tolerance for the convergence of the bisection intervals. */
107 /* An interval [LEFT,RIGHT] has converged if */
108 /* RIGHT-LEFT.LT.MAX( RTOL1*GAP, RTOL2*MAX(|LEFT|,|RIGHT|) ) */
109 /* where GAP is the (estimated) distance to the nearest */
110 /* eigenvalue. */
111 
112 /* OFFSET (input) INTEGER */
113 /* Offset for the arrays W, WGAP and WERR, i.e., the IFIRST-OFFSET */
114 /* through ILAST-OFFSET elements of these arrays are to be used. */
115 
116 /* W (input/output) DOUBLE PRECISION array, dimension (N) */
117 /* On input, W( IFIRST-OFFSET ) through W( ILAST-OFFSET ) are */
118 /* estimates of the eigenvalues of L D L^T indexed IFIRST throug */
119 /* ILAST. */
120 /* On output, these estimates are refined. */
121 
122 /* WGAP (input/output) DOUBLE PRECISION array, dimension (N-1) */
123 /* On input, the (estimated) gaps between consecutive */
124 /* eigenvalues of L D L^T, i.e., WGAP(I-OFFSET) is the gap between */
125 /* eigenvalues I and I+1. Note that if IFIRST.EQ.ILAST */
126 /* then WGAP(IFIRST-OFFSET) must be set to ZERO. */
127 /* On output, these gaps are refined. */
128 
129 /* WERR (input/output) DOUBLE PRECISION array, dimension (N) */
130 /* On input, WERR( IFIRST-OFFSET ) through WERR( ILAST-OFFSET ) are */
131 /* the errors in the estimates of the corresponding elements in W. */
132 /* On output, these errors are refined. */
133 
134 /* WORK (workspace) DOUBLE PRECISION array, dimension (2*N) */
135 /* Workspace. */
136 
137 /* IWORK (workspace) INTEGER array, dimension (2*N) */
138 /* Workspace. */
139 
140 /* PIVMIN (input) DOUBLE PRECISION */
141 /* The minimum pivot in the Sturm sequence. */
142 
143 /* SPDIAM (input) DOUBLE PRECISION */
144 /* The spectral diameter of the matrix. */
145 
146 /* TWIST (input) INTEGER */
147 /* The twist index for the twisted factorization that is used */
148 /* for the negcount. */
149 /* TWIST = N: Compute negcount from L D L^T - LAMBDA I = L+ D+ L+^T */
150 /* TWIST = 1: Compute negcount from L D L^T - LAMBDA I = U- D- U-^T */
151 /* TWIST = R: Compute negcount from L D L^T - LAMBDA I = N(r) D(r) N(r) */
152 
153 /* INFO (output) INTEGER */
154 /* Error flag. */
155 
156 /* Further Details */
157 /* =============== */
158 
159 /* Based on contributions by */
160 /* Beresford Parlett, University of California, Berkeley, USA */
161 /* Jim Demmel, University of California, Berkeley, USA */
162 /* Inderjit Dhillon, University of Texas, Austin, USA */
163 /* Osni Marques, LBNL/NERSC, USA */
164 /* Christof Voemel, University of California, Berkeley, USA */
165 
166 /* ===================================================================== */
167 
168 /* .. Parameters .. */
169 /* .. */
170 /* .. Local Scalars .. */
171 /* .. */
172 /* .. External Functions .. */
173 
174 /* .. */
175 /* .. Intrinsic Functions .. */
176 /* .. */
177 /* .. Executable Statements .. */
178 
179  /* Parameter adjustments */
180  --iwork;
181  --work;
182  --werr;
183  --wgap;
184  --w;
185  --lld;
186  --d__;
187 
188  /* Function Body */
189  *info = 0;
190 
191  maxitr = (integer) ((template_blas_log(*spdiam + *pivmin) - template_blas_log(*pivmin)) / template_blas_log(2.)) +
192  2;
193  mnwdth = *pivmin * 2.;
194 
195  r__ = *twist;
196  if (r__ < 1 || r__ > *n) {
197  r__ = *n;
198  }
199 
200 /* Initialize unconverged intervals in [ WORK(2*I-1), WORK(2*I) ]. */
201 /* The Sturm Count, Count( WORK(2*I-1) ) is arranged to be I-1, while */
202 /* Count( WORK(2*I) ) is stored in IWORK( 2*I ). The integer IWORK( 2*I-1 ) */
203 /* for an unconverged interval is set to the index of the next unconverged */
204 /* interval, and is -1 or 0 for a converged interval. Thus a linked */
205 /* list of unconverged intervals is set up. */
206 
207  i1 = *ifirst;
208 /* The number of unconverged intervals */
209  nint = 0;
210 /* The last unconverged interval found */
211  prev = 0;
212  rgap = wgap[i1 - *offset];
213  i__1 = *ilast;
214  for (i__ = i1; i__ <= i__1; ++i__) {
215  k = i__ << 1;
216  ii = i__ - *offset;
217  left = w[ii] - werr[ii];
218  right = w[ii] + werr[ii];
219  lgap = rgap;
220  rgap = wgap[ii];
221  gap = minMACRO(lgap,rgap);
222 /* Make sure that [LEFT,RIGHT] contains the desired eigenvalue */
223 /* Compute negcount from dstqds facto L+D+L+^T = L D L^T - LEFT */
224 
225 /* Do while( NEGCNT(LEFT).GT.I-1 ) */
226 
227  back = werr[ii];
228 L20:
229  negcnt = template_lapack_laneg(n, &d__[1], &lld[1], &left, pivmin, &r__);
230  if (negcnt > i__ - 1) {
231  left -= back;
232  back *= 2.;
233  goto L20;
234  }
235 
236 /* Do while( NEGCNT(RIGHT).LT.I ) */
237 /* Compute negcount from dstqds facto L+D+L+^T = L D L^T - RIGHT */
238 
239  back = werr[ii];
240 L50:
241  negcnt = template_lapack_laneg(n, &d__[1], &lld[1], &right, pivmin, &r__);
242  if (negcnt < i__) {
243  right += back;
244  back *= 2.;
245  goto L50;
246  }
247  width = (d__1 = left - right, absMACRO(d__1)) * .5;
248 /* Computing MAX */
249  d__1 = absMACRO(left), d__2 = absMACRO(right);
250  tmp = maxMACRO(d__1,d__2);
251 /* Computing MAX */
252  d__1 = *rtol1 * gap, d__2 = *rtol2 * tmp;
253  cvrgd = maxMACRO(d__1,d__2);
254  if (width <= cvrgd || width <= mnwdth) {
255 /* This interval has already converged and does not need refinement. */
256 /* (Note that the gaps might change through refining the */
257 /* eigenvalues, however, they can only get bigger.) */
258 /* Remove it from the list. */
259  iwork[k - 1] = -1;
260 /* Make sure that I1 always points to the first unconverged interval */
261  if (i__ == i1 && i__ < *ilast) {
262  i1 = i__ + 1;
263  }
264  if (prev >= i1 && i__ <= *ilast) {
265  iwork[(prev << 1) - 1] = i__ + 1;
266  }
267  } else {
268 /* unconverged interval found */
269  prev = i__;
270  ++nint;
271  iwork[k - 1] = i__ + 1;
272  iwork[k] = negcnt;
273  }
274  work[k - 1] = left;
275  work[k] = right;
276 /* L75: */
277  }
278 
279 /* Do while( NINT.GT.0 ), i.e. there are still unconverged intervals */
280 /* and while (ITER.LT.MAXITR) */
281 
282  iter = 0;
283 L80:
284  prev = i1 - 1;
285  i__ = i1;
286  olnint = nint;
287  i__1 = olnint;
288  for (ip = 1; ip <= i__1; ++ip) {
289  k = i__ << 1;
290  ii = i__ - *offset;
291  rgap = wgap[ii];
292  lgap = rgap;
293  if (ii > 1) {
294  lgap = wgap[ii - 1];
295  }
296  gap = minMACRO(lgap,rgap);
297  next = iwork[k - 1];
298  left = work[k - 1];
299  right = work[k];
300  mid = (left + right) * .5;
301 /* semiwidth of interval */
302  width = right - mid;
303 /* Computing MAX */
304  d__1 = absMACRO(left), d__2 = absMACRO(right);
305  tmp = maxMACRO(d__1,d__2);
306 /* Computing MAX */
307  d__1 = *rtol1 * gap, d__2 = *rtol2 * tmp;
308  cvrgd = maxMACRO(d__1,d__2);
309  if (width <= cvrgd || width <= mnwdth || iter == maxitr) {
310 /* reduce number of unconverged intervals */
311  --nint;
312 /* Mark interval as converged. */
313  iwork[k - 1] = 0;
314  if (i1 == i__) {
315  i1 = next;
316  } else {
317 /* Prev holds the last unconverged interval previously examined */
318  if (prev >= i1) {
319  iwork[(prev << 1) - 1] = next;
320  }
321  }
322  i__ = next;
323  goto L100;
324  }
325  prev = i__;
326 
327 /* Perform one bisection step */
328 
329  negcnt = template_lapack_laneg(n, &d__[1], &lld[1], &mid, pivmin, &r__);
330  if (negcnt <= i__ - 1) {
331  work[k - 1] = mid;
332  } else {
333  work[k] = mid;
334  }
335  i__ = next;
336 L100:
337  ;
338  }
339  ++iter;
340 /* do another loop if there are still unconverged intervals */
341 /* However, in the last iteration, all intervals are accepted */
342 /* since this is the best we can do. */
343  if (nint > 0 && iter <= maxitr) {
344  goto L80;
345  }
346 
347 
348 /* At this point, all the intervals have converged */
349  i__1 = *ilast;
350  for (i__ = *ifirst; i__ <= i__1; ++i__) {
351  k = i__ << 1;
352  ii = i__ - *offset;
353 /* All intervals marked by '0' have been refined. */
354  if (iwork[k - 1] == 0) {
355  w[ii] = (work[k - 1] + work[k]) * .5;
356  werr[ii] = work[k] - w[ii];
357  }
358 /* L110: */
359  }
360 
361  i__1 = *ilast;
362  for (i__ = *ifirst + 1; i__ <= i__1; ++i__) {
363  k = i__ << 1;
364  ii = i__ - *offset;
365 /* Computing MAX */
366  d__1 = 0., d__2 = w[ii] - werr[ii] - w[ii - 1] - werr[ii - 1];
367  wgap[ii - 1] = maxMACRO(d__1,d__2);
368 /* L111: */
369  }
370  return 0;
371 
372 /* End of DLARRB */
373 
374 } /* dlarrb_ */
375 
376 #endif
minMACRO
#define minMACRO(a, b)
Definition: template_blas_common.h:46
absMACRO
#define absMACRO(x)
Definition: template_blas_common.h:47
mat::left
@ left
Definition: Matrix.h:75
template_lapack_larrb
int template_lapack_larrb(integer *n, Treal *d__, Treal *lld, integer *ifirst, integer *ilast, Treal *rtol1, Treal *rtol2, integer *offset, Treal *w, Treal *wgap, Treal *werr, Treal *work, integer *iwork, Treal *pivmin, Treal *spdiam, integer *twist, integer *info)
Definition: template_lapack_larrb.h:45
mat::right
@ right
Definition: Matrix.h:75
template_blas_log
Treal template_blas_log(Treal x)
template_lapack_laneg.h
integer
int integer
Definition: template_blas_common.h:40
template_lapack_laneg
int template_lapack_laneg(integer *n, Treal *d__, Treal *lld, Treal *sigma, Treal *pivmin, integer *r__)
Definition: template_lapack_laneg.h:45
maxMACRO
#define maxMACRO(a, b)
Definition: template_blas_common.h:45