ergo
template_lapack_larrd.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_LARRD_HEADER
38#define TEMPLATE_LAPACK_LARRD_HEADER
39
40template<class Treal>
41int 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 }
644L70:
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
static const real gu
Definition fun-pz81.c:68
Treal template_blas_log(Treal x)
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_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