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