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