ergo
template_lapack_larrd.h
Go to the documentation of this file.
1 /* Ergo, version 3.4, a program for linear scaling electronic structure
2  * calculations.
3  * Copyright (C) 2014 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_LARRD_HEADER
36 #define TEMPLATE_LAPACK_LARRD_HEADER
37 
38 template<class Treal>
39 int template_lapack_larrd(const char *range, const char *order, const integer *n, Treal
40  *vl, Treal *vu, integer *il, integer *iu, Treal *gers,
41  Treal *reltol, Treal *d__, Treal *e, Treal *e2,
42  Treal *pivmin, integer *nsplit, integer *isplit, integer *m,
43  Treal *w, Treal *werr, Treal *wl, Treal *wu,
44  integer *iblock, integer *indexw, Treal *work, integer *iwork,
45  integer *info)
46 {
47  /* System generated locals */
48  integer i__1, i__2, i__3;
49  Treal d__1, d__2;
50 
51 
52  /* Local variables */
53  integer i__, j, ib, ie, je, nb;
54  Treal gl;
55  integer im, in;
56  Treal gu;
57  integer iw, jee;
58  Treal eps;
59  integer nwl;
60  Treal wlu = 0;
61  Treal wul = 0;
62  integer nwu;
63  Treal tmp1, tmp2;
64  integer iend, jblk, ioff, iout, itmp1, itmp2, jdisc;
65  integer iinfo;
66  Treal atoli;
67  integer iwoff, itmax;
68  Treal wkill, rtoli, uflow, tnorm;
69  integer ibegin;
70  integer irange, idiscl, idumma[1];
71  integer idiscu;
72  logical ncnvrg, toofew;
73 
74 
75 /* -- LAPACK auxiliary routine (version 3.2.1) -- */
76 /* -- LAPACK is a software package provided by Univ. of Tennessee, -- */
77 /* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- */
78 /* -- April 2009 -- */
79 
80 /* .. Scalar Arguments .. */
81 /* .. */
82 /* .. Array Arguments .. */
83 /* .. */
84 
85 /* Purpose */
86 /* ======= */
87 
88 /* DLARRD computes the eigenvalues of a symmetric tridiagonal */
89 /* matrix T to suitable accuracy. This is an auxiliary code to be */
90 /* called from DSTEMR. */
91 /* The user may ask for all eigenvalues, all eigenvalues */
92 /* in the half-open interval (VL, VU], or the IL-th through IU-th */
93 /* eigenvalues. */
94 
95 /* To avoid overflow, the matrix must be scaled so that its */
96 /* largest element is no greater than overflow**(1/2) * */
97 /* underflow**(1/4) in absolute value, and for greatest */
98 /* accuracy, it should not be much smaller than that. */
99 
100 /* See W. Kahan "Accurate Eigenvalues of a Symmetric Tridiagonal */
101 /* Matrix", Report CS41, Computer Science Dept., Stanford */
102 /* University, July 21, 1966. */
103 
104 /* Arguments */
105 /* ========= */
106 
107 /* RANGE (input) CHARACTER */
108 /* = 'A': ("All") all eigenvalues will be found. */
109 /* = 'V': ("Value") all eigenvalues in the half-open interval */
110 /* (VL, VU] will be found. */
111 /* = 'I': ("Index") the IL-th through IU-th eigenvalues (of the */
112 /* entire matrix) will be found. */
113 
114 /* ORDER (input) CHARACTER */
115 /* = 'B': ("By Block") the eigenvalues will be grouped by */
116 /* split-off block (see IBLOCK, ISPLIT) and */
117 /* ordered from smallest to largest within */
118 /* the block. */
119 /* = 'E': ("Entire matrix") */
120 /* the eigenvalues for the entire matrix */
121 /* will be ordered from smallest to */
122 /* largest. */
123 
124 /* N (input) INTEGER */
125 /* The order of the tridiagonal matrix T. N >= 0. */
126 
127 /* VL (input) DOUBLE PRECISION */
128 /* VU (input) DOUBLE PRECISION */
129 /* If RANGE='V', the lower and upper bounds of the interval to */
130 /* be searched for eigenvalues. Eigenvalues less than or equal */
131 /* to VL, or greater than VU, will not be returned. VL < VU. */
132 /* Not referenced if RANGE = 'A' or 'I'. */
133 
134 /* IL (input) INTEGER */
135 /* IU (input) INTEGER */
136 /* If RANGE='I', the indices (in ascending order) of the */
137 /* smallest and largest eigenvalues to be returned. */
138 /* 1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0. */
139 /* Not referenced if RANGE = 'A' or 'V'. */
140 
141 /* GERS (input) DOUBLE PRECISION array, dimension (2*N) */
142 /* The N Gerschgorin intervals (the i-th Gerschgorin interval */
143 /* is (GERS(2*i-1), GERS(2*i)). */
144 
145 /* RELTOL (input) DOUBLE PRECISION */
146 /* The minimum relative width of an interval. When an interval */
147 /* is narrower than RELTOL times the larger (in */
148 /* magnitude) endpoint, then it is considered to be */
149 /* sufficiently small, i.e., converged. Note: this should */
150 /* always be at least radix*machine epsilon. */
151 
152 /* D (input) DOUBLE PRECISION array, dimension (N) */
153 /* The n diagonal elements of the tridiagonal matrix T. */
154 
155 /* E (input) DOUBLE PRECISION array, dimension (N-1) */
156 /* The (n-1) off-diagonal elements of the tridiagonal matrix T. */
157 
158 /* E2 (input) DOUBLE PRECISION array, dimension (N-1) */
159 /* The (n-1) squared off-diagonal elements of the tridiagonal matrix T. */
160 
161 /* PIVMIN (input) DOUBLE PRECISION */
162 /* The minimum pivot allowed in the Sturm sequence for T. */
163 
164 /* NSPLIT (input) INTEGER */
165 /* The number of diagonal blocks in the matrix T. */
166 /* 1 <= NSPLIT <= N. */
167 
168 /* ISPLIT (input) INTEGER array, dimension (N) */
169 /* The splitting points, at which T breaks up into submatrices. */
170 /* The first submatrix consists of rows/columns 1 to ISPLIT(1), */
171 /* the second of rows/columns ISPLIT(1)+1 through ISPLIT(2), */
172 /* etc., and the NSPLIT-th consists of rows/columns */
173 /* ISPLIT(NSPLIT-1)+1 through ISPLIT(NSPLIT)=N. */
174 /* (Only the first NSPLIT elements will actually be used, but */
175 /* since the user cannot know a priori what value NSPLIT will */
176 /* have, N words must be reserved for ISPLIT.) */
177 
178 /* M (output) INTEGER */
179 /* The actual number of eigenvalues found. 0 <= M <= N. */
180 /* (See also the description of INFO=2,3.) */
181 
182 /* W (output) DOUBLE PRECISION array, dimension (N) */
183 /* On exit, the first M elements of W will contain the */
184 /* eigenvalue approximations. DLARRD computes an interval */
185 /* I_j = (a_j, b_j] that includes eigenvalue j. The eigenvalue */
186 /* approximation is given as the interval midpoint */
187 /* W(j)= ( a_j + b_j)/2. The corresponding error is bounded by */
188 /* WERR(j) = abs( a_j - b_j)/2 */
189 
190 /* WERR (output) DOUBLE PRECISION array, dimension (N) */
191 /* The error bound on the corresponding eigenvalue approximation */
192 /* in W. */
193 
194 /* WL (output) DOUBLE PRECISION */
195 /* WU (output) DOUBLE PRECISION */
196 /* The interval (WL, WU] contains all the wanted eigenvalues. */
197 /* If RANGE='V', then WL=VL and WU=VU. */
198 /* If RANGE='A', then WL and WU are the global Gerschgorin bounds */
199 /* on the spectrum. */
200 /* If RANGE='I', then WL and WU are computed by DLAEBZ from the */
201 /* index range specified. */
202 
203 /* IBLOCK (output) INTEGER array, dimension (N) */
204 /* At each row/column j where E(j) is zero or small, the */
205 /* matrix T is considered to split into a block diagonal */
206 /* matrix. On exit, if INFO = 0, IBLOCK(i) specifies to which */
207 /* block (from 1 to the number of blocks) the eigenvalue W(i) */
208 /* belongs. (DLARRD may use the remaining N-M elements as */
209 /* workspace.) */
210 
211 /* INDEXW (output) INTEGER array, dimension (N) */
212 /* The indices of the eigenvalues within each block (submatrix); */
213 /* for example, INDEXW(i)= j and IBLOCK(i)=k imply that the */
214 /* i-th eigenvalue W(i) is the j-th eigenvalue in block k. */
215 
216 /* WORK (workspace) DOUBLE PRECISION array, dimension (4*N) */
217 
218 /* IWORK (workspace) INTEGER array, dimension (3*N) */
219 
220 /* INFO (output) INTEGER */
221 /* = 0: successful exit */
222 /* < 0: if INFO = -i, the i-th argument had an illegal value */
223 /* > 0: some or all of the eigenvalues failed to converge or */
224 /* were not computed: */
225 /* =1 or 3: Bisection failed to converge for some */
226 /* eigenvalues; these eigenvalues are flagged by a */
227 /* negative block number. The effect is that the */
228 /* eigenvalues may not be as accurate as the */
229 /* absolute and relative tolerances. This is */
230 /* generally caused by unexpectedly inaccurate */
231 /* arithmetic. */
232 /* =2 or 3: RANGE='I' only: Not all of the eigenvalues */
233 /* IL:IU were found. */
234 /* Effect: M < IU+1-IL */
235 /* Cause: non-monotonic arithmetic, causing the */
236 /* Sturm sequence to be non-monotonic. */
237 /* Cure: recalculate, using RANGE='A', and pick */
238 /* out eigenvalues IL:IU. In some cases, */
239 /* increasing the PARAMETER "FUDGE" may */
240 /* make things work. */
241 /* = 4: RANGE='I', and the Gershgorin interval */
242 /* initially used was too small. No eigenvalues */
243 /* were computed. */
244 /* Probable cause: your machine has sloppy */
245 /* floating-point arithmetic. */
246 /* Cure: Increase the PARAMETER "FUDGE", */
247 /* recompile, and try again. */
248 
249 /* Internal Parameters */
250 /* =================== */
251 
252 /* FUDGE DOUBLE PRECISION, default = 2 */
253 /* A "fudge factor" to widen the Gershgorin intervals. Ideally, */
254 /* a value of 1 should work, but on machines with sloppy */
255 /* arithmetic, this needs to be larger. The default for */
256 /* publicly released versions should be large enough to handle */
257 /* the worst machine around. Note that this has no effect */
258 /* on accuracy of the solution. */
259 
260 /* Based on contributions by */
261 /* W. Kahan, University of California, Berkeley, USA */
262 /* Beresford Parlett, University of California, Berkeley, USA */
263 /* Jim Demmel, University of California, Berkeley, USA */
264 /* Inderjit Dhillon, University of Texas, Austin, USA */
265 /* Osni Marques, LBNL/NERSC, USA */
266 /* Christof Voemel, University of California, Berkeley, USA */
267 
268 /* ===================================================================== */
269 
270 /* .. Parameters .. */
271 /* .. */
272 /* .. Local Scalars .. */
273 /* .. */
274 /* .. Local Arrays .. */
275 /* .. */
276 /* .. External Functions .. */
277 /* .. */
278 /* .. External Subroutines .. */
279 /* .. */
280 /* .. Intrinsic Functions .. */
281 /* .. */
282 /* .. Executable Statements .. */
283 
284 /* Table of constant values */
285 
286  integer c__1 = 1;
287  integer c_n1 = -1;
288  integer c__3 = 3;
289  integer c__2 = 2;
290  integer c__0 = 0;
291 
292  /* Parameter adjustments */
293  --iwork;
294  --work;
295  --indexw;
296  --iblock;
297  --werr;
298  --w;
299  --isplit;
300  --e2;
301  --e;
302  --d__;
303  --gers;
304 
305  /* Function Body */
306  *info = 0;
307 
308 /* Decode RANGE */
309 
310  if (template_blas_lsame(range, "A")) {
311  irange = 1;
312  } else if (template_blas_lsame(range, "V")) {
313  irange = 2;
314  } else if (template_blas_lsame(range, "I")) {
315  irange = 3;
316  } else {
317  irange = 0;
318  }
319 
320 /* Check for Errors */
321 
322  if (irange <= 0) {
323  *info = -1;
324  } else if (! (template_blas_lsame(order, "B") || template_blas_lsame(order,
325  "E"))) {
326  *info = -2;
327  } else if (*n < 0) {
328  *info = -3;
329  } else if (irange == 2) {
330  if (*vl >= *vu) {
331  *info = -5;
332  }
333  } else if (irange == 3 && (*il < 1 || *il > maxMACRO(1,*n))) {
334  *info = -6;
335  } else if (irange == 3 && (*iu < minMACRO(*n,*il) || *iu > *n)) {
336  *info = -7;
337  }
338 
339  if (*info != 0) {
340  return 0;
341  }
342 /* Initialize error flags */
343  *info = 0;
344  ncnvrg = FALSE_;
345  toofew = FALSE_;
346 /* Quick return if possible */
347  *m = 0;
348  if (*n == 0) {
349  return 0;
350  }
351 /* Simplification: */
352  if (irange == 3 && *il == 1 && *iu == *n) {
353  irange = 1;
354  }
355 /* Get machine constants */
356  eps = template_lapack_lamch("P",(Treal)0);
357  uflow = template_lapack_lamch("U",(Treal)0);
358 /* Special Case when N=1 */
359 /* Treat case of 1x1 matrix for quick return */
360  if (*n == 1) {
361  if ( irange == 1 || ( irange == 2 && d__[1] > *vl && d__[1] <= *vu ) ||
362  ( irange == 3 && *il == 1 && *iu == 1 ) ) {
363  *m = 1;
364  w[1] = d__[1];
365 /* The computation error of the eigenvalue is zero */
366  werr[1] = 0.;
367  iblock[1] = 1;
368  indexw[1] = 1;
369  }
370  return 0;
371  }
372 /* NB is the minimum vector length for vector bisection, or 0 */
373 /* if only scalar is to be done. */
374  nb = template_lapack_ilaenv(&c__1, "DSTEBZ", " ", n, &c_n1, &c_n1, &c_n1, (ftnlen)6, (ftnlen)1);
375  if (nb <= 1) {
376  nb = 0;
377  }
378 /* Find global spectral radius */
379  gl = d__[1];
380  gu = d__[1];
381  i__1 = *n;
382  for (i__ = 1; i__ <= i__1; ++i__) {
383 /* Computing MIN */
384  d__1 = gl, d__2 = gers[(i__ << 1) - 1];
385  gl = minMACRO(d__1,d__2);
386 /* Computing MAX */
387  d__1 = gu, d__2 = gers[i__ * 2];
388  gu = maxMACRO(d__1,d__2);
389 /* L5: */
390  }
391 /* Compute global Gerschgorin bounds and spectral diameter */
392 /* Computing MAX */
393  d__1 = absMACRO(gl), d__2 = absMACRO(gu);
394  tnorm = maxMACRO(d__1,d__2);
395  gl = gl - tnorm * 2. * eps * *n - *pivmin * 4.;
396  gu = gu + tnorm * 2. * eps * *n + *pivmin * 4.;
397 /* [JAN/28/2009] remove the line below since SPDIAM variable not use */
398 /* SPDIAM = GU - GL */
399 /* Input arguments for DLAEBZ: */
400 /* The relative tolerance. An interval (a,b] lies within */
401 /* "relative tolerance" if b-a < RELTOL*max(|a|,|b|), */
402  rtoli = *reltol;
403 /* Set the absolute tolerance for interval convergence to zero to force */
404 /* interval convergence based on relative size of the interval. */
405 /* This is dangerous because intervals might not converge when RELTOL is */
406 /* small. But at least a very small number should be selected so that for */
407 /* strongly graded matrices, the code can get relatively accurate */
408 /* eigenvalues. */
409  atoli = uflow * 4. + *pivmin * 4.;
410  if (irange == 3) {
411 /* RANGE='I': Compute an interval containing eigenvalues */
412 /* IL through IU. The initial interval [GL,GU] from the global */
413 /* Gerschgorin bounds GL and GU is refined by DLAEBZ. */
414  itmax = (integer) ((template_blas_log(tnorm + *pivmin) - template_blas_log(*pivmin)) / template_blas_log(2.)) +
415  2;
416  work[*n + 1] = gl;
417  work[*n + 2] = gl;
418  work[*n + 3] = gu;
419  work[*n + 4] = gu;
420  work[*n + 5] = gl;
421  work[*n + 6] = gu;
422  iwork[1] = -1;
423  iwork[2] = -1;
424  iwork[3] = *n + 1;
425  iwork[4] = *n + 1;
426  iwork[5] = *il - 1;
427  iwork[6] = *iu;
428 
429  template_lapack_laebz(&c__3, &itmax, n, &c__2, &c__2, &nb, &atoli, &rtoli, pivmin, &
430  d__[1], &e[1], &e2[1], &iwork[5], &work[*n + 1], &work[*n + 5]
431 , &iout, &iwork[1], &w[1], &iblock[1], &iinfo);
432  if (iinfo != 0) {
433  *info = iinfo;
434  return 0;
435  }
436 /* On exit, output intervals may not be ordered by ascending negcount */
437  if (iwork[6] == *iu) {
438  *wl = work[*n + 1];
439  wlu = work[*n + 3];
440  nwl = iwork[1];
441  *wu = work[*n + 4];
442  wul = work[*n + 2];
443  nwu = iwork[4];
444  } else {
445  *wl = work[*n + 2];
446  wlu = work[*n + 4];
447  nwl = iwork[2];
448  *wu = work[*n + 3];
449  wul = work[*n + 1];
450  nwu = iwork[3];
451  }
452 /* On exit, the interval [WL, WLU] contains a value with negcount NWL, */
453 /* and [WUL, WU] contains a value with negcount NWU. */
454  if (nwl < 0 || nwl >= *n || nwu < 1 || nwu > *n) {
455  *info = 4;
456  return 0;
457  }
458  } else if (irange == 2) {
459  *wl = *vl;
460  *wu = *vu;
461  } else if (irange == 1) {
462  *wl = gl;
463  *wu = gu;
464  }
465 /* Find Eigenvalues -- Loop Over blocks and recompute NWL and NWU. */
466 /* NWL accumulates the number of eigenvalues .le. WL, */
467 /* NWU accumulates the number of eigenvalues .le. WU */
468  *m = 0;
469  iend = 0;
470  *info = 0;
471  nwl = 0;
472  nwu = 0;
473 
474  i__1 = *nsplit;
475  for (jblk = 1; jblk <= i__1; ++jblk) {
476  ioff = iend;
477  ibegin = ioff + 1;
478  iend = isplit[jblk];
479  in = iend - ioff;
480 
481  if (in == 1) {
482 /* 1x1 block */
483  if (*wl >= d__[ibegin] - *pivmin) {
484  ++nwl;
485  }
486  if (*wu >= d__[ibegin] - *pivmin) {
487  ++nwu;
488  }
489  if (irange == 1 || ( *wl < d__[ibegin] - *pivmin && *wu >= d__[
490  ibegin] - *pivmin ) ) {
491  ++(*m);
492  w[*m] = d__[ibegin];
493  werr[*m] = 0.;
494 /* The gap for a single block doesn't matter for the later */
495 /* algorithm and is assigned an arbitrary large value */
496  iblock[*m] = jblk;
497  indexw[*m] = 1;
498  }
499 /* Disabled 2x2 case because of a failure on the following matrix */
500 /* RANGE = 'I', IL = IU = 4 */
501 /* Original Tridiagonal, d = [ */
502 /* -0.150102010615740E+00 */
503 /* -0.849897989384260E+00 */
504 /* -0.128208148052635E-15 */
505 /* 0.128257718286320E-15 */
506 /* ]; */
507 /* e = [ */
508 /* -0.357171383266986E+00 */
509 /* -0.180411241501588E-15 */
510 /* -0.175152352710251E-15 */
511 /* ]; */
512 
513 /* ELSE IF( IN.EQ.2 ) THEN */
514 /* * 2x2 block */
515 /* DISC = SQRT( (HALF*(D(IBEGIN)-D(IEND)))**2 + E(IBEGIN)**2 ) */
516 /* TMP1 = HALF*(D(IBEGIN)+D(IEND)) */
517 /* L1 = TMP1 - DISC */
518 /* IF( WL.GE. L1-PIVMIN ) */
519 /* $ NWL = NWL + 1 */
520 /* IF( WU.GE. L1-PIVMIN ) */
521 /* $ NWU = NWU + 1 */
522 /* IF( IRANGE.EQ.ALLRNG .OR. ( WL.LT.L1-PIVMIN .AND. WU.GE. */
523 /* $ L1-PIVMIN ) ) THEN */
524 /* M = M + 1 */
525 /* W( M ) = L1 */
526 /* * The uncertainty of eigenvalues of a 2x2 matrix is very small */
527 /* WERR( M ) = EPS * ABS( W( M ) ) * TWO */
528 /* IBLOCK( M ) = JBLK */
529 /* INDEXW( M ) = 1 */
530 /* ENDIF */
531 /* L2 = TMP1 + DISC */
532 /* IF( WL.GE. L2-PIVMIN ) */
533 /* $ NWL = NWL + 1 */
534 /* IF( WU.GE. L2-PIVMIN ) */
535 /* $ NWU = NWU + 1 */
536 /* IF( IRANGE.EQ.ALLRNG .OR. ( WL.LT.L2-PIVMIN .AND. WU.GE. */
537 /* $ L2-PIVMIN ) ) THEN */
538 /* M = M + 1 */
539 /* W( M ) = L2 */
540 /* * The uncertainty of eigenvalues of a 2x2 matrix is very small */
541 /* WERR( M ) = EPS * ABS( W( M ) ) * TWO */
542 /* IBLOCK( M ) = JBLK */
543 /* INDEXW( M ) = 2 */
544 /* ENDIF */
545  } else {
546 /* General Case - block of size IN >= 2 */
547 /* Compute local Gerschgorin interval and use it as the initial */
548 /* interval for DLAEBZ */
549  gu = d__[ibegin];
550  gl = d__[ibegin];
551  tmp1 = 0.;
552  i__2 = iend;
553  for (j = ibegin; j <= i__2; ++j) {
554 /* Computing MIN */
555  d__1 = gl, d__2 = gers[(j << 1) - 1];
556  gl = minMACRO(d__1,d__2);
557 /* Computing MAX */
558  d__1 = gu, d__2 = gers[j * 2];
559  gu = maxMACRO(d__1,d__2);
560 /* L40: */
561  }
562 /* [JAN/28/2009] */
563 /* change SPDIAM by TNORM in lines 2 and 3 thereafter */
564 /* line 1: remove computation of SPDIAM (not useful anymore) */
565 /* SPDIAM = GU - GL */
566 /* GL = GL - FUDGE*SPDIAM*EPS*IN - FUDGE*PIVMIN */
567 /* GU = GU + FUDGE*SPDIAM*EPS*IN + FUDGE*PIVMIN */
568  gl = gl - tnorm * 2. * eps * in - *pivmin * 2.;
569  gu = gu + tnorm * 2. * eps * in + *pivmin * 2.;
570 
571  if (irange > 1) {
572  if (gu < *wl) {
573 /* the local block contains none of the wanted eigenvalues */
574  nwl += in;
575  nwu += in;
576  goto L70;
577  }
578 /* refine search interval if possible, only range (WL,WU] matters */
579  gl = maxMACRO(gl,*wl);
580  gu = minMACRO(gu,*wu);
581  if (gl >= gu) {
582  goto L70;
583  }
584  }
585 /* Find negcount of initial interval boundaries GL and GU */
586  work[*n + 1] = gl;
587  work[*n + in + 1] = gu;
588  template_lapack_laebz(&c__1, &c__0, &in, &in, &c__1, &nb, &atoli, &rtoli,
589  pivmin, &d__[ibegin], &e[ibegin], &e2[ibegin], idumma, &
590  work[*n + 1], &work[*n + (in << 1) + 1], &im, &iwork[1], &
591  w[*m + 1], &iblock[*m + 1], &iinfo);
592  if (iinfo != 0) {
593  *info = iinfo;
594  return 0;
595  }
596 
597  nwl += iwork[1];
598  nwu += iwork[in + 1];
599  iwoff = *m - iwork[1];
600 /* Compute Eigenvalues */
601  itmax = (integer) ((template_blas_log(gu - gl + *pivmin) - template_blas_log(*pivmin)) / template_blas_log(
602  2.)) + 2;
603  template_lapack_laebz(&c__2, &itmax, &in, &in, &c__1, &nb, &atoli, &rtoli,
604  pivmin, &d__[ibegin], &e[ibegin], &e2[ibegin], idumma, &
605  work[*n + 1], &work[*n + (in << 1) + 1], &iout, &iwork[1],
606  &w[*m + 1], &iblock[*m + 1], &iinfo);
607  if (iinfo != 0) {
608  *info = iinfo;
609  return 0;
610  }
611 
612 /* Copy eigenvalues into W and IBLOCK */
613 /* Use -JBLK for block number for unconverged eigenvalues. */
614 /* Loop over the number of output intervals from DLAEBZ */
615  i__2 = iout;
616  for (j = 1; j <= i__2; ++j) {
617 /* eigenvalue approximation is middle point of interval */
618  tmp1 = (work[j + *n] + work[j + in + *n]) * .5;
619 /* semi length of error interval */
620  tmp2 = (d__1 = work[j + *n] - work[j + in + *n], absMACRO(d__1)) *
621  .5;
622  if (j > iout - iinfo) {
623 /* Flag non-convergence. */
624  ncnvrg = TRUE_;
625  ib = -jblk;
626  } else {
627  ib = jblk;
628  }
629  i__3 = iwork[j + in] + iwoff;
630  for (je = iwork[j] + 1 + iwoff; je <= i__3; ++je) {
631  w[je] = tmp1;
632  werr[je] = tmp2;
633  indexw[je] = je - iwoff;
634  iblock[je] = ib;
635 /* L50: */
636  }
637 /* L60: */
638  }
639 
640  *m += im;
641  }
642 L70:
643  ;
644  }
645 /* If RANGE='I', then (WL,WU) contains eigenvalues NWL+1,...,NWU */
646 /* If NWL+1 < IL or NWU > IU, discard extra eigenvalues. */
647  if (irange == 3) {
648  idiscl = *il - 1 - nwl;
649  idiscu = nwu - *iu;
650 
651  if (idiscl > 0) {
652  im = 0;
653  i__1 = *m;
654  for (je = 1; je <= i__1; ++je) {
655 /* Remove some of the smallest eigenvalues from the left so that */
656 /* at the end IDISCL =0. Move all eigenvalues up to the left. */
657  if (w[je] <= wlu && idiscl > 0) {
658  --idiscl;
659  } else {
660  ++im;
661  w[im] = w[je];
662  werr[im] = werr[je];
663  indexw[im] = indexw[je];
664  iblock[im] = iblock[je];
665  }
666 /* L80: */
667  }
668  *m = im;
669  }
670  if (idiscu > 0) {
671 /* Remove some of the largest eigenvalues from the right so that */
672 /* at the end IDISCU =0. Move all eigenvalues up to the left. */
673  im = *m + 1;
674  for (je = *m; je >= 1; --je) {
675  if (w[je] >= wul && idiscu > 0) {
676  --idiscu;
677  } else {
678  --im;
679  w[im] = w[je];
680  werr[im] = werr[je];
681  indexw[im] = indexw[je];
682  iblock[im] = iblock[je];
683  }
684 /* L81: */
685  }
686  jee = 0;
687  i__1 = *m;
688  for (je = im; je <= i__1; ++je) {
689  ++jee;
690  w[jee] = w[je];
691  werr[jee] = werr[je];
692  indexw[jee] = indexw[je];
693  iblock[jee] = iblock[je];
694 /* L82: */
695  }
696  *m = *m - im + 1;
697  }
698  if (idiscl > 0 || idiscu > 0) {
699 /* Code to deal with effects of bad arithmetic. (If N(w) is */
700 /* monotone non-decreasing, this should never happen.) */
701 /* Some low eigenvalues to be discarded are not in (WL,WLU], */
702 /* or high eigenvalues to be discarded are not in (WUL,WU] */
703 /* so just kill off the smallest IDISCL/largest IDISCU */
704 /* eigenvalues, by marking the corresponding IBLOCK = 0 */
705  if (idiscl > 0) {
706  wkill = *wu;
707  i__1 = idiscl;
708  for (jdisc = 1; jdisc <= i__1; ++jdisc) {
709  iw = 0;
710  i__2 = *m;
711  for (je = 1; je <= i__2; ++je) {
712  if (iblock[je] != 0 && (w[je] < wkill || iw == 0)) {
713  iw = je;
714  wkill = w[je];
715  }
716 /* L90: */
717  }
718  iblock[iw] = 0;
719 /* L100: */
720  }
721  }
722  if (idiscu > 0) {
723  wkill = *wl;
724  i__1 = idiscu;
725  for (jdisc = 1; jdisc <= i__1; ++jdisc) {
726  iw = 0;
727  i__2 = *m;
728  for (je = 1; je <= i__2; ++je) {
729  if (iblock[je] != 0 && (w[je] >= wkill || iw == 0)) {
730  iw = je;
731  wkill = w[je];
732  }
733 /* L110: */
734  }
735  iblock[iw] = 0;
736 /* L120: */
737  }
738  }
739 /* Now erase all eigenvalues with IBLOCK set to zero */
740  im = 0;
741  i__1 = *m;
742  for (je = 1; je <= i__1; ++je) {
743  if (iblock[je] != 0) {
744  ++im;
745  w[im] = w[je];
746  werr[im] = werr[je];
747  indexw[im] = indexw[je];
748  iblock[im] = iblock[je];
749  }
750 /* L130: */
751  }
752  *m = im;
753  }
754  if (idiscl < 0 || idiscu < 0) {
755  toofew = TRUE_;
756  }
757  }
758 
759  if ( ( irange == 1 && *m != *n ) || ( irange == 3 && *m != *iu - *il + 1 ) ) {
760  toofew = TRUE_;
761  }
762 /* If ORDER='B', do nothing the eigenvalues are already sorted by */
763 /* block. */
764 /* If ORDER='E', sort the eigenvalues from smallest to largest */
765  if (template_blas_lsame(order, "E") && *nsplit > 1) {
766  i__1 = *m - 1;
767  for (je = 1; je <= i__1; ++je) {
768  ie = 0;
769  tmp1 = w[je];
770  i__2 = *m;
771  for (j = je + 1; j <= i__2; ++j) {
772  if (w[j] < tmp1) {
773  ie = j;
774  tmp1 = w[j];
775  }
776 /* L140: */
777  }
778  if (ie != 0) {
779  tmp2 = werr[ie];
780  itmp1 = iblock[ie];
781  itmp2 = indexw[ie];
782  w[ie] = w[je];
783  werr[ie] = werr[je];
784  iblock[ie] = iblock[je];
785  indexw[ie] = indexw[je];
786  w[je] = tmp1;
787  werr[je] = tmp2;
788  iblock[je] = itmp1;
789  indexw[je] = itmp2;
790  }
791 /* L150: */
792  }
793  }
794 
795  *info = 0;
796  if (ncnvrg) {
797  ++(*info);
798  }
799  if (toofew) {
800  *info += 2;
801  }
802  return 0;
803 
804 /* End of DLARRD */
805 
806 } /* dlarrd_ */
807 
808 #endif
static const real gu
Definition: fun-pz81.c:71
#define absMACRO(x)
Definition: template_blas_common.h:45
int template_lapack_laebz(const integer *ijob, const integer *nitmax, const integer *n, const integer *mmax, const integer *minp, const integer *nbmin, const Treal *abstol, const Treal *reltol, const Treal *pivmin, const Treal *d__, const Treal *e, const Treal *e2, integer *nval, Treal *ab, Treal *c__, integer *mout, integer *nab, Treal *work, integer *iwork, integer *info)
Definition: template_lapack_laebz.h:40
int integer
Definition: template_blas_common.h:38
integer template_lapack_ilaenv(const integer *ispec, const char *name__, const char *opts, const integer *n1, const integer *n2, const integer *n3, const integer *n4, ftnlen name_len, ftnlen opts_len)
Definition: template_lapack_common.cc:279
#define maxMACRO(a, b)
Definition: template_blas_common.h:43
#define minMACRO(a, b)
Definition: template_blas_common.h:44
Treal template_blas_log(Treal x)
int template_lapack_larrd(const char *range, const char *order, const integer *n, Treal *vl, Treal *vu, integer *il, integer *iu, Treal *gers, Treal *reltol, Treal *d__, Treal *e, Treal *e2, Treal *pivmin, integer *nsplit, integer *isplit, integer *m, Treal *w, Treal *werr, Treal *wl, Treal *wu, integer *iblock, integer *indexw, Treal *work, integer *iwork, integer *info)
Definition: template_lapack_larrd.h:39
Treal template_lapack_lamch(const char *cmach, Treal dummyReal)
Definition: template_lapack_lamch.h:199
bool logical
Definition: template_blas_common.h:39
#define TRUE_
Definition: template_lapack_common.h:40
#define FALSE_
Definition: template_lapack_common.h:41
int ftnlen
Definition: template_blas_common.h:40
logical template_blas_lsame(const char *ca, const char *cb)
Definition: template_blas_common.cc:44