36 #ifndef MAT_PURIFICATION
37 #define MAT_PURIFICATION
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.");
void setActualThresh(Treal const thr)
Definition: PuriStepInfo.h:85
Tmatrix & X
Definition: Purification.h:61
Treal getEigAccLoss() const
Definition: PuriStepInfo.h:140
Interval< Treal > getHomoF() const
Returns the best interval containing the homo eigenvalue (not transformed to [0, 1]) ...
Definition: PuriInfo.h:112
Treal upp() const
Definition: Interval.h:143
void computeExactValues(Tmatrix const &X, Tmatrix const &X2)
Definition: PuriStepInfo.h:152
void purify()
Definition: Purification.h:186
void setTimeThresh(float t)
Definition: PuriStepInfo.h:164
void improveEigInterval(Interval< Treal > const eInt)
Improve eigenvalue bounds and delta if possible.
Definition: PuriStepInfo.h:436
void tic()
Definition: matInclude.cc:77
static Interval intersect(Interval const &A, Interval const &B)
Definition: Interval.h:51
Interval< Treal > getLumoF() const
Returns the best interval containing the lumo eigenvalue (not transformed to [0, 1]) ...
Definition: PuriInfo.h:116
PuriInfo< Treal, VectorType, TdebugPolicy > & info
Definition: Purification.h:65
Interval< Treal > const & getLumo() const
Definition: PuriStepInfo.h:118
void improveHomoLumo(Interval< Treal > const homoInt, Interval< Treal > const lumoInt)
Definition: PuriStepInfo.h:310
void increase(Treal const value)
Increases interval with value in both directions.
Definition: Interval.h:131
bool empty() const
Definition: Interval.h:49
void setTimeXmX2Norm(float t)
Definition: PuriStepInfo.h:166
Definition: matInclude.h:158
float toc()
Definition: matInclude.cc:81
Definition: Interval.h:44
void setTimeTotal(float t)
Definition: PuriStepInfo.h:167
Definition: matInclude.h:135
void setChosenThresh(Treal const thr)
Definition: PuriStepInfo.h:83
void setTraceX(Treal const trX)
Definition: PuriStepInfo.h:91
void setTimeSquare(float t)
Definition: PuriStepInfo.h:165
Treal low() const
Definition: Interval.h:142
Tmatrix::VectorType VectorType
Definition: Purification.h:45
Definition: Purification.h:43
PuriStepInfo< Treal, Tvector, TdebugPolicy > & getNext()
Definition: PuriInfo.h:84
Interval< Treal > const & getHomo() const
Definition: PuriStepInfo.h:115
normType const normXmX2
Definition: Purification.h:64
void stepComputeInfo(PuriStepInfo< Treal, VectorType, TdebugPolicy > ¤tStep)
Definition: Purification.h:203
void setTraceX2(Treal const trX2)
Definition: PuriStepInfo.h:92
void setNnzX(size_t const nzX)
Definition: PuriStepInfo.h:135
void setPoly()
Definition: PuriStepInfo.h:292
void setEigVecPtr(Tvector *eigVecPtr_)
Improves homo and lumo bounds if the new ones are better.
Definition: PuriStepInfo.h:107
Treal getOptimalThresh() const
Definition: PuriInfo.h:363
Definition: matInclude.h:135
Tmatrix X2
Definition: Purification.h:62
void step()
Definition: Purification.h:121
void checkIntervals(const char *descriptionString) const
Definition: PuriStepInfo.h:147
void setNnzX2(size_t const nzX2)
Definition: PuriStepInfo.h:137
int niter
Definition: Purification.h:66
Purification(Tmatrix &M, normType const normXmX2, PuriInfo< Treal, VectorType, TdebugPolicy > &info)
Constructor.
Definition: Purification.h:75
void computeEigAccLoss()
Computes a probable upper bound of the accuracy that is lost in the eigenvalues of X * X because of l...
Definition: PuriStepInfo.h:467
void setXmX2EuclNorm(Interval< Treal > const XmX2eucl)
Sets XmX2EuclNorm bounds.
Definition: PuriStepInfo.h:98
Interval< Treal > const & getEigInterval() const
Definition: PuriStepInfo.h:121
void setMemUsageBeforeTrunc()
Definition: PuriStepInfo.h:158
Treal template_blas_sqrt(Treal x)
normType const normTruncation
Definition: Purification.h:63
normType
Definition: matInclude.h:135
Interval< Treal > getEigFInterval() const
Definition: PuriInfo.h:81