ergo
template_lapack_stemr.h
Go to the documentation of this file.
1 /* Ergo, version 3.8, a program for linear scaling electronic structure
2  * calculations.
3  * Copyright (C) 2019 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 
40 template<class Treal>
41 int 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 */
636 L36:
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;
657 L39:
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
template_blas_swap
int template_blas_swap(const integer *n, Treal *dx, const integer *incx, Treal *dy, const integer *incy)
Definition: template_blas_swap.h:42
template_blas_sqrt
Treal template_blas_sqrt(Treal x)
template_lapack_stemr
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
template_lapack_laev2
int template_lapack_laev2(Treal *a, Treal *b, Treal *c__, Treal *rt1, Treal *rt2, Treal *cs1, Treal *sn1)
Definition: template_lapack_laev2.h:42
template_lapack_lamch
Treal template_lapack_lamch(const char *cmach, Treal dummyReal)
Definition: template_lapack_lamch.h:202
template_blas_scal
int template_blas_scal(const integer *n, const Treal *da, Treal *dx, const integer *incx)
Definition: template_blas_scal.h:43
minMACRO
#define minMACRO(a, b)
Definition: template_blas_common.h:46
absMACRO
#define absMACRO(x)
Definition: template_blas_common.h:47
template_lapack_lanst
Treal template_lapack_lanst(const char *norm, const integer *n, const Treal *d__, const Treal *e)
Definition: template_lapack_lanst.h:42
template_lapack_larrc
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
logical
bool logical
Definition: template_blas_common.h:41
template_lapack_larrj
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
template_lapack_larrv
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
template_blas_erbla
int template_blas_erbla(const char *srname, integer *info)
Definition: template_blas_common.cc:146
template_lapack_larrr
int template_lapack_larrr(const integer *n, Treal *d__, Treal *e, integer *info)
Definition: template_lapack_larrr.h:41
template_blas_copy
int template_blas_copy(const integer *n, const Treal *dx, const integer *incx, Treal *dy, const integer *incy)
Definition: template_blas_copy.h:42
template_blas_lsame
logical template_blas_lsame(const char *ca, const char *cb)
Definition: template_blas_common.cc:46
integer
int integer
Definition: template_blas_common.h:40
template_lapack_lae2
int template_lapack_lae2(const Treal *a, const Treal *b, const Treal *c__, Treal *rt1, Treal *rt2)
Definition: template_lapack_lae2.h:42
FALSE_
#define FALSE_
Definition: template_lapack_common.h:43
maxMACRO
#define maxMACRO(a, b)
Definition: template_blas_common.h:45
template_lapack_larre
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
template_lapack_lasrt
int template_lapack_lasrt(const char *id, const integer *n, Treal *d__, integer *info)
Definition: template_lapack_lasrt.h:42