36 #ifndef MAT_PURIFICATION
37 #define MAT_PURIFICATION
40 #define __ID__ "$Id: Purification.h 4437 2012-07-05 09:01:18Z elias $"
42 template<
typename Treal,
typename Tmatrix,
typename TdebugPolicy>
73 template<
typename Treal,
typename Tmatrix,
typename TdebugPolicy>
78 :X(M), X2(M), normTruncation( infop.getNormForTruncation() ), normXmX2(normXmX2_inp), info(infop),niter(0) {
89 X.add_identity(-lmax);
90 X *= ((Treal)1.0 / (lmin - lmax));
94 homo = (homo - lmax) / (lmin - lmax);
95 lumo = (lumo - lmax) / (lmin - lmax);
107 throw Failure(
"Purification<Treal, Tmatrix, TdebugPolicy>::"
108 "Purification(Tmatrix&, normType const, "
109 "PuriInfo<Treal, VectorType, TdebugPolicy>&): "
110 "HOMO or LUMO empty in purification constructor. ");
120 template<
typename Treal,
typename Tmatrix,
typename TdebugPolicy>
138 if (info(niter).getPoly()) {
141 X2 += (Treal)2.0 * X;
155 Treal chosenThresh = info.getOptimalThresh();
163 stepComputeInfo(currentStep);
165 currentStep.
getHomo().low() > 0.9 &&
166 currentStep.
getLumo().upp() < 0.1 &&
167 info(niter-5).getHomo().low() > 0.9 &&
168 info(niter-5).getLumo().upp() < 0.1 &&
169 ((1 - currentStep.
getHomo().low()) >
170 (1 - info(niter-5).getHomo().low()) -
173 info(niter-5).getLumo().upp() -
175 throw Failure(
"Purification<Treal, Tmatrix, TdebugPolicy>"
176 "::step(): HOMO-LUMO gap has not increased in the "
177 "last five iterations. Desired accuracy too high for"
178 "current floating point precision?");
185 template<
typename Treal,
typename Tmatrix,
typename TdebugPolicy>
188 while (niter < info.getMaxSteps() - 1 && !info.converged() ) {
191 if ( info.converged() ) {
194 info.forceCorrectOccupation();
195 info.improveCorrectOccupation();
197 info.improveHomoLumoF();
201 template<
typename Treal,
typename Tmatrix,
typename TdebugPolicy>
204 Treal
const XmX2ENIsSmallValue = 0.207106781186547;
205 Treal
const ONE = 1.0;
224 typename Tmatrix::VectorType * eigVecPtr =
new typename Tmatrix::VectorType;
228 if (info.ShouldComputeXmX2EuclNormAccurately(diffAcc)) {
230 XmX2EN = Tmatrix::euclDiffIfSmall(X, X2,
235 XmX2EN = Tmatrix::diffIfSmall(X, X2,
240 XmX2EN = Tmatrix::diffIfSmall(X, X2,
246 if (!eigVecPtr->is_empty())
258 currentStep.
checkIntervals(
"Purification::stepComputeInfo after setXmX2EuclNorm.");
268 throw Failure(
"Purification<Treal, Tmatrix, TdebugPolicy>"
270 "Eigenvalues moved to far from [0, 1] interval.");