ergo
template_lapack_stebz.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_STEBZ_HEADER
38 #define TEMPLATE_LAPACK_STEBZ_HEADER
39 
40 
41 template<class Treal>
42 int template_lapack_stebz(const char *range, const char *order, const integer *n, const Treal
43  *vl, const Treal *vu, const integer *il, const integer *iu, const Treal *abstol,
44  const Treal *d__, const Treal *e, integer *m, integer *nsplit,
45  Treal *w, integer *iblock, integer *isplit, Treal *work,
46  integer *iwork, integer *info)
47 {
48 /* -- LAPACK routine (version 3.0) --
49  Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
50  Courant Institute, Argonne National Lab, and Rice University
51  June 30, 1999
52 
53 
54  Purpose
55  =======
56 
57  DSTEBZ computes the eigenvalues of a symmetric tridiagonal
58  matrix T. The user may ask for all eigenvalues, all eigenvalues
59  in the half-open interval (VL, VU], or the IL-th through IU-th
60  eigenvalues.
61 
62  To avoid overflow, the matrix must be scaled so that its
63  largest element is no greater than overflow**(1/2) *
64  underflow**(1/4) in absolute value, and for greatest
65  accuracy, it should not be much smaller than that.
66 
67  See W. Kahan "Accurate Eigenvalues of a Symmetric Tridiagonal
68  Matrix", Report CS41, Computer Science Dept., Stanford
69  University, July 21, 1966.
70 
71  Arguments
72  =========
73 
74  RANGE (input) CHARACTER
75  = 'A': ("All") all eigenvalues will be found.
76  = 'V': ("Value") all eigenvalues in the half-open interval
77  (VL, VU] will be found.
78  = 'I': ("Index") the IL-th through IU-th eigenvalues (of the
79  entire matrix) will be found.
80 
81  ORDER (input) CHARACTER
82  = 'B': ("By Block") the eigenvalues will be grouped by
83  split-off block (see IBLOCK, ISPLIT) and
84  ordered from smallest to largest within
85  the block.
86  = 'E': ("Entire matrix")
87  the eigenvalues for the entire matrix
88  will be ordered from smallest to
89  largest.
90 
91  N (input) INTEGER
92  The order of the tridiagonal matrix T. N >= 0.
93 
94  VL (input) DOUBLE PRECISION
95  VU (input) DOUBLE PRECISION
96  If RANGE='V', the lower and upper bounds of the interval to
97  be searched for eigenvalues. Eigenvalues less than or equal
98  to VL, or greater than VU, will not be returned. VL < VU.
99  Not referenced if RANGE = 'A' or 'I'.
100 
101  IL (input) INTEGER
102  IU (input) INTEGER
103  If RANGE='I', the indices (in ascending order) of the
104  smallest and largest eigenvalues to be returned.
105  1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
106  Not referenced if RANGE = 'A' or 'V'.
107 
108  ABSTOL (input) DOUBLE PRECISION
109  The absolute tolerance for the eigenvalues. An eigenvalue
110  (or cluster) is considered to be located if it has been
111  determined to lie in an interval whose width is ABSTOL or
112  less. If ABSTOL is less than or equal to zero, then ULP*|T|
113  will be used, where |T| means the 1-norm of T.
114 
115  Eigenvalues will be computed most accurately when ABSTOL is
116  set to twice the underflow threshold 2*DLAMCH('S'), not zero.
117 
118  D (input) DOUBLE PRECISION array, dimension (N)
119  The n diagonal elements of the tridiagonal matrix T.
120 
121  E (input) DOUBLE PRECISION array, dimension (N-1)
122  The (n-1) off-diagonal elements of the tridiagonal matrix T.
123 
124  M (output) INTEGER
125  The actual number of eigenvalues found. 0 <= M <= N.
126  (See also the description of INFO=2,3.)
127 
128  NSPLIT (output) INTEGER
129  The number of diagonal blocks in the matrix T.
130  1 <= NSPLIT <= N.
131 
132  W (output) DOUBLE PRECISION array, dimension (N)
133  On exit, the first M elements of W will contain the
134  eigenvalues. (DSTEBZ may use the remaining N-M elements as
135  workspace.)
136 
137  IBLOCK (output) INTEGER array, dimension (N)
138  At each row/column j where E(j) is zero or small, the
139  matrix T is considered to split into a block diagonal
140  matrix. On exit, if INFO = 0, IBLOCK(i) specifies to which
141  block (from 1 to the number of blocks) the eigenvalue W(i)
142  belongs. (DSTEBZ may use the remaining N-M elements as
143  workspace.)
144 
145  ISPLIT (output) INTEGER array, dimension (N)
146  The splitting points, at which T breaks up into submatrices.
147  The first submatrix consists of rows/columns 1 to ISPLIT(1),
148  the second of rows/columns ISPLIT(1)+1 through ISPLIT(2),
149  etc., and the NSPLIT-th consists of rows/columns
150  ISPLIT(NSPLIT-1)+1 through ISPLIT(NSPLIT)=N.
151  (Only the first NSPLIT elements will actually be used, but
152  since the user cannot know a priori what value NSPLIT will
153  have, N words must be reserved for ISPLIT.)
154 
155  WORK (workspace) DOUBLE PRECISION array, dimension (4*N)
156 
157  IWORK (workspace) INTEGER array, dimension (3*N)
158 
159  INFO (output) INTEGER
160  = 0: successful exit
161  < 0: if INFO = -i, the i-th argument had an illegal value
162  > 0: some or all of the eigenvalues failed to converge or
163  were not computed:
164  =1 or 3: Bisection failed to converge for some
165  eigenvalues; these eigenvalues are flagged by a
166  negative block number. The effect is that the
167  eigenvalues may not be as accurate as the
168  absolute and relative tolerances. This is
169  generally caused by unexpectedly inaccurate
170  arithmetic.
171  =2 or 3: RANGE='I' only: Not all of the eigenvalues
172  IL:IU were found.
173  Effect: M < IU+1-IL
174  Cause: non-monotonic arithmetic, causing the
175  Sturm sequence to be non-monotonic.
176  Cure: recalculate, using RANGE='A', and pick
177  out eigenvalues IL:IU. In some cases,
178  increasing the PARAMETER "FUDGE" may
179  make things work.
180  = 4: RANGE='I', and the Gershgorin interval
181  initially used was too small. No eigenvalues
182  were computed.
183  Probable cause: your machine has sloppy
184  floating-point arithmetic.
185  Cure: Increase the PARAMETER "FUDGE",
186  recompile, and try again.
187 
188  Internal Parameters
189  ===================
190 
191  RELFAC DOUBLE PRECISION, default = 2.0e0
192  The relative tolerance. An interval (a,b] lies within
193  "relative tolerance" if b-a < RELFAC*ulp*max(|a|,|b|),
194  where "ulp" is the machine precision (distance from 1 to
195  the next larger floating point number.)
196 
197  FUDGE DOUBLE PRECISION, default = 2
198  A "fudge factor" to widen the Gershgorin intervals. Ideally,
199  a value of 1 should work, but on machines with sloppy
200  arithmetic, this needs to be larger. The default for
201  publicly released versions should be large enough to handle
202  the worst machine around. Note that this has no effect
203  on accuracy of the solution.
204 
205  =====================================================================
206 
207 
208  Parameter adjustments */
209  /* Table of constant values */
210  integer c__1 = 1;
211  integer c_n1 = -1;
212  integer c__3 = 3;
213  integer c__2 = 2;
214  integer c__0 = 0;
215 
216  /* System generated locals */
217  integer i__1, i__2, i__3;
218  Treal d__1, d__2, d__3, d__4, d__5;
219  /* Local variables */
220  integer iend, ioff, iout, itmp1, j, jdisc;
221  integer iinfo;
222  Treal atoli;
223  integer iwoff;
224  Treal bnorm;
225  integer itmax;
226  Treal wkill, rtoli, tnorm;
227  integer ib, jb, ie, je, nb;
228  Treal gl;
229  integer im, in;
230  integer ibegin;
231  Treal gu;
232  integer iw;
233  Treal wl;
234  integer irange, idiscl;
235  Treal safemn, wu;
236  integer idumma[1];
237  integer idiscu, iorder;
238  logical ncnvrg;
239  Treal pivmin;
240  logical toofew;
241  integer nwl;
242  Treal ulp, wlu, wul;
243  integer nwu;
244  Treal tmp1, tmp2;
245 
246 
247  --iwork;
248  --work;
249  --isplit;
250  --iblock;
251  --w;
252  --e;
253  --d__;
254 
255  /* Initialization added by Elias to get rid of compiler warnings. */
256  wlu = wul = 0;
257  /* Function Body */
258  *info = 0;
259 
260 /* Decode RANGE */
261 
262  if (template_blas_lsame(range, "A")) {
263  irange = 1;
264  } else if (template_blas_lsame(range, "V")) {
265  irange = 2;
266  } else if (template_blas_lsame(range, "I")) {
267  irange = 3;
268  } else {
269  irange = 0;
270  }
271 
272 /* Decode ORDER */
273 
274  if (template_blas_lsame(order, "B")) {
275  iorder = 2;
276  } else if (template_blas_lsame(order, "E")) {
277  iorder = 1;
278  } else {
279  iorder = 0;
280  }
281 
282 /* Check for Errors */
283 
284  if (irange <= 0) {
285  *info = -1;
286  } else if (iorder <= 0) {
287  *info = -2;
288  } else if (*n < 0) {
289  *info = -3;
290  } else if (irange == 2) {
291  if (*vl >= *vu) {
292  *info = -5;
293  }
294  } else if (irange == 3 && (*il < 1 || *il > maxMACRO(1,*n))) {
295  *info = -6;
296  } else if (irange == 3 && (*iu < minMACRO(*n,*il) || *iu > *n)) {
297  *info = -7;
298  }
299 
300  if (*info != 0) {
301  i__1 = -(*info);
302  template_blas_erbla("STEBZ ", &i__1);
303  return 0;
304  }
305 
306 /* Initialize error flags */
307 
308  *info = 0;
309  ncnvrg = FALSE_;
310  toofew = FALSE_;
311 
312 /* Quick return if possible */
313 
314  *m = 0;
315  if (*n == 0) {
316  return 0;
317  }
318 
319 /* Simplifications: */
320 
321  if (irange == 3 && *il == 1 && *iu == *n) {
322  irange = 1;
323  }
324 
325 /* Get machine constants
326  NB is the minimum vector length for vector bisection, or 0
327  if only scalar is to be done. */
328 
329  safemn = template_lapack_lamch("S", (Treal)0);
330  ulp = template_lapack_lamch("P", (Treal)0);
331  rtoli = ulp * 2.;
332  nb = template_lapack_ilaenv(&c__1, "DSTEBZ", " ", n, &c_n1, &c_n1, &c_n1, (ftnlen)6, (
333  ftnlen)1);
334  if (nb <= 1) {
335  nb = 0;
336  }
337 
338 /* Special Case when N=1 */
339 
340  if (*n == 1) {
341  *nsplit = 1;
342  isplit[1] = 1;
343  if (irange == 2 && (*vl >= d__[1] || *vu < d__[1])) {
344  *m = 0;
345  } else {
346  w[1] = d__[1];
347  iblock[1] = 1;
348  *m = 1;
349  }
350  return 0;
351  }
352 
353 /* Compute Splitting Points */
354 
355  *nsplit = 1;
356  work[*n] = 0.;
357  pivmin = 1.;
358 
359 /* DIR$ NOVECTOR */
360  i__1 = *n;
361  for (j = 2; j <= i__1; ++j) {
362 /* Computing 2nd power */
363  d__1 = e[j - 1];
364  tmp1 = d__1 * d__1;
365 /* Computing 2nd power */
366  d__2 = ulp;
367  if ((d__1 = d__[j] * d__[j - 1], absMACRO(d__1)) * (d__2 * d__2) + safemn
368  > tmp1) {
369  isplit[*nsplit] = j - 1;
370  ++(*nsplit);
371  work[j - 1] = 0.;
372  } else {
373  work[j - 1] = tmp1;
374  pivmin = maxMACRO(pivmin,tmp1);
375  }
376 /* L10: */
377  }
378  isplit[*nsplit] = *n;
379  pivmin *= safemn;
380 
381 /* Compute Interval and ATOLI */
382 
383  if (irange == 3) {
384 
385 /* RANGE='I': Compute the interval containing eigenvalues
386  IL through IU.
387 
388  Compute Gershgorin interval for entire (split) matrix
389  and use it as the initial interval */
390 
391  gu = d__[1];
392  gl = d__[1];
393  tmp1 = 0.;
394 
395  i__1 = *n - 1;
396  for (j = 1; j <= i__1; ++j) {
397  tmp2 = template_blas_sqrt(work[j]);
398 /* Computing MAX */
399  d__1 = gu, d__2 = d__[j] + tmp1 + tmp2;
400  gu = maxMACRO(d__1,d__2);
401 /* Computing MIN */
402  d__1 = gl, d__2 = d__[j] - tmp1 - tmp2;
403  gl = minMACRO(d__1,d__2);
404  tmp1 = tmp2;
405 /* L20: */
406  }
407 
408 /* Computing MAX */
409  d__1 = gu, d__2 = d__[*n] + tmp1;
410  gu = maxMACRO(d__1,d__2);
411 /* Computing MIN */
412  d__1 = gl, d__2 = d__[*n] - tmp1;
413  gl = minMACRO(d__1,d__2);
414 /* Computing MAX */
415  d__1 = absMACRO(gl), d__2 = absMACRO(gu);
416  tnorm = maxMACRO(d__1,d__2);
417  gl = gl - tnorm * 2. * ulp * *n - pivmin * 4.;
418  gu = gu + tnorm * 2. * ulp * *n + pivmin * 2.;
419 
420 /* Compute Iteration parameters */
421 
422  itmax = (integer) ((template_blas_log(tnorm + pivmin) - template_blas_log(pivmin)) / template_blas_log(2.)) + 2;
423  if (*abstol <= 0.) {
424  atoli = ulp * tnorm;
425  } else {
426  atoli = *abstol;
427  }
428 
429  work[*n + 1] = gl;
430  work[*n + 2] = gl;
431  work[*n + 3] = gu;
432  work[*n + 4] = gu;
433  work[*n + 5] = gl;
434  work[*n + 6] = gu;
435  iwork[1] = -1;
436  iwork[2] = -1;
437  iwork[3] = *n + 1;
438  iwork[4] = *n + 1;
439  iwork[5] = *il - 1;
440  iwork[6] = *iu;
441 
442  template_lapack_laebz(&c__3, &itmax, n, &c__2, &c__2, &nb, &atoli, &rtoli, &pivmin,
443  &d__[1], &e[1], &work[1], &iwork[5], &work[*n + 1], &work[*n
444  + 5], &iout, &iwork[1], &w[1], &iblock[1], &iinfo);
445 
446  if (iwork[6] == *iu) {
447  wl = work[*n + 1];
448  wlu = work[*n + 3];
449  nwl = iwork[1];
450  wu = work[*n + 4];
451  wul = work[*n + 2];
452  nwu = iwork[4];
453  } else {
454  wl = work[*n + 2];
455  wlu = work[*n + 4];
456  nwl = iwork[2];
457  wu = work[*n + 3];
458  wul = work[*n + 1];
459  nwu = iwork[3];
460  }
461 
462  if (nwl < 0 || nwl >= *n || nwu < 1 || nwu > *n) {
463  *info = 4;
464  return 0;
465  }
466  } else {
467 
468 /* RANGE='A' or 'V' -- Set ATOLI
469 
470  Computing MAX */
471  d__3 = absMACRO(d__[1]) + absMACRO(e[1]), d__4 = (d__1 = d__[*n], absMACRO(d__1)) + (
472  d__2 = e[*n - 1], absMACRO(d__2));
473  tnorm = maxMACRO(d__3,d__4);
474 
475  i__1 = *n - 1;
476  for (j = 2; j <= i__1; ++j) {
477 /* Computing MAX */
478  d__4 = tnorm, d__5 = (d__1 = d__[j], absMACRO(d__1)) + (d__2 = e[j - 1]
479  , absMACRO(d__2)) + (d__3 = e[j], absMACRO(d__3));
480  tnorm = maxMACRO(d__4,d__5);
481 /* L30: */
482  }
483 
484  if (*abstol <= 0.) {
485  atoli = ulp * tnorm;
486  } else {
487  atoli = *abstol;
488  }
489 
490  if (irange == 2) {
491  wl = *vl;
492  wu = *vu;
493  } else {
494  wl = 0.;
495  wu = 0.;
496  }
497  }
498 
499 /* Find Eigenvalues -- Loop Over Blocks and recompute NWL and NWU.
500  NWL accumulates the number of eigenvalues .le. WL,
501  NWU accumulates the number of eigenvalues .le. WU */
502 
503  *m = 0;
504  iend = 0;
505  *info = 0;
506  nwl = 0;
507  nwu = 0;
508 
509  i__1 = *nsplit;
510  for (jb = 1; jb <= i__1; ++jb) {
511  ioff = iend;
512  ibegin = ioff + 1;
513  iend = isplit[jb];
514  in = iend - ioff;
515 
516  if (in == 1) {
517 
518 /* Special Case -- IN=1 */
519 
520  if (irange == 1 || wl >= d__[ibegin] - pivmin) {
521  ++nwl;
522  }
523  if (irange == 1 || wu >= d__[ibegin] - pivmin) {
524  ++nwu;
525  }
526  if (irange == 1 || ( wl < d__[ibegin] - pivmin && wu >= d__[ibegin]
527  - pivmin ) ) {
528  ++(*m);
529  w[*m] = d__[ibegin];
530  iblock[*m] = jb;
531  }
532  } else {
533 
534 /* General Case -- IN > 1
535 
536  Compute Gershgorin Interval
537  and use it as the initial interval */
538 
539  gu = d__[ibegin];
540  gl = d__[ibegin];
541  tmp1 = 0.;
542 
543  i__2 = iend - 1;
544  for (j = ibegin; j <= i__2; ++j) {
545  tmp2 = (d__1 = e[j], absMACRO(d__1));
546 /* Computing MAX */
547  d__1 = gu, d__2 = d__[j] + tmp1 + tmp2;
548  gu = maxMACRO(d__1,d__2);
549 /* Computing MIN */
550  d__1 = gl, d__2 = d__[j] - tmp1 - tmp2;
551  gl = minMACRO(d__1,d__2);
552  tmp1 = tmp2;
553 /* L40: */
554  }
555 
556 /* Computing MAX */
557  d__1 = gu, d__2 = d__[iend] + tmp1;
558  gu = maxMACRO(d__1,d__2);
559 /* Computing MIN */
560  d__1 = gl, d__2 = d__[iend] - tmp1;
561  gl = minMACRO(d__1,d__2);
562 /* Computing MAX */
563  d__1 = absMACRO(gl), d__2 = absMACRO(gu);
564  bnorm = maxMACRO(d__1,d__2);
565  gl = gl - bnorm * 2. * ulp * in - pivmin * 2.;
566  gu = gu + bnorm * 2. * ulp * in + pivmin * 2.;
567 
568 /* Compute ATOLI for the current submatrix */
569 
570  if (*abstol <= 0.) {
571 /* Computing MAX */
572  d__1 = absMACRO(gl), d__2 = absMACRO(gu);
573  atoli = ulp * maxMACRO(d__1,d__2);
574  } else {
575  atoli = *abstol;
576  }
577 
578  if (irange > 1) {
579  if (gu < wl) {
580  nwl += in;
581  nwu += in;
582  goto L70;
583  }
584  gl = maxMACRO(gl,wl);
585  gu = minMACRO(gu,wu);
586  if (gl >= gu) {
587  goto L70;
588  }
589  }
590 
591 /* Set Up Initial Interval */
592 
593  work[*n + 1] = gl;
594  work[*n + in + 1] = gu;
595  template_lapack_laebz(&c__1, &c__0, &in, &in, &c__1, &nb, &atoli, &rtoli, &
596  pivmin, &d__[ibegin], &e[ibegin], &work[ibegin], idumma, &
597  work[*n + 1], &work[*n + (in << 1) + 1], &im, &iwork[1], &
598  w[*m + 1], &iblock[*m + 1], &iinfo);
599 
600  nwl += iwork[1];
601  nwu += iwork[in + 1];
602  iwoff = *m - iwork[1];
603 
604 /* Compute Eigenvalues */
605 
606  itmax = (integer) ((template_blas_log(gu - gl + pivmin) - template_blas_log(pivmin)) / template_blas_log(2.)
607  ) + 2;
608  template_lapack_laebz(&c__2, &itmax, &in, &in, &c__1, &nb, &atoli, &rtoli, &
609  pivmin, &d__[ibegin], &e[ibegin], &work[ibegin], idumma, &
610  work[*n + 1], &work[*n + (in << 1) + 1], &iout, &iwork[1],
611  &w[*m + 1], &iblock[*m + 1], &iinfo);
612 
613 /* Copy Eigenvalues Into W and IBLOCK
614  Use -JB for block number for unconverged eigenvalues. */
615 
616  i__2 = iout;
617  for (j = 1; j <= i__2; ++j) {
618  tmp1 = (work[j + *n] + work[j + in + *n]) * .5;
619 
620 /* Flag non-convergence. */
621 
622  if (j > iout - iinfo) {
623  ncnvrg = TRUE_;
624  ib = -jb;
625  } else {
626  ib = jb;
627  }
628  i__3 = iwork[j + in] + iwoff;
629  for (je = iwork[j] + 1 + iwoff; je <= i__3; ++je) {
630  w[je] = tmp1;
631  iblock[je] = ib;
632 /* L50: */
633  }
634 /* L60: */
635  }
636 
637  *m += im;
638  }
639 L70:
640  ;
641  }
642 
643 /* If RANGE='I', then (WL,WU) contains eigenvalues NWL+1,...,NWU
644  If NWL+1 < IL or NWU > IU, discard extra eigenvalues. */
645 
646  if (irange == 3) {
647  im = 0;
648  idiscl = *il - 1 - nwl;
649  idiscu = nwu - *iu;
650 
651  if (idiscl > 0 || idiscu > 0) {
652  i__1 = *m;
653  for (je = 1; je <= i__1; ++je) {
654  if (w[je] <= wlu && idiscl > 0) {
655  --idiscl;
656  } else if (w[je] >= wul && idiscu > 0) {
657  --idiscu;
658  } else {
659  ++im;
660  w[im] = w[je];
661  iblock[im] = iblock[je];
662  }
663 /* L80: */
664  }
665  *m = im;
666  }
667  if (idiscl > 0 || idiscu > 0) {
668 
669 /* Code to deal with effects of bad arithmetic:
670  Some low eigenvalues to be discarded are not in (WL,WLU],
671  or high eigenvalues to be discarded are not in (WUL,WU]
672  so just kill off the smallest IDISCL/largest IDISCU
673  eigenvalues, by simply finding the smallest/largest
674  eigenvalue(s).
675 
676  (If N(w) is monotone non-decreasing, this should never
677  happen.) */
678 
679  if (idiscl > 0) {
680  wkill = wu;
681  i__1 = idiscl;
682  for (jdisc = 1; jdisc <= i__1; ++jdisc) {
683  iw = 0;
684  i__2 = *m;
685  for (je = 1; je <= i__2; ++je) {
686  if (iblock[je] != 0 && (w[je] < wkill || iw == 0)) {
687  iw = je;
688  wkill = w[je];
689  }
690 /* L90: */
691  }
692  iblock[iw] = 0;
693 /* L100: */
694  }
695  }
696  if (idiscu > 0) {
697 
698  wkill = wl;
699  i__1 = idiscu;
700  for (jdisc = 1; jdisc <= i__1; ++jdisc) {
701  iw = 0;
702  i__2 = *m;
703  for (je = 1; je <= i__2; ++je) {
704  if (iblock[je] != 0 && (w[je] > wkill || iw == 0)) {
705  iw = je;
706  wkill = w[je];
707  }
708 /* L110: */
709  }
710  iblock[iw] = 0;
711 /* L120: */
712  }
713  }
714  im = 0;
715  i__1 = *m;
716  for (je = 1; je <= i__1; ++je) {
717  if (iblock[je] != 0) {
718  ++im;
719  w[im] = w[je];
720  iblock[im] = iblock[je];
721  }
722 /* L130: */
723  }
724  *m = im;
725  }
726  if (idiscl < 0 || idiscu < 0) {
727  toofew = TRUE_;
728  }
729  }
730 
731 /* If ORDER='B', do nothing -- the eigenvalues are already sorted
732  by block.
733  If ORDER='E', sort the eigenvalues from smallest to largest */
734 
735  if (iorder == 1 && *nsplit > 1) {
736  i__1 = *m - 1;
737  for (je = 1; je <= i__1; ++je) {
738  ie = 0;
739  tmp1 = w[je];
740  i__2 = *m;
741  for (j = je + 1; j <= i__2; ++j) {
742  if (w[j] < tmp1) {
743  ie = j;
744  tmp1 = w[j];
745  }
746 /* L140: */
747  }
748 
749  if (ie != 0) {
750  itmp1 = iblock[ie];
751  w[ie] = w[je];
752  iblock[ie] = iblock[je];
753  w[je] = tmp1;
754  iblock[je] = itmp1;
755  }
756 /* L150: */
757  }
758  }
759 
760  *info = 0;
761  if (ncnvrg) {
762  ++(*info);
763  }
764  if (toofew) {
765  *info += 2;
766  }
767  return 0;
768 
769 /* End of DSTEBZ */
770 
771 } /* dstebz_ */
772 
773 #endif
static const real gu
Definition: fun-pz81.c:68
#define absMACRO(x)
Definition: template_blas_common.h:47
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
int integer
Definition: template_blas_common.h:40
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
#define maxMACRO(a, b)
Definition: template_blas_common.h:45
#define minMACRO(a, b)
Definition: template_blas_common.h:46
int template_lapack_stebz(const char *range, const char *order, const integer *n, const Treal *vl, const Treal *vu, const integer *il, const integer *iu, const Treal *abstol, const Treal *d__, const Treal *e, integer *m, integer *nsplit, Treal *w, integer *iblock, integer *isplit, Treal *work, integer *iwork, integer *info)
Definition: template_lapack_stebz.h:42
Treal template_blas_log(Treal x)
int template_blas_erbla(const char *srname, integer *info)
Definition: template_blas_common.cc:146
Treal template_lapack_lamch(const char *cmach, Treal dummyReal)
Definition: template_lapack_lamch.h:202
bool logical
Definition: template_blas_common.h:41
#define TRUE_
Definition: template_lapack_common.h:42
#define FALSE_
Definition: template_lapack_common.h:43
int ftnlen
Definition: template_blas_common.h:42
Treal template_blas_sqrt(Treal x)
logical template_blas_lsame(const char *ca, const char *cb)
Definition: template_blas_common.cc:46