ergo
template_lapack_stemr.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_STEMR_HEADER
38#define TEMPLATE_LAPACK_STEMR_HEADER
39
40template<class Treal>
41int template_lapack_stemr(const char *jobz, const char *range, const integer *n, Treal *
42 d__, Treal *e, const Treal *vl, const Treal *vu, const integer *il,
43 const integer *iu, integer *m, Treal *w, Treal *z__, const integer *ldz,
44 const integer *nzc, integer *isuppz, logical *tryrac, Treal *work,
45 integer *lwork, integer *iwork, integer *liwork, integer *info)
46{
47 /* System generated locals */
48 integer z_dim1, z_offset, i__1, i__2;
49 Treal d__1, d__2;
50
51 /* Builtin functions */
52
53 /* Local variables */
54 integer i__, j;
55 Treal r1, r2;
56 integer jj;
57 Treal cs = 0;
58 integer in;
59 Treal sn = 0, wl, wu;
60 integer iil, iiu;
61 Treal eps, tmp;
62 integer indd, iend, jblk, wend;
63 Treal rmin, rmax;
64 integer itmp;
65 Treal tnrm;
66 integer inde2, itmp2;
67 Treal rtol1, rtol2;
68 Treal scale;
69 integer indgp;
70 integer iinfo, iindw, ilast;
71 integer lwmin;
72 logical wantz;
73 logical alleig;
74 integer ibegin;
75 logical indeig;
76 integer iindbl;
77 logical valeig;
78 integer wbegin;
79 Treal safmin;
80 Treal bignum;
81 integer inderr, iindwk, indgrs, offset;
82 Treal thresh;
83 integer iinspl, ifirst, indwrk, liwmin, nzcmin;
84 Treal pivmin;
85 integer nsplit;
86 Treal smlnum;
87 logical lquery, zquery;
88
89
90/* -- LAPACK computational routine (version 3.2) -- */
91/* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
92/* November 2006 */
93
94/* .. Scalar Arguments .. */
95/* .. */
96/* .. Array Arguments .. */
97/* .. */
98
99/* Purpose */
100/* ======= */
101
102/* DSTEMR computes selected eigenvalues and, optionally, eigenvectors */
103/* of a real symmetric tridiagonal matrix T. Any such unreduced matrix has */
104/* a well defined set of pairwise different real eigenvalues, the corresponding */
105/* real eigenvectors are pairwise orthogonal. */
106
107/* The spectrum may be computed either completely or partially by specifying */
108/* either an interval (VL,VU] or a range of indices IL:IU for the desired */
109/* eigenvalues. */
110
111/* Depending on the number of desired eigenvalues, these are computed either */
112/* by bisection or the dqds algorithm. Numerically orthogonal eigenvectors are */
113/* computed by the use of various suitable L D L^T factorizations near clusters */
114/* of close eigenvalues (referred to as RRRs, Relatively Robust */
115/* Representations). An informal sketch of the algorithm follows. */
116
117/* For each unreduced block (submatrix) of T, */
118/* (a) Compute T - sigma I = L D L^T, so that L and D */
119/* define all the wanted eigenvalues to high relative accuracy. */
120/* This means that small relative changes in the entries of D and L */
121/* cause only small relative changes in the eigenvalues and */
122/* eigenvectors. The standard (unfactored) representation of the */
123/* tridiagonal matrix T does not have this property in general. */
124/* (b) Compute the eigenvalues to suitable accuracy. */
125/* If the eigenvectors are desired, the algorithm attains full */
126/* accuracy of the computed eigenvalues only right before */
127/* the corresponding vectors have to be computed, see steps c) and d). */
128/* (c) For each cluster of close eigenvalues, select a new */
129/* shift close to the cluster, find a new factorization, and refine */
130/* the shifted eigenvalues to suitable accuracy. */
131/* (d) For each eigenvalue with a large enough relative separation compute */
132/* the corresponding eigenvector by forming a rank revealing twisted */
133/* factorization. Go back to (c) for any clusters that remain. */
134
135/* For more details, see: */
136/* - Inderjit S. Dhillon and Beresford N. Parlett: "Multiple representations */
137/* to compute orthogonal eigenvectors of symmetric tridiagonal matrices," */
138/* Linear Algebra and its Applications, 387(1), pp. 1-28, August 2004. */
139/* - Inderjit Dhillon and Beresford Parlett: "Orthogonal Eigenvectors and */
140/* Relative Gaps," SIAM Journal on Matrix Analysis and Applications, Vol. 25, */
141/* 2004. Also LAPACK Working Note 154. */
142/* - Inderjit Dhillon: "A new O(n^2) algorithm for the symmetric */
143/* tridiagonal eigenvalue/eigenvector problem", */
144/* Computer Science Division Technical Report No. UCB/CSD-97-971, */
145/* UC Berkeley, May 1997. */
146
147/* Notes: */
148/* 1.DSTEMR works only on machines which follow IEEE-754 */
149/* floating-point standard in their handling of infinities and NaNs. */
150/* This permits the use of efficient inner loops avoiding a check for */
151/* zero divisors. */
152
153/* Arguments */
154/* ========= */
155
156/* JOBZ (input) CHARACTER*1 */
157/* = 'N': Compute eigenvalues only; */
158/* = 'V': Compute eigenvalues and eigenvectors. */
159
160/* RANGE (input) CHARACTER*1 */
161/* = 'A': all eigenvalues will be found. */
162/* = 'V': all eigenvalues in the half-open interval (VL,VU] */
163/* will be found. */
164/* = 'I': the IL-th through IU-th eigenvalues will be found. */
165
166/* N (input) INTEGER */
167/* The order of the matrix. N >= 0. */
168
169/* D (input/output) DOUBLE PRECISION array, dimension (N) */
170/* On entry, the N diagonal elements of the tridiagonal matrix */
171/* T. On exit, D is overwritten. */
172
173/* E (input/output) DOUBLE PRECISION array, dimension (N) */
174/* On entry, the (N-1) subdiagonal elements of the tridiagonal */
175/* matrix T in elements 1 to N-1 of E. E(N) need not be set on */
176/* input, but is used internally as workspace. */
177/* On exit, E is overwritten. */
178
179/* VL (input) DOUBLE PRECISION */
180/* VU (input) DOUBLE PRECISION */
181/* If RANGE='V', the lower and upper bounds of the interval to */
182/* be searched for eigenvalues. VL < VU. */
183/* Not referenced if RANGE = 'A' or 'I'. */
184
185/* IL (input) INTEGER */
186/* IU (input) INTEGER */
187/* If RANGE='I', the indices (in ascending order) of the */
188/* smallest and largest eigenvalues to be returned. */
189/* 1 <= IL <= IU <= N, if N > 0. */
190/* Not referenced if RANGE = 'A' or 'V'. */
191
192/* M (output) INTEGER */
193/* The total number of eigenvalues found. 0 <= M <= N. */
194/* If RANGE = 'A', M = N, and if RANGE = 'I', M = IU-IL+1. */
195
196/* W (output) DOUBLE PRECISION array, dimension (N) */
197/* The first M elements contain the selected eigenvalues in */
198/* ascending order. */
199
200/* Z (output) DOUBLE PRECISION array, dimension (LDZ, max(1,M) ) */
201/* If JOBZ = 'V', and if INFO = 0, then the first M columns of Z */
202/* contain the orthonormal eigenvectors of the matrix T */
203/* corresponding to the selected eigenvalues, with the i-th */
204/* column of Z holding the eigenvector associated with W(i). */
205/* If JOBZ = 'N', then Z is not referenced. */
206/* Note: the user must ensure that at least max(1,M) columns are */
207/* supplied in the array Z; if RANGE = 'V', the exact value of M */
208/* is not known in advance and can be computed with a workspace */
209/* query by setting NZC = -1, see below. */
210
211/* LDZ (input) INTEGER */
212/* The leading dimension of the array Z. LDZ >= 1, and if */
213/* JOBZ = 'V', then LDZ >= max(1,N). */
214
215/* NZC (input) INTEGER */
216/* The number of eigenvectors to be held in the array Z. */
217/* If RANGE = 'A', then NZC >= max(1,N). */
218/* If RANGE = 'V', then NZC >= the number of eigenvalues in (VL,VU]. */
219/* If RANGE = 'I', then NZC >= IU-IL+1. */
220/* If NZC = -1, then a workspace query is assumed; the */
221/* routine calculates the number of columns of the array Z that */
222/* are needed to hold the eigenvectors. */
223/* This value is returned as the first entry of the Z array, and */
224/* no error message related to NZC is issued by XERBLA. */
225
226/* ISUPPZ (output) INTEGER ARRAY, dimension ( 2*max(1,M) ) */
227/* The support of the eigenvectors in Z, i.e., the indices */
228/* indicating the nonzero elements in Z. The i-th computed eigenvector */
229/* is nonzero only in elements ISUPPZ( 2*i-1 ) through */
230/* ISUPPZ( 2*i ). This is relevant in the case when the matrix */
231/* is split. ISUPPZ is only accessed when JOBZ is 'V' and N > 0. */
232
233/* TRYRAC (input/output) LOGICAL */
234/* If TRYRAC.EQ..TRUE., indicates that the code should check whether */
235/* the tridiagonal matrix defines its eigenvalues to high relative */
236/* accuracy. If so, the code uses relative-accuracy preserving */
237/* algorithms that might be (a bit) slower depending on the matrix. */
238/* If the matrix does not define its eigenvalues to high relative */
239/* accuracy, the code can uses possibly faster algorithms. */
240/* If TRYRAC.EQ..FALSE., the code is not required to guarantee */
241/* relatively accurate eigenvalues and can use the fastest possible */
242/* techniques. */
243/* On exit, a .TRUE. TRYRAC will be set to .FALSE. if the matrix */
244/* does not define its eigenvalues to high relative accuracy. */
245
246/* WORK (workspace/output) DOUBLE PRECISION array, dimension (LWORK) */
247/* On exit, if INFO = 0, WORK(1) returns the optimal */
248/* (and minimal) LWORK. */
249
250/* LWORK (input) INTEGER */
251/* The dimension of the array WORK. LWORK >= max(1,18*N) */
252/* if JOBZ = 'V', and LWORK >= max(1,12*N) if JOBZ = 'N'. */
253/* If LWORK = -1, then a workspace query is assumed; the routine */
254/* only calculates the optimal size of the WORK array, returns */
255/* this value as the first entry of the WORK array, and no error */
256/* message related to LWORK is issued by XERBLA. */
257
258/* IWORK (workspace/output) INTEGER array, dimension (LIWORK) */
259/* On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK. */
260
261/* LIWORK (input) INTEGER */
262/* The dimension of the array IWORK. LIWORK >= max(1,10*N) */
263/* if the eigenvectors are desired, and LIWORK >= max(1,8*N) */
264/* if only the eigenvalues are to be computed. */
265/* If LIWORK = -1, then a workspace query is assumed; the */
266/* routine only calculates the optimal size of the IWORK array, */
267/* returns this value as the first entry of the IWORK array, and */
268/* no error message related to LIWORK is issued by XERBLA. */
269
270/* INFO (output) INTEGER */
271/* On exit, INFO */
272/* = 0: successful exit */
273/* < 0: if INFO = -i, the i-th argument had an illegal value */
274/* > 0: if INFO = 1X, internal error in DLARRE, */
275/* if INFO = 2X, internal error in DLARRV. */
276/* Here, the digit X = ABS( IINFO ) < 10, where IINFO is */
277/* the nonzero error code returned by DLARRE or */
278/* DLARRV, respectively. */
279
280
281/* Further Details */
282/* =============== */
283
284/* Based on contributions by */
285/* Beresford Parlett, University of California, Berkeley, USA */
286/* Jim Demmel, University of California, Berkeley, USA */
287/* Inderjit Dhillon, University of Texas, Austin, USA */
288/* Osni Marques, LBNL/NERSC, USA */
289/* Christof Voemel, University of California, Berkeley, USA */
290
291/* ===================================================================== */
292
293/* .. Parameters .. */
294/* .. */
295/* .. Local Scalars .. */
296/* .. */
297/* .. */
298/* .. External Functions .. */
299/* .. */
300/* .. External Subroutines .. */
301/* .. */
302/* .. Intrinsic Functions .. */
303/* .. */
304/* .. Executable Statements .. */
305
306/* Test the input parameters. */
307
308 /* Parameter adjustments */
309 /* Table of constant values */
310 integer c__1 = 1;
311 Treal c_b18 = .001;
312
313 --d__;
314 --e;
315 --w;
316 z_dim1 = *ldz;
317 z_offset = 1 + z_dim1;
318 z__ -= z_offset;
319 --isuppz;
320 --work;
321 --iwork;
322
323 /* Function Body */
324 wantz = template_blas_lsame(jobz, "V");
325 alleig = template_blas_lsame(range, "A");
326 valeig = template_blas_lsame(range, "V");
327 indeig = template_blas_lsame(range, "I");
328
329 lquery = *lwork == -1 || *liwork == -1;
330 zquery = *nzc == -1;
331/* DSTEMR needs WORK of size 6*N, IWORK of size 3*N. */
332/* In addition, DLARRE needs WORK of size 6*N, IWORK of size 5*N. */
333/* Furthermore, DLARRV needs WORK of size 12*N, IWORK of size 7*N. */
334 if (wantz) {
335 lwmin = *n * 18;
336 liwmin = *n * 10;
337 } else {
338/* need less workspace if only the eigenvalues are wanted */
339 lwmin = *n * 12;
340 liwmin = *n << 3;
341 }
342 wl = 0.;
343 wu = 0.;
344 iil = 0;
345 iiu = 0;
346 if (valeig) {
347/* We do not reference VL, VU in the cases RANGE = 'I','A' */
348/* The interval (WL, WU] contains all the wanted eigenvalues. */
349/* It is either given by the user or computed in DLARRE. */
350 wl = *vl;
351 wu = *vu;
352 } else if (indeig) {
353/* We do not reference IL, IU in the cases RANGE = 'V','A' */
354 iil = *il;
355 iiu = *iu;
356 }
357
358 *info = 0;
359 if (! (wantz || template_blas_lsame(jobz, "N"))) {
360 *info = -1;
361 } else if (! (alleig || valeig || indeig)) {
362 *info = -2;
363 } else if (*n < 0) {
364 *info = -3;
365 } else if (valeig && *n > 0 && wu <= wl) {
366 *info = -7;
367 } else if (indeig && (iil < 1 || iil > *n)) {
368 *info = -8;
369 } else if (indeig && (iiu < iil || iiu > *n)) {
370 *info = -9;
371 } else if (*ldz < 1 || ( wantz && *ldz < *n ) ) {
372 *info = -13;
373 } else if (*lwork < lwmin && ! lquery) {
374 *info = -17;
375 } else if (*liwork < liwmin && ! lquery) {
376 *info = -19;
377 }
378
379/* Get machine constants. */
380
381 safmin = template_lapack_lamch("Safe minimum", (Treal)0);
382 eps = template_lapack_lamch("Precision", (Treal)0);
383 smlnum = safmin / eps;
384 bignum = 1. / smlnum;
385 rmin = template_blas_sqrt(smlnum);
386/* Computing MIN */
387 d__1 = template_blas_sqrt(bignum), d__2 = 1. / template_blas_sqrt(template_blas_sqrt(safmin));
388 rmax = minMACRO(d__1,d__2);
389
390 if (*info == 0) {
391 work[1] = (Treal) lwmin;
392 iwork[1] = liwmin;
393
394 if (wantz && alleig) {
395 nzcmin = *n;
396 } else if (wantz && valeig) {
397 template_lapack_larrc("T", n, vl, vu, &d__[1], &e[1], &safmin, &nzcmin, &itmp, &
398 itmp2, info);
399 } else if (wantz && indeig) {
400 nzcmin = iiu - iil + 1;
401 } else {
402/* WANTZ .EQ. FALSE. */
403 nzcmin = 0;
404 }
405 if (zquery && *info == 0) {
406 z__[z_dim1 + 1] = (Treal) nzcmin;
407 } else if (*nzc < nzcmin && ! zquery) {
408 *info = -14;
409 }
410 }
411 if (*info != 0) {
412
413 i__1 = -(*info);
414 template_blas_erbla("DSTEMR", &i__1);
415
416 return 0;
417 } else if (lquery || zquery) {
418 return 0;
419 }
420
421/* Handle N = 0, 1, and 2 cases immediately */
422
423 *m = 0;
424 if (*n == 0) {
425 return 0;
426 }
427
428 if (*n == 1) {
429 if (alleig || indeig) {
430 *m = 1;
431 w[1] = d__[1];
432 } else {
433 if (wl < d__[1] && wu >= d__[1]) {
434 *m = 1;
435 w[1] = d__[1];
436 }
437 }
438 if (wantz && ! zquery) {
439 z__[z_dim1 + 1] = 1.;
440 isuppz[1] = 1;
441 isuppz[2] = 1;
442 }
443 return 0;
444 }
445
446 if (*n == 2) {
447 if (! wantz) {
448 template_lapack_lae2(&d__[1], &e[1], &d__[2], &r1, &r2);
449 } else if (wantz && ! zquery) {
450 template_lapack_laev2(&d__[1], &e[1], &d__[2], &r1, &r2, &cs, &sn);
451 }
452 if (alleig || ( valeig && r2 > wl && r2 <= wu ) || ( indeig && iil == 1 ) ) {
453 ++(*m);
454 w[*m] = r2;
455 if (wantz && ! zquery) {
456 z__[*m * z_dim1 + 1] = -sn;
457 z__[*m * z_dim1 + 2] = cs;
458/* Note: At most one of SN and CS can be zero. */
459 if (sn != 0.) {
460 if (cs != 0.) {
461 isuppz[(*m << 1) - 1] = 1;
462 isuppz[(*m << 1) - 1] = 2;
463 } else {
464 isuppz[(*m << 1) - 1] = 1;
465 isuppz[(*m << 1) - 1] = 1;
466 }
467 } else {
468 isuppz[(*m << 1) - 1] = 2;
469 isuppz[*m * 2] = 2;
470 }
471 }
472 }
473 if (alleig || ( valeig && r1 > wl && r1 <= wu ) || ( indeig && iiu == 2 ) ) {
474 ++(*m);
475 w[*m] = r1;
476 if (wantz && ! zquery) {
477 z__[*m * z_dim1 + 1] = cs;
478 z__[*m * z_dim1 + 2] = sn;
479/* Note: At most one of SN and CS can be zero. */
480 if (sn != 0.) {
481 if (cs != 0.) {
482 isuppz[(*m << 1) - 1] = 1;
483 isuppz[(*m << 1) - 1] = 2;
484 } else {
485 isuppz[(*m << 1) - 1] = 1;
486 isuppz[(*m << 1) - 1] = 1;
487 }
488 } else {
489 isuppz[(*m << 1) - 1] = 2;
490 isuppz[*m * 2] = 2;
491 }
492 }
493 }
494 return 0;
495 }
496/* Continue with general N */
497 indgrs = 1;
498 inderr = (*n << 1) + 1;
499 indgp = *n * 3 + 1;
500 indd = (*n << 2) + 1;
501 inde2 = *n * 5 + 1;
502 indwrk = *n * 6 + 1;
503
504 iinspl = 1;
505 iindbl = *n + 1;
506 iindw = (*n << 1) + 1;
507 iindwk = *n * 3 + 1;
508
509/* Scale matrix to allowable range, if necessary. */
510/* The allowable range is related to the PIVMIN parameter; see the */
511/* comments in DLARRD. The preference for scaling small values */
512/* up is heuristic; we expect users' matrices not to be close to the */
513/* RMAX threshold. */
514
515 scale = 1.;
516 tnrm = template_lapack_lanst("M", n, &d__[1], &e[1]);
517 if (tnrm > 0. && tnrm < rmin) {
518 scale = rmin / tnrm;
519 } else if (tnrm > rmax) {
520 scale = rmax / tnrm;
521 }
522 if (scale != 1.) {
523 template_blas_scal(n, &scale, &d__[1], &c__1);
524 i__1 = *n - 1;
525 template_blas_scal(&i__1, &scale, &e[1], &c__1);
526 tnrm *= scale;
527 if (valeig) {
528/* If eigenvalues in interval have to be found, */
529/* scale (WL, WU] accordingly */
530 wl *= scale;
531 wu *= scale;
532 }
533 }
534
535/* Compute the desired eigenvalues of the tridiagonal after splitting */
536/* into smaller subblocks if the corresponding off-diagonal elements */
537/* are small */
538/* THRESH is the splitting parameter for DLARRE */
539/* A negative THRESH forces the old splitting criterion based on the */
540/* size of the off-diagonal. A positive THRESH switches to splitting */
541/* which preserves relative accuracy. */
542
543 if (*tryrac) {
544/* Test whether the matrix warrants the more expensive relative approach. */
545 template_lapack_larrr(n, &d__[1], &e[1], &iinfo);
546 } else {
547/* The user does not care about relative accurately eigenvalues */
548 iinfo = -1;
549 }
550/* Set the splitting criterion */
551 if (iinfo == 0) {
552 thresh = eps;
553 } else {
554 thresh = -eps;
555/* relative accuracy is desired but T does not guarantee it */
556 *tryrac = FALSE_;
557 }
558
559 if (*tryrac) {
560/* Copy original diagonal, needed to guarantee relative accuracy */
561 template_blas_copy(n, &d__[1], &c__1, &work[indd], &c__1);
562 }
563/* Store the squares of the offdiagonal values of T */
564 i__1 = *n - 1;
565 for (j = 1; j <= i__1; ++j) {
566/* Computing 2nd power */
567 d__1 = e[j];
568 work[inde2 + j - 1] = d__1 * d__1;
569/* L5: */
570 }
571/* Set the tolerance parameters for bisection */
572 if (! wantz) {
573/* DLARRE computes the eigenvalues to full precision. */
574 rtol1 = eps * 4.;
575 rtol2 = eps * 4.;
576 } else {
577/* DLARRE computes the eigenvalues to less than full precision. */
578/* DLARRV will refine the eigenvalue approximations, and we can */
579/* need less accurate initial bisection in DLARRE. */
580/* Note: these settings do only affect the subset case and DLARRE */
581 rtol1 = template_blas_sqrt(eps);
582/* Computing MAX */
583 d__1 = template_blas_sqrt(eps) * .005, d__2 = eps * 4.;
584 rtol2 = maxMACRO(d__1,d__2);
585 }
586 template_lapack_larre(range, n, &wl, &wu, &iil, &iiu, &d__[1], &e[1], &work[inde2], &
587 rtol1, &rtol2, &thresh, &nsplit, &iwork[iinspl], m, &w[1], &work[
588 inderr], &work[indgp], &iwork[iindbl], &iwork[iindw], &work[
589 indgrs], &pivmin, &work[indwrk], &iwork[iindwk], &iinfo);
590 if (iinfo != 0) {
591 *info = absMACRO(iinfo) + 10;
592 return 0;
593 }
594/* Note that if RANGE .NE. 'V', DLARRE computes bounds on the desired */
595/* part of the spectrum. All desired eigenvalues are contained in */
596/* (WL,WU] */
597 if (wantz) {
598
599/* Compute the desired eigenvectors corresponding to the computed */
600/* eigenvalues */
601
602 template_lapack_larrv(n, &wl, &wu, &d__[1], &e[1], &pivmin, &iwork[iinspl], m, &
603 c__1, m, &c_b18, &rtol1, &rtol2, &w[1], &work[inderr], &work[
604 indgp], &iwork[iindbl], &iwork[iindw], &work[indgrs], &z__[
605 z_offset], ldz, &isuppz[1], &work[indwrk], &iwork[iindwk], &
606 iinfo);
607 if (iinfo != 0) {
608 *info = absMACRO(iinfo) + 20;
609 return 0;
610 }
611 } else {
612/* DLARRE computes eigenvalues of the (shifted) root representation */
613/* DLARRV returns the eigenvalues of the unshifted matrix. */
614/* However, if the eigenvectors are not desired by the user, we need */
615/* to apply the corresponding shifts from DLARRE to obtain the */
616/* eigenvalues of the original matrix. */
617 i__1 = *m;
618 for (j = 1; j <= i__1; ++j) {
619 itmp = iwork[iindbl + j - 1];
620 w[j] += e[iwork[iinspl + itmp - 1]];
621/* L20: */
622 }
623 }
624
625 if (*tryrac) {
626/* Refine computed eigenvalues so that they are relatively accurate */
627/* with respect to the original matrix T. */
628 ibegin = 1;
629 wbegin = 1;
630 i__1 = iwork[iindbl + *m - 1];
631 for (jblk = 1; jblk <= i__1; ++jblk) {
632 iend = iwork[iinspl + jblk - 1];
633 in = iend - ibegin + 1;
634 wend = wbegin - 1;
635/* check if any eigenvalues have to be refined in this block */
636L36:
637 if (wend < *m) {
638 if (iwork[iindbl + wend] == jblk) {
639 ++wend;
640 goto L36;
641 }
642 }
643 if (wend < wbegin) {
644 ibegin = iend + 1;
645 goto L39;
646 }
647 offset = iwork[iindw + wbegin - 1] - 1;
648 ifirst = iwork[iindw + wbegin - 1];
649 ilast = iwork[iindw + wend - 1];
650 rtol2 = eps * 4.;
651 template_lapack_larrj(&in, &work[indd + ibegin - 1], &work[inde2 + ibegin - 1],
652 &ifirst, &ilast, &rtol2, &offset, &w[wbegin], &work[
653 inderr + wbegin - 1], &work[indwrk], &iwork[iindwk], &
654 pivmin, &tnrm, &iinfo);
655 ibegin = iend + 1;
656 wbegin = wend + 1;
657L39:
658 ;
659 }
660 }
661
662/* If matrix was scaled, then rescale eigenvalues appropriately. */
663
664 if (scale != 1.) {
665 d__1 = 1. / scale;
666 template_blas_scal(m, &d__1, &w[1], &c__1);
667 }
668
669/* If eigenvalues are not in increasing order, then sort them, */
670/* possibly along with eigenvectors. */
671
672 if (nsplit > 1) {
673 if (! wantz) {
674 template_lapack_lasrt("I", m, &w[1], &iinfo);
675 if (iinfo != 0) {
676 *info = 3;
677 return 0;
678 }
679 } else {
680 i__1 = *m - 1;
681 for (j = 1; j <= i__1; ++j) {
682 i__ = 0;
683 tmp = w[j];
684 i__2 = *m;
685 for (jj = j + 1; jj <= i__2; ++jj) {
686 if (w[jj] < tmp) {
687 i__ = jj;
688 tmp = w[jj];
689 }
690/* L50: */
691 }
692 if (i__ != 0) {
693 w[i__] = w[j];
694 w[j] = tmp;
695 if (wantz) {
696 template_blas_swap(n, &z__[i__ * z_dim1 + 1], &c__1, &z__[j *
697 z_dim1 + 1], &c__1);
698 itmp = isuppz[(i__ << 1) - 1];
699 isuppz[(i__ << 1) - 1] = isuppz[(j << 1) - 1];
700 isuppz[(j << 1) - 1] = itmp;
701 itmp = isuppz[i__ * 2];
702 isuppz[i__ * 2] = isuppz[j * 2];
703 isuppz[j * 2] = itmp;
704 }
705 }
706/* L60: */
707 }
708 }
709 }
710
711
712 work[1] = (Treal) lwmin;
713 iwork[1] = liwmin;
714 return 0;
715
716/* End of DSTEMR */
717
718} /* dstemr_ */
719
720
721#endif
Treal template_blas_sqrt(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
#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
int template_blas_copy(const integer *n, const Treal *dx, const integer *incx, Treal *dy, const integer *incy)
Definition template_blas_copy.h:42
int template_blas_scal(const integer *n, const Treal *da, Treal *dx, const integer *incx)
Definition template_blas_scal.h:43
int template_blas_swap(const integer *n, Treal *dx, const integer *incx, Treal *dy, const integer *incy)
Definition template_blas_swap.h:42
#define FALSE_
Definition template_lapack_common.h:43
int template_lapack_lae2(const Treal *a, const Treal *b, const Treal *c__, Treal *rt1, Treal *rt2)
Definition template_lapack_lae2.h:42
int template_lapack_laev2(Treal *a, Treal *b, Treal *c__, Treal *rt1, Treal *rt2, Treal *cs1, Treal *sn1)
Definition template_lapack_laev2.h:42
Treal template_lapack_lamch(const char *cmach, Treal dummyReal)
Definition template_lapack_lamch.h:202
Treal template_lapack_lanst(const char *norm, const integer *n, const Treal *d__, const Treal *e)
Definition template_lapack_lanst.h:42
int template_lapack_larrc(const char *jobt, const integer *n, const Treal *vl, const Treal *vu, Treal *d__, Treal *e, Treal *pivmin, integer *eigcnt, integer *lcnt, integer *rcnt, integer *info)
Definition template_lapack_larrc.h:41
int template_lapack_larre(const char *range, const integer *n, Treal *vl, Treal *vu, integer *il, integer *iu, Treal *d__, Treal *e, Treal *e2, Treal *rtol1, Treal *rtol2, Treal *spltol, integer *nsplit, integer *isplit, integer *m, Treal *w, Treal *werr, Treal *wgap, integer *iblock, integer *indexw, Treal *gers, Treal *pivmin, Treal *work, integer *iwork, integer *info)
Definition template_lapack_larre.h:46
int template_lapack_larrj(integer *n, Treal *d__, Treal *e2, integer *ifirst, integer *ilast, Treal *rtol, integer *offset, Treal *w, Treal *werr, Treal *work, integer *iwork, Treal *pivmin, Treal *spdiam, integer *info)
Definition template_lapack_larrj.h:41
int template_lapack_larrr(const integer *n, Treal *d__, Treal *e, integer *info)
Definition template_lapack_larrr.h:41
int template_lapack_larrv(const integer *n, Treal *vl, Treal *vu, Treal *d__, Treal *l, Treal *pivmin, integer *isplit, integer *m, integer *dol, integer *dou, Treal *minrgp, Treal *rtol1, Treal *rtol2, Treal *w, Treal *werr, Treal *wgap, integer *iblock, integer *indexw, Treal *gers, Treal *z__, const integer *ldz, integer *isuppz, Treal *work, integer *iwork, integer *info)
Definition template_lapack_larrv.h:45
int template_lapack_lasrt(const char *id, const integer *n, Treal *d__, integer *info)
Definition template_lapack_lasrt.h:42
int template_lapack_stemr(const char *jobz, const char *range, const integer *n, Treal *d__, Treal *e, const Treal *vl, const Treal *vu, const integer *il, const integer *iu, integer *m, Treal *w, Treal *z__, const integer *ldz, const integer *nzc, integer *isuppz, logical *tryrac, Treal *work, integer *lwork, integer *iwork, integer *liwork, integer *info)
Definition template_lapack_stemr.h:41