63 template<
typename Treal,
typename Tmatrix>
70 const Treal trunc = 0,
83 Treal
homo(Treal tol = 1e-15
88 Treal
lumo(Treal tol = 1e-15
98 void print_data(
int const start,
int const stop)
const;
147 template<
typename Treal,
typename Tmatrix>
150 Fun(
int const*
const p,
int const pl, Treal
const s)
151 :pol(p), pollength(pl), shift(s) {
152 assert(shift <= 1 && shift >= 0);
153 assert(pollength >= 0);
155 Treal
eval(Treal
const x)
const {
157 for (
int ind = 0; ind < pollength; ind++ )
158 y = 2 * pol[ind] * y + (1 - 2 * pol[ind]) * y * y;
173 template<
typename Treal,
typename Tmatrix>
176 const Treal trunc,
const int maxmm)
177 :X(F), D(DM), n(size), nocc(noc), frob_trunc(trunc), maxmul(maxmm),
178 lmin(0), lmax(0), nmul(0), nmul_firstpart(0),
179 idemerror(0), tracediff(0), polys(0) {
184 X.add_identity(-
lmax);
195 template<
typename Treal,
typename Tmatrix>
198 assert(nmul_firstpart == 0);
199 Treal delta, beta, trD2;
204 if (nmul >= maxmul) {
207 "Purification reached maxmul"
208 " without convergence", maxmul);
210 if (tracediff[nmul] > 0) {
215 D = -ONE * X * X + TWO * D;
218 D.frob_thresh(frob_trunc);
219 idemerror[nmul] = Tmatrix::frob_diff(D, X);
221 tracediff[nmul] = D.trace() - nocc;
225 if (idemerror[nmul - 1] < 1 / (Treal)4 &&
228 trD2 = (tracediff[nmul] + nocc -
229 2 * polys[nmul - 1] * (tracediff[nmul - 1] + nocc)) /
230 (1 - 2 * polys[nmul - 1]);
233 while (ind > 0 && polys[ind] == polys[ind - 1]) {
234 delta = delta + frob_trunc;
239 }
while((trD2 + beta * (1 + delta) * n - (1 + delta + beta) *
240 (tracediff[nmul - 1] + nocc)) /
241 ((1 + 2 * delta) * (delta + beta)) < n - nocc - 1 ||
242 (trD2 - delta * (1 - beta) * n - (1 - delta - beta) *
243 (tracediff[nmul - 1] + nocc)) /
244 ((1 + 2 * delta) * (delta + beta)) < nocc - 1);
252 if (tracediff[nmul - 1] > 0) {
255 D = -ONE * X * X + TWO * D;
262 D.frob_thresh(frob_trunc);
263 idemerror[nmul] = Tmatrix::frob_diff(D, X);
265 tracediff[nmul] = D.trace() - nocc;
267 nmul_firstpart = nmul;
271 if (nmul + 1 >= maxmul) {
274 "Purification reached maxmul"
275 " without convergence", maxmul);
277 if (tracediff[nmul] > 0) {
279 idemerror[nmul] = Tmatrix::frob_diff(D, X);
283 tracediff[nmul] = D.trace() - nocc;
285 D = -ONE * X * X + TWO * D;
286 idemerror[nmul] = Tmatrix::frob_diff(D, X);
289 tracediff[nmul] = D.trace() - nocc;
293 X = -ONE * D * D + TWO * X;
294 idemerror[nmul] = Tmatrix::frob_diff(D, X);
297 tracediff[nmul] = X.trace() - nocc;
300 idemerror[nmul] = Tmatrix::frob_diff(D, X);
303 tracediff[nmul] = D.trace() - nocc;
305 D.frob_thresh(frob_trunc);
307 }
while (idemerror[nmul - 1] > frob_trunc);
314 template<
typename Treal,
typename Tmatrix>
316 Fun const fermifun(polys, nmul, 0.5);
317 Treal chempot =
bisection(fermifun, (Treal)0, (Treal)1, tol);
318 return (lmin - lmax) * chempot + lmax;
321 template<
typename Treal,
typename Tmatrix>
325 for (
int mul = nmul_firstpart; mul < nmul; mul++) {
326 if (idemerror[mul] < 1.0 / 4) {
328 tmp =
bisection(homofun, (Treal)0, (Treal)1, tol);
333 homo = tmp > homo ? tmp : homo;
336 return (lmin - lmax) * homo + lmax;
339 template<
typename Treal,
typename Tmatrix>
343 for (
int mul = nmul_firstpart; mul < nmul; mul++) {
344 if (idemerror[mul] < 1.0 / 4) {
346 tmp =
bisection(lumofun, (Treal)0, (Treal)1, tol);
351 lumo = tmp < lumo ? tmp : lumo;
354 return (lmin - lmax) * lumo + lmax;
357 template<
typename Treal,
typename Tmatrix>
359 for (
int ind = start; ind < stop; ind ++) {
360 std::cout <<
"Iteration: " << ind
361 <<
" Idempotency error: " << idemerror[ind]
362 <<
" Tracediff: " << tracediff[ind]
363 <<
" Poly: " << polys[ind]