ergo
template_lapack_larrc.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_LARRC_HEADER
38 #define TEMPLATE_LAPACK_LARRC_HEADER
39 
40 template<class Treal>
41 int template_lapack_larrc(const char *jobt, const integer *n, const Treal *vl,
42  const Treal *vu, Treal *d__, Treal *e, Treal *pivmin,
43  integer *eigcnt, integer *lcnt, integer *rcnt, integer *info)
44 {
45  /* System generated locals */
46  integer i__1;
47  Treal d__1;
48 
49  /* Local variables */
50  integer i__;
51  Treal sl, su, tmp, tmp2;
52  logical matt;
53  Treal lpivot, rpivot;
54 
55 
56 /* -- LAPACK auxiliary routine (version 3.2) -- */
57 /* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
58 /* November 2006 */
59 
60 /* .. Scalar Arguments .. */
61 /* .. */
62 /* .. Array Arguments .. */
63 /* .. */
64 
65 /* Purpose */
66 /* ======= */
67 
68 /* Find the number of eigenvalues of the symmetric tridiagonal matrix T */
69 /* that are in the interval (VL,VU] if JOBT = 'T', and of L D L^T */
70 /* if JOBT = 'L'. */
71 
72 /* Arguments */
73 /* ========= */
74 
75 /* JOBT (input) CHARACTER*1 */
76 /* = 'T': Compute Sturm count for matrix T. */
77 /* = 'L': Compute Sturm count for matrix L D L^T. */
78 
79 /* N (input) INTEGER */
80 /* The order of the matrix. N > 0. */
81 
82 /* VL (input) DOUBLE PRECISION */
83 /* VU (input) DOUBLE PRECISION */
84 /* The lower and upper bounds for the eigenvalues. */
85 
86 /* D (input) DOUBLE PRECISION array, dimension (N) */
87 /* JOBT = 'T': The N diagonal elements of the tridiagonal matrix T. */
88 /* JOBT = 'L': The N diagonal elements of the diagonal matrix D. */
89 
90 /* E (input) DOUBLE PRECISION array, dimension (N) */
91 /* JOBT = 'T': The N-1 offdiagonal elements of the matrix T. */
92 /* JOBT = 'L': The N-1 offdiagonal elements of the matrix L. */
93 
94 /* PIVMIN (input) DOUBLE PRECISION */
95 /* The minimum pivot in the Sturm sequence for T. */
96 
97 /* EIGCNT (output) INTEGER */
98 /* The number of eigenvalues of the symmetric tridiagonal matrix T */
99 /* that are in the interval (VL,VU] */
100 
101 /* LCNT (output) INTEGER */
102 /* RCNT (output) INTEGER */
103 /* The left and right negcounts of the interval. */
104 
105 /* INFO (output) INTEGER */
106 
107 /* Further Details */
108 /* =============== */
109 
110 /* Based on contributions by */
111 /* Beresford Parlett, University of California, Berkeley, USA */
112 /* Jim Demmel, University of California, Berkeley, USA */
113 /* Inderjit Dhillon, University of Texas, Austin, USA */
114 /* Osni Marques, LBNL/NERSC, USA */
115 /* Christof Voemel, University of California, Berkeley, USA */
116 
117 /* ===================================================================== */
118 
119 /* .. Parameters .. */
120 /* .. */
121 /* .. Local Scalars .. */
122 /* .. */
123 /* .. External Functions .. */
124 /* .. */
125 /* .. Executable Statements .. */
126 
127  /* Parameter adjustments */
128  --e;
129  --d__;
130 
131  /* Function Body */
132  *info = 0;
133  *lcnt = 0;
134  *rcnt = 0;
135  *eigcnt = 0;
136  matt = template_blas_lsame(jobt, "T");
137  if (matt) {
138 /* Sturm sequence count on T */
139  lpivot = d__[1] - *vl;
140  rpivot = d__[1] - *vu;
141  if (lpivot <= 0.) {
142  ++(*lcnt);
143  }
144  if (rpivot <= 0.) {
145  ++(*rcnt);
146  }
147  i__1 = *n - 1;
148  for (i__ = 1; i__ <= i__1; ++i__) {
149 /* Computing 2nd power */
150  d__1 = e[i__];
151  tmp = d__1 * d__1;
152  lpivot = d__[i__ + 1] - *vl - tmp / lpivot;
153  rpivot = d__[i__ + 1] - *vu - tmp / rpivot;
154  if (lpivot <= 0.) {
155  ++(*lcnt);
156  }
157  if (rpivot <= 0.) {
158  ++(*rcnt);
159  }
160 /* L10: */
161  }
162  } else {
163 /* Sturm sequence count on L D L^T */
164  sl = -(*vl);
165  su = -(*vu);
166  i__1 = *n - 1;
167  for (i__ = 1; i__ <= i__1; ++i__) {
168  lpivot = d__[i__] + sl;
169  rpivot = d__[i__] + su;
170  if (lpivot <= 0.) {
171  ++(*lcnt);
172  }
173  if (rpivot <= 0.) {
174  ++(*rcnt);
175  }
176  tmp = e[i__] * d__[i__] * e[i__];
177 
178  tmp2 = tmp / lpivot;
179  if (tmp2 == 0.) {
180  sl = tmp - *vl;
181  } else {
182  sl = sl * tmp2 - *vl;
183  }
184 
185  tmp2 = tmp / rpivot;
186  if (tmp2 == 0.) {
187  su = tmp - *vu;
188  } else {
189  su = su * tmp2 - *vu;
190  }
191 /* L20: */
192  }
193  lpivot = d__[*n] + sl;
194  rpivot = d__[*n] + su;
195  if (lpivot <= 0.) {
196  ++(*lcnt);
197  }
198  if (rpivot <= 0.) {
199  ++(*rcnt);
200  }
201  }
202  *eigcnt = *rcnt - *lcnt;
203  return 0;
204 
205 /* end of DLARRC */
206 
207 } /* dlarrc_ */
208 
209 #endif
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_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