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