36 #ifndef MAT_PURISTEPINFO
37 #define MAT_PURISTEPINFO
41 #define __ID__ "$Id: PuriStepInfo.h 4437 2012-07-05 09:01:18Z elias $"
54 template<
typename Treal,
typename Tvector,
typename TdebugPolicy>
57 explicit PuriStepInfo(
int nn = -1,
int noc = -1, Treal eigvalConvCrit = 0.0)
80 return homoLumo_converged || XmX2norm_converged;
151 template<
typename Tmatrix>
185 return homo.length() < accuracyLimit;
192 return lumo.length() < accuracyLimit;
269 template<
typename Treal,
typename Tvector,
typename TdebugPolicy>
280 <<
" n0: "<<psi.
getN0()
281 <<
" n1: "<<psi.
getN1()
291 template<
typename Treal,
typename Tvector,
typename TdebugPolicy>
295 if (homo.low() + lumo.upp() > 1)
308 template<
typename Treal,
typename Tvector,
typename TdebugPolicy>
312 checkIntervals(
"PuriStepInfo::improveHomoLumo A0");
313 homo.intersect(homoInt);
314 checkIntervals(
"PuriStepInfo::improveHomoLumo A1");
315 lumo.intersect(lumoInt);
316 checkIntervals(
"PuriStepInfo::improveHomoLumo A2");
319 if (homo.low() > 0.5 && lumo.upp() < 0.5)
320 this->setCorrectOccupation();
322 if (correctOccupation && 1.0 - XmX2EuclNorm.upp() * 4.0 > 0) {
324 sqrtInt((XmX2EuclNorm * (Treal)(-4.0)) + (Treal)1.0);
336 checkIntervals(
"PuriStepInfo::improveHomoLumo B");
339 if (!eigInterval.cover((1 +
template_blas_sqrt(1.0 + 4.0 * XmX2EuclNorm.low())) / 2) &&
343 (tmp + (Treal)1.0) / (Treal)2.0;
345 ((tmp * (Treal)(-1.0)) + (Treal)1.0) / (Treal)2.0;
352 lumo.intersect(lumoTmp);
358 homo.intersect(homoTmp);
360 checkIntervals(
"PuriStepInfo::improveHomoLumo C");
365 template<
typename Treal,
typename Tvector,
typename TdebugPolicy>
420 this->improveHomoLumo(homoTmp, lumoTmp);
423 template<
typename Treal,
typename Tvector,
typename TdebugPolicy>
426 if (actualThresh >= gap.
length())
429 Treal error = actualThresh / (gap.
length() - actualThresh);
430 return error < 1.0 ? error : (Treal)1.0;
434 template<
typename Treal,
typename Tvector,
typename TdebugPolicy>
437 eigInterval.intersect(eInt);
438 Treal delta1 = eigInterval.upp() - 1;
439 Treal delta0 = -eigInterval.low();
440 delta = delta1 > delta0 ? delta1 : delta0;
445 template<
typename Treal,
typename Tvector,
typename TdebugPolicy>
449 if (XmX2EuclNorm.upp() < 1 / (Treal)4)
451 n1 = (traceX2 - delta * (1 - beta) * n -
452 (1 - delta - beta) * traceX) /
453 ((1 + 2 * delta) * (delta + beta));
454 n0 = (traceX2 + beta * (1 + delta) * n -
455 (1 + delta + beta) * traceX) /
456 ((1 + 2 * delta) * (delta + beta));
459 correctOccupation = 1;
466 template<
typename Treal,
typename Tvector,
typename TdebugPolicy>
468 Treal nnzPerRowX = nnzX / (Treal)n;
469 Treal maxAbsErrPerElX2 = getRelPrecision<Treal>() * nnzPerRowX;