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