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