61 template<
typename Treal,
typename Tmatrix>
68 const Treal trunc = 0,
81 Treal
homo(Treal tol = 1e-15
86 Treal
lumo(Treal tol = 1e-15
96 void print_data(
int const start,
int const stop)
const;
145 template<
typename Treal,
typename Tmatrix>
148 Fun(
int const*
const p,
int const pl, Treal
const s)
149 :pol(p), pollength(pl), shift(s) {
150 assert(shift <= 1 && shift >= 0);
151 assert(pollength >= 0);
153 Treal
eval(Treal
const x)
const {
155 for (
int ind = 0; ind < pollength; ind++ )
156 y = 2 * pol[ind] * y + (1 - 2 * pol[ind]) * y * y;
171 template<
typename Treal,
typename Tmatrix>
174 const Treal trunc,
const int maxmm)
175 :X(F), D(DM), n(size), nocc(noc), frob_trunc(trunc), maxmul(maxmm),
176 lmin(0), lmax(0), nmul(0), nmul_firstpart(0),
177 idemerror(0), tracediff(0), polys(0) {
182 X.add_identity(-
lmax);
193 template<
typename Treal,
typename Tmatrix>
196 assert(nmul_firstpart == 0);
197 Treal delta, beta, trD2;
202 if (nmul >= maxmul) {
205 "Purification reached maxmul"
206 " without convergence", maxmul);
208 if (tracediff[nmul] > 0) {
213 D = -ONE * X * X + TWO * D;
216 D.frob_thresh(frob_trunc);
217 idemerror[nmul] = Tmatrix::frob_diff(D, X);
219 tracediff[nmul] = D.trace() - nocc;
223 if (idemerror[nmul - 1] < 1 / (Treal)4 &&
226 trD2 = (tracediff[nmul] + nocc -
227 2 * polys[nmul - 1] * (tracediff[nmul - 1] + nocc)) /
228 (1 - 2 * polys[nmul - 1]);
231 while (ind > 0 && polys[ind] == polys[ind - 1]) {
232 delta = delta + frob_trunc;
237 }
while((trD2 + beta * (1 + delta) * n - (1 + delta + beta) *
238 (tracediff[nmul - 1] + nocc)) /
239 ((1 + 2 * delta) * (delta + beta)) < n - nocc - 1 ||
240 (trD2 - delta * (1 - beta) * n - (1 - delta - beta) *
241 (tracediff[nmul - 1] + nocc)) /
242 ((1 + 2 * delta) * (delta + beta)) < nocc - 1);
250 if (tracediff[nmul - 1] > 0) {
253 D = -ONE * X * X + TWO * D;
260 D.frob_thresh(frob_trunc);
261 idemerror[nmul] = Tmatrix::frob_diff(D, X);
263 tracediff[nmul] = D.trace() - nocc;
265 nmul_firstpart = nmul;
269 if (nmul + 1 >= maxmul) {
272 "Purification reached maxmul"
273 " without convergence", maxmul);
275 if (tracediff[nmul] > 0) {
277 idemerror[nmul] = Tmatrix::frob_diff(D, X);
281 tracediff[nmul] = D.trace() - nocc;
283 D = -ONE * X * X + TWO * D;
284 idemerror[nmul] = Tmatrix::frob_diff(D, X);
287 tracediff[nmul] = D.trace() - nocc;
291 X = -ONE * D * D + TWO * X;
292 idemerror[nmul] = Tmatrix::frob_diff(D, X);
295 tracediff[nmul] = X.trace() - nocc;
298 idemerror[nmul] = Tmatrix::frob_diff(D, X);
301 tracediff[nmul] = D.trace() - nocc;
303 D.frob_thresh(frob_trunc);
305 }
while (idemerror[nmul - 1] > frob_trunc);
312 template<
typename Treal,
typename Tmatrix>
314 Fun const fermifun(polys, nmul, 0.5);
315 Treal chempot =
bisection(fermifun, (Treal)0, (Treal)1, tol);
316 return (lmin - lmax) * chempot + lmax;
319 template<
typename Treal,
typename Tmatrix>
323 for (
int mul = nmul_firstpart; mul < nmul; mul++) {
324 if (idemerror[mul] < 1.0 / 4) {
326 tmp =
bisection(homofun, (Treal)0, (Treal)1, tol);
331 homo = tmp > homo ? tmp : homo;
334 return (lmin - lmax) * homo + lmax;
337 template<
typename Treal,
typename Tmatrix>
341 for (
int mul = nmul_firstpart; mul < nmul; mul++) {
342 if (idemerror[mul] < 1.0 / 4) {
344 tmp =
bisection(lumofun, (Treal)0, (Treal)1, tol);
349 lumo = tmp < lumo ? tmp : lumo;
352 return (lmin - lmax) * lumo + lmax;
355 template<
typename Treal,
typename Tmatrix>
357 for (
int ind = start; ind < stop; ind ++) {
358 std::cout <<
"Iteration: " << ind
359 <<
" Idempotency error: " << idemerror[ind]
360 <<
" Tracediff: " << tracediff[ind]
361 <<
" Poly: " << polys[ind]