ergo
template_lapack_stebz.h
Go to the documentation of this file.
1/* Ergo, version 3.8.2, a program for linear scaling electronic structure
2 * calculations.
3 * Copyright (C) 2023 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
41template<class Treal>
42int 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 }
639L70:
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
Treal template_blas_sqrt(Treal x)
Treal template_blas_log(Treal x)
int template_blas_erbla(const char *srname, integer *info)
Definition: template_blas_common.cc:146
logical template_blas_lsame(const char *ca, const char *cb)
Definition: template_blas_common.cc:46
int integer
Definition: template_blas_common.h:40
#define absMACRO(x)
Definition: template_blas_common.h:47
int ftnlen
Definition: template_blas_common.h:42
#define minMACRO(a, b)
Definition: template_blas_common.h:46
#define maxMACRO(a, b)
Definition: template_blas_common.h:45
bool logical
Definition: template_blas_common.h:41
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 TRUE_
Definition: template_lapack_common.h:42
#define FALSE_
Definition: template_lapack_common.h:43
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
Treal template_lapack_lamch(const char *cmach, Treal dummyReal)
Definition: template_lapack_lamch.h:202
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