ergo
template_lapack_stevx.h
Go to the documentation of this file.
1 /* Ergo, version 3.4, a program for linear scaling electronic structure
2  * calculations.
3  * Copyright (C) 2014 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_STEVX_HEADER
36 #define TEMPLATE_LAPACK_STEVX_HEADER
37 
38 
39 template<class Treal>
40 int template_lapack_stevx(const char *jobz, const char *range, const integer *n, Treal *
41  d__, Treal *e, const Treal *vl, const Treal *vu, const integer *il,
42  const integer *iu, const Treal *abstol, integer *m, Treal *w,
43  Treal *z__, const integer *ldz, Treal *work, integer *iwork,
44  integer *ifail, integer *info)
45 {
46 /* -- LAPACK driver routine (version 3.0) --
47  Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
48  Courant Institute, Argonne National Lab, and Rice University
49  June 30, 1999
50 
51 
52  Purpose
53  =======
54 
55  DSTEVX computes selected eigenvalues and, optionally, eigenvectors
56  of a real symmetric tridiagonal matrix A. Eigenvalues and
57  eigenvectors can be selected by specifying either a range of values
58  or a range of indices for the desired eigenvalues.
59 
60  Arguments
61  =========
62 
63  JOBZ (input) CHARACTER*1
64  = 'N': Compute eigenvalues only;
65  = 'V': Compute eigenvalues and eigenvectors.
66 
67  RANGE (input) CHARACTER*1
68  = 'A': all eigenvalues will be found.
69  = 'V': all eigenvalues in the half-open interval (VL,VU]
70  will be found.
71  = 'I': the IL-th through IU-th eigenvalues will be found.
72 
73  N (input) INTEGER
74  The order of the matrix. N >= 0.
75 
76  D (input/output) DOUBLE PRECISION array, dimension (N)
77  On entry, the n diagonal elements of the tridiagonal matrix
78  A.
79  On exit, D may be multiplied by a constant factor chosen
80  to avoid over/underflow in computing the eigenvalues.
81 
82  E (input/output) DOUBLE PRECISION array, dimension (N)
83  On entry, the (n-1) subdiagonal elements of the tridiagonal
84  matrix A in elements 1 to N-1 of E; E(N) need not be set.
85  On exit, E may be multiplied by a constant factor chosen
86  to avoid over/underflow in computing the eigenvalues.
87 
88  VL (input) DOUBLE PRECISION
89  VU (input) DOUBLE PRECISION
90  If RANGE='V', the lower and upper bounds of the interval to
91  be searched for eigenvalues. VL < VU.
92  Not referenced if RANGE = 'A' or 'I'.
93 
94  IL (input) INTEGER
95  IU (input) INTEGER
96  If RANGE='I', the indices (in ascending order) of the
97  smallest and largest eigenvalues to be returned.
98  1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
99  Not referenced if RANGE = 'A' or 'V'.
100 
101  ABSTOL (input) DOUBLE PRECISION
102  The absolute error tolerance for the eigenvalues.
103  An approximate eigenvalue is accepted as converged
104  when it is determined to lie in an interval [a,b]
105  of width less than or equal to
106 
107  ABSTOL + EPS * max( |a|,|b| ) ,
108 
109  where EPS is the machine precision. If ABSTOL is less
110  than or equal to zero, then EPS*|T| will be used in
111  its place, where |T| is the 1-norm of the tridiagonal
112  matrix.
113 
114  Eigenvalues will be computed most accurately when ABSTOL is
115  set to twice the underflow threshold 2*DLAMCH('S'), not zero.
116  If this routine returns with INFO>0, indicating that some
117  eigenvectors did not converge, try setting ABSTOL to
118  2*DLAMCH('S').
119 
120  See "Computing Small Singular Values of Bidiagonal Matrices
121  with Guaranteed High Relative Accuracy," by Demmel and
122  Kahan, LAPACK Working Note #3.
123 
124  M (output) INTEGER
125  The total number of eigenvalues found. 0 <= M <= N.
126  If RANGE = 'A', M = N, and if RANGE = 'I', M = IU-IL+1.
127 
128  W (output) DOUBLE PRECISION array, dimension (N)
129  The first M elements contain the selected eigenvalues in
130  ascending order.
131 
132  Z (output) DOUBLE PRECISION array, dimension (LDZ, max(1,M) )
133  If JOBZ = 'V', then if INFO = 0, the first M columns of Z
134  contain the orthonormal eigenvectors of the matrix A
135  corresponding to the selected eigenvalues, with the i-th
136  column of Z holding the eigenvector associated with W(i).
137  If an eigenvector fails to converge (INFO > 0), then that
138  column of Z contains the latest approximation to the
139  eigenvector, and the index of the eigenvector is returned
140  in IFAIL. If JOBZ = 'N', then Z is not referenced.
141  Note: the user must ensure that at least max(1,M) columns are
142  supplied in the array Z; if RANGE = 'V', the exact value of M
143  is not known in advance and an upper bound must be used.
144 
145  LDZ (input) INTEGER
146  The leading dimension of the array Z. LDZ >= 1, and if
147  JOBZ = 'V', LDZ >= max(1,N).
148 
149  WORK (workspace) DOUBLE PRECISION array, dimension (5*N)
150 
151  IWORK (workspace) INTEGER array, dimension (5*N)
152 
153  IFAIL (output) INTEGER array, dimension (N)
154  If JOBZ = 'V', then if INFO = 0, the first M elements of
155  IFAIL are zero. If INFO > 0, then IFAIL contains the
156  indices of the eigenvectors that failed to converge.
157  If JOBZ = 'N', then IFAIL is not referenced.
158 
159  INFO (output) INTEGER
160  = 0: successful exit
161  < 0: if INFO = -i, the i-th argument had an illegal value
162  > 0: if INFO = i, then i eigenvectors failed to converge.
163  Their indices are stored in array IFAIL.
164 
165  =====================================================================
166 
167 
168  Test the input parameters.
169 
170  Parameter adjustments */
171  /* Table of constant values */
172  integer c__1 = 1;
173 
174  /* System generated locals */
175  integer z_dim1, z_offset, i__1, i__2;
176  Treal d__1, d__2;
177  /* Local variables */
178  integer imax;
179  Treal rmin, rmax, tnrm;
180  integer itmp1, i__, j;
181  Treal sigma;
182  char order[1];
183  logical wantz;
184  integer jj;
185  logical alleig, indeig;
186  integer iscale, indibl;
187  logical valeig;
188  Treal safmin;
189  Treal bignum;
190  integer indisp;
191  integer indiwo;
192  integer indwrk;
193  integer nsplit;
194  Treal smlnum, eps, vll, vuu, tmp1;
195 #define z___ref(a_1,a_2) z__[(a_2)*z_dim1 + a_1]
196 
197 
198  --d__;
199  --e;
200  --w;
201  z_dim1 = *ldz;
202  z_offset = 1 + z_dim1 * 1;
203  z__ -= z_offset;
204  --work;
205  --iwork;
206  --ifail;
207 
208  /* Function Body */
209  wantz = template_blas_lsame(jobz, "V");
210  alleig = template_blas_lsame(range, "A");
211  valeig = template_blas_lsame(range, "V");
212  indeig = template_blas_lsame(range, "I");
213 
214  *info = 0;
215  if (! (wantz || template_blas_lsame(jobz, "N"))) {
216  *info = -1;
217  } else if (! (alleig || valeig || indeig)) {
218  *info = -2;
219  } else if (*n < 0) {
220  *info = -3;
221  } else {
222  if (valeig) {
223  if (*n > 0 && *vu <= *vl) {
224  *info = -7;
225  }
226  } else if (indeig) {
227  if (*il < 1 || *il > maxMACRO(1,*n)) {
228  *info = -8;
229  } else if (*iu < minMACRO(*n,*il) || *iu > *n) {
230  *info = -9;
231  }
232  }
233  }
234  if (*info == 0) {
235  if (*ldz < 1 || (wantz && *ldz < *n) ) {
236  *info = -14;
237  }
238  }
239 
240  if (*info != 0) {
241  i__1 = -(*info);
242  template_blas_erbla("STEVX ", &i__1);
243  return 0;
244  }
245 
246 /* Quick return if possible */
247 
248  *m = 0;
249  if (*n == 0) {
250  return 0;
251  }
252 
253  if (*n == 1) {
254  if (alleig || indeig) {
255  *m = 1;
256  w[1] = d__[1];
257  } else {
258  if (*vl < d__[1] && *vu >= d__[1]) {
259  *m = 1;
260  w[1] = d__[1];
261  }
262  }
263  if (wantz) {
264  z___ref(1, 1) = 1.;
265  }
266  return 0;
267  }
268 
269 /* Get machine constants. */
270 
271  safmin = template_lapack_lamch("Safe minimum", (Treal)0);
272  eps = template_lapack_lamch("Precision", (Treal)0);
273  smlnum = safmin / eps;
274  bignum = 1. / smlnum;
275  rmin = template_blas_sqrt(smlnum);
276 /* Computing MIN */
277  d__1 = template_blas_sqrt(bignum), d__2 = 1. / template_blas_sqrt(template_blas_sqrt(safmin));
278  rmax = minMACRO(d__1,d__2);
279 
280 /* Scale matrix to allowable range, if necessary. */
281 
282  iscale = 0;
283  if (valeig) {
284  vll = *vl;
285  vuu = *vu;
286  } else {
287  vll = 0.;
288  vuu = 0.;
289  }
290  tnrm = template_lapack_lanst("M", n, &d__[1], &e[1]);
291  if (tnrm > 0. && tnrm < rmin) {
292  iscale = 1;
293  sigma = rmin / tnrm;
294  } else if (tnrm > rmax) {
295  iscale = 1;
296  sigma = rmax / tnrm;
297  }
298  if (iscale == 1) {
299  template_blas_scal(n, &sigma, &d__[1], &c__1);
300  i__1 = *n - 1;
301  template_blas_scal(&i__1, &sigma, &e[1], &c__1);
302  if (valeig) {
303  vll = *vl * sigma;
304  vuu = *vu * sigma;
305  }
306  }
307 
308 /* If all eigenvalues are desired and ABSTOL is less than zero, then
309  call DSTERF or SSTEQR. If this fails for some eigenvalue, then
310  try DSTEBZ. */
311 
312  if ((alleig || (indeig && *il == 1 && *iu == *n) ) && *abstol <= 0.) {
313  template_blas_copy(n, &d__[1], &c__1, &w[1], &c__1);
314  i__1 = *n - 1;
315  template_blas_copy(&i__1, &e[1], &c__1, &work[1], &c__1);
316  indwrk = *n + 1;
317  if (! wantz) {
318  template_lapack_sterf(n, &w[1], &work[1], info);
319  } else {
320  template_lapack_steqr("I", n, &w[1], &work[1], &z__[z_offset], ldz, &work[
321  indwrk], info);
322  if (*info == 0) {
323  i__1 = *n;
324  for (i__ = 1; i__ <= i__1; ++i__) {
325  ifail[i__] = 0;
326 /* L10: */
327  }
328  }
329  }
330  if (*info == 0) {
331  *m = *n;
332  goto L20;
333  }
334  *info = 0;
335  }
336 
337 /* Otherwise, call DSTEBZ and, if eigenvectors are desired, SSTEIN. */
338 
339  if (wantz) {
340  *(unsigned char *)order = 'B';
341  } else {
342  *(unsigned char *)order = 'E';
343  }
344  indwrk = 1;
345  indibl = 1;
346  indisp = indibl + *n;
347  indiwo = indisp + *n;
348  template_lapack_stebz(range, order, n, &vll, &vuu, il, iu, abstol, &d__[1], &e[1], m, &
349  nsplit, &w[1], &iwork[indibl], &iwork[indisp], &work[indwrk], &
350  iwork[indiwo], info);
351 
352  if (wantz) {
353  template_lapack_stein(n, &d__[1], &e[1], m, &w[1], &iwork[indibl], &iwork[indisp], &
354  z__[z_offset], ldz, &work[indwrk], &iwork[indiwo], &ifail[1],
355  info);
356  }
357 
358 /* If matrix was scaled, then rescale eigenvalues appropriately. */
359 
360 L20:
361  if (iscale == 1) {
362  if (*info == 0) {
363  imax = *m;
364  } else {
365  imax = *info - 1;
366  }
367  d__1 = 1. / sigma;
368  template_blas_scal(&imax, &d__1, &w[1], &c__1);
369  }
370 
371 /* If eigenvalues are not in order, then sort them, along with
372  eigenvectors. */
373 
374  if (wantz) {
375  i__1 = *m - 1;
376  for (j = 1; j <= i__1; ++j) {
377  i__ = 0;
378  tmp1 = w[j];
379  i__2 = *m;
380  for (jj = j + 1; jj <= i__2; ++jj) {
381  if (w[jj] < tmp1) {
382  i__ = jj;
383  tmp1 = w[jj];
384  }
385 /* L30: */
386  }
387 
388  if (i__ != 0) {
389  itmp1 = iwork[indibl + i__ - 1];
390  w[i__] = w[j];
391  iwork[indibl + i__ - 1] = iwork[indibl + j - 1];
392  w[j] = tmp1;
393  iwork[indibl + j - 1] = itmp1;
394  template_blas_swap(n, &z___ref(1, i__), &c__1, &z___ref(1, j), &c__1);
395  if (*info != 0) {
396  itmp1 = ifail[i__];
397  ifail[i__] = ifail[j];
398  ifail[j] = itmp1;
399  }
400  }
401 /* L40: */
402  }
403  }
404 
405  return 0;
406 
407 /* End of DSTEVX */
408 
409 } /* dstevx_ */
410 
411 #undef z___ref
412 
413 
414 #endif
int template_blas_scal(const integer *n, const Treal *da, Treal *dx, const integer *incx)
Definition: template_blas_scal.h:41
int integer
Definition: template_blas_common.h:38
#define maxMACRO(a, b)
Definition: template_blas_common.h:43
#define minMACRO(a, b)
Definition: template_blas_common.h:44
int template_lapack_stebz(const char *range, const char *order, const integer *n, const Treal *vl, const Treal *vu, const integer *il, const integer *iu, const Treal *abstol, const Treal *d__, const Treal *e, integer *m, integer *nsplit, Treal *w, integer *iblock, integer *isplit, Treal *work, integer *iwork, integer *info)
Definition: template_lapack_stebz.h:40
int template_lapack_steqr(const char *compz, const integer *n, Treal *d__, Treal *e, Treal *z__, const integer *ldz, Treal *work, integer *info)
Definition: template_lapack_steqr.h:40
int template_blas_erbla(const char *srname, integer *info)
Definition: template_blas_common.cc:144
int template_blas_swap(const integer *n, Treal *dx, const integer *incx, Treal *dy, const integer *incy)
Definition: template_blas_swap.h:40
#define z___ref(a_1, a_2)
int template_lapack_stein(const integer *n, const Treal *d__, const Treal *e, const integer *m, const Treal *w, const integer *iblock, const integer *isplit, Treal *z__, const integer *ldz, Treal *work, integer *iwork, integer *ifail, integer *info)
Definition: template_lapack_stein.h:40
Treal template_lapack_lamch(const char *cmach, Treal dummyReal)
Definition: template_lapack_lamch.h:199
bool logical
Definition: template_blas_common.h:39
int template_blas_copy(const integer *n, const Treal *dx, const integer *incx, Treal *dy, const integer *incy)
Definition: template_blas_copy.h:40
int template_lapack_sterf(const integer *n, Treal *d__, Treal *e, integer *info)
Definition: template_lapack_sterf.h:41
Treal template_lapack_lanst(const char *norm, const integer *n, const Treal *d__, const Treal *e)
Definition: template_lapack_lanst.h:40
Treal template_blas_sqrt(Treal x)
logical template_blas_lsame(const char *ca, const char *cb)
Definition: template_blas_common.cc:44
int template_lapack_stevx(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, const Treal *abstol, integer *m, Treal *w, Treal *z__, const integer *ldz, Treal *work, integer *iwork, integer *ifail, integer *info)
Definition: template_lapack_stevx.h:40