46 template<
typename Treal,
typename Tvector,
typename TdebugPolicy>
53 Treal toleratedEigenvalError,
54 Treal toleratedSubspaceError,
63 for (
int ind = 0; ind <
maxSteps; ++ind)
121 resultVector.resize(
nSteps);
122 for(
int i = 0; i <
nSteps; i++)
123 resultVector[i] =
step[i].getPoly();
127 resultVector.resize(
nSteps);
128 for(
int i = 0; i <
nSteps; i++)
129 resultVector[i] =
step[i].getChosenThresh();
145 void mTimings(std::ostream & file)
const;
146 void mInfo(std::ostream & file)
const;
147 void mMemUsage(std::ostream & file)
const;
154 for(
int ind = 0; ind <
nSteps; ++ind)
164 for(
int ind = 0; ind <
nSteps; ++ind)
174 Tvector & eigVecHOMO)
const;
220 template<
typename Treal,
typename Tvector,
typename TdebugPolicy>
222 if ( step[nSteps-1].getCorrectOccupation() )
224 step[nSteps-1].setCorrectOccupation();
225 correct_occupation_was_forced_flag =
true;
229 template<
typename Treal,
typename Tvector,
typename TdebugPolicy>
231 if (step[nSteps-1].getCorrectOccupation()) {
235 for (
int ind = 0; ind < nSteps; ind++) {
236 Treal XmX2Eucl = step[ind].getXmX2EuclNorm().upp();
237 if ( XmX2Eucl < 1 / (Treal)4) {
241 while (!middleInt.
empty() && i < nSteps - 1) {
242 middleInt.
puriStep(step[i].getPoly());
243 middleInt.
decrease(step[i+1].getActualThresh());
248 if (middleInt.
cover(0.5))
249 step[ind].setCorrectOccupation();
255 template<
typename Treal,
typename Tvector,
typename TdebugPolicy>
257 for (
int ind = nSteps - 2; ind >= 0; ind--) {
258 step[ind].exchangeInfoWithNext(step[ind + 1]);
260 for (
int ind = 0; ind < nSteps - 1; ind++) {
261 step[ind].exchangeInfoWithNext(step[ind + 1]);
265 template<
typename Treal,
typename Tvector,
typename TdebugPolicy>
268 for (
int ind = 0; ind < end; ind++) {
269 error += step[ind].subspaceError();
274 template<
typename Treal,
typename Tvector,
typename TdebugPolicy>
281 Treal tolError = tolSubspaceError - subspaceError(nSteps - 1);
282 Treal initialAltThresh = 1;
293 Treal lastThresh = 0;
294 if (nSteps == 0 || nSteps == 1) {
295 Treal lmax = eigFInterval.upp();
296 Treal lmin = eigFInterval.low();
299 (homoF.upp() - lmax) / (lmin - lmax));
303 step[nSteps - 2].getHomo().low());
304 initialAltThresh = getThreshIncreasingGap(initialGap);
305 lastThresh = step[nSteps - 2].getChosenThresh();
306 initialAltThresh = initialAltThresh > lastThresh / 10 ?
307 initialAltThresh : lastThresh / 10;
308 initialGap.
puriStep(step[nSteps - 2].getPoly());
310 if (initialGap.
empty())
314 if (1.0 - initialGap.
upp() < tolEigenvalError &&
315 initialGap.
low() < tolEigenvalError) {
321 Treal thetaPerStep = 0.0;
325 Treal currThresh = 0;
326 int stepsLeftPrev = -2;
328 while (stepsLeft > stepsLeftPrev) {
330 altThresh = initialAltThresh;
331 currThresh = thetaPerStep * gap.
length() / (1 + thetaPerStep);
332 currThresh = currThresh < altThresh ? currThresh : altThresh;
334 lastThresh = currThresh;
335 stepsLeftPrev = stepsLeft;
337 while (!gap.
empty() &&
338 (1.0 - gap.
upp() > tolEigenvalError ||
339 gap.
low() > tolEigenvalError)) {
340 altThresh = getThreshIncreasingGap(gap);
341 altThresh = altThresh > lastThresh / 10 ? altThresh : lastThresh / 10;
345 currThresh = thetaPerStep * gap.
length() / (1 + thetaPerStep);
346 currThresh = currThresh < altThresh ? currThresh : altThresh;
348 lastThresh = currThresh;
351 thetaPerStep = tolError / (stepsLeft + 1);
355 thetaPerStep = tolError / (stepsLeftPrev + 1);
356 thresh = thetaPerStep * initialGap.
length() / (1 + thetaPerStep);
362 template<
typename Treal,
typename Tvector,
typename TdebugPolicy>
366 estimateStepsLeft(stepsLeft, thresh);
368 step[nSteps - 1].setEstimatedStepsLeft(stepsLeft);
370 thresh = tolSubspaceError / 100;
374 step[nSteps - 2].getHomo().low());
375 if (!middleGap.
empty()) {
377 Treal altThresh = getThreshIncreasingGap(middleGap);
379 Treal lastThresh = step[nSteps - 2].getChosenThresh();
380 altThresh = altThresh > lastThresh / 10 ? altThresh : lastThresh / 10;
381 thresh = thresh < altThresh ? thresh : altThresh;
384 Treal tmp = 1 - xmX2Eucl.
upp() * 4.0;
387 thresh = thresh < altThresh ? thresh : altThresh;
397 template<
typename Treal,
typename Tvector,
typename TdebugPolicy>
401 Treal delta = middleGap.
upp() - x;
439 template<
typename Treal,
typename Tvector,
typename TdebugPolicy>
443 if (nSteps == 0 || nSteps == 1) {
447 Treal ep = 0.207106781186547;
458 homoCopy.
puriStep(step[nSteps - 2].getPoly());
459 lumoCopy.
puriStep(step[nSteps - 2].getPoly());
464 howAccurate = step[nSteps - 1].getChosenThresh() / 100;
465 Treal highestPossibleAccuracy = 10.0 * step[nSteps - 1].getEigAccLoss();
468 howAccurate = howAccurate > highestPossibleAccuracy ?
469 howAccurate : highestPossibleAccuracy;
471 if (homoCopy.
length() > 0.2 || lumoCopy.
length() > 0.2) {
473 if (step[nSteps - 2].getN0() / (n - nocc) > 0.5 &&
474 step[nSteps - 2].getN1() / nocc > 0.5 &&
475 step[nSteps - 2].getXmX2EuclNorm().midPoint() < ep)
482 bool computeHomo =
true;
483 bool computeLumo =
true;
484 if (homoIsAccuratelyKnown() ||
485 homoCopy.
upp() < 0.5 ||
488 if (lumoIsAccuratelyKnown() ||
489 lumoCopy.
low() > 0.5 ||
492 return computeHomo || computeLumo;
497 template<
typename Treal,
typename Tvector,
typename TdebugPolicy>
500 Treal lmax = eigFInterval.upp();
501 Treal lmin = eigFInterval.low();
504 homoF.intersect(altHomo);
505 lumoF.intersect(altLumo);
508 template<
typename Treal,
typename Tvector,
typename TdebugPolicy>
511 for (
int ind = 0; ind < nSteps; ++ind)
512 accTime += (
double)step[ind].getTimeSquare();
515 template<
typename Treal,
typename Tvector,
typename TdebugPolicy>
518 for (
int ind = 0; ind < nSteps; ++ind)
519 accTime += (
double)step[ind].getTimeThresh();
522 template<
typename Treal,
typename Tvector,
typename TdebugPolicy>
525 for (
int ind = 0; ind < nSteps; ++ind)
526 accTime += (
double)step[ind].getTimeXmX2Norm();
529 template<
typename Treal,
typename Tvector,
typename TdebugPolicy>
532 for (
int ind = 0; ind < nSteps; ++ind)
533 accTime += (
double)step[ind].getTimeTotal();
538 template<
typename Treal,
typename Tvector,
typename TdebugPolicy>
542 for (
int ind = 0; ind < nSteps; ++ind) {
544 step[ind].getTimeSquare()<<
" "<<
545 step[ind].getTimeThresh()<<
" "<<
546 step[ind].getTimeXmX2Norm()<<
" "<<
547 step[ind].getTimeTotal()<<
" "<<
548 step[ind].getTimeXX2Write()<<
" "<<
549 step[ind].getTimeXX2Read()<<
" "<<
552 s<<
"];\n"<<
"figure; bar(puriTime(:,1:3),'stacked')"<<std::endl<<
553 "legend('Matrix Square', 'Truncation', '||X-X^2||'),"<<
554 " xlabel('Iteration'), ylabel('Time (seconds)')"<<std::endl;
555 s<<
"figure; plot(puriTime(:,4),'-x')"<<std::endl;
558 template<
typename Treal,
typename Tvector,
typename TdebugPolicy>
561 s<<
"% PURIFICATION INFO IN MATLAB/OCTAVE FILE FORMAT"<<std::endl;
564 s<<
"n = "<<n<<
";"<<std::endl;
565 s<<
"nocc = "<<nocc<<
";"<<std::endl;
566 s<<
"tolSubspaceError = "<<tolSubspaceError<<
";"<<std::endl;
567 s<<
"tolEigenvalError = "<<tolEigenvalError<<
";"<<std::endl;
568 s<<
"correct_occupation_was_forced_flag = "
569 <<correct_occupation_was_forced_flag<<
";"<<std::endl;
570 s<<
"% Each row of the following matrix contains purification info \n"
571 <<
"% for one step.\n"
572 <<
"% The columns are arranged as: \n"
573 <<
"% ind, n0, n1, trace(X), trace(X*X), poly, actualThresh, delta, "
574 <<
"correctOcc, XmX2low, XmX2upp, HOMOmid, LUMOmid, "
575 <<
"nnz(X), nnz(X*X), chosenThresh, estimatedStepsLeft "
578 for (
int ind = 0; ind < nSteps; ++ind) {
581 step[ind].getN0()<<
" "<<
582 step[ind].getN1()<<
" "<<
583 step[ind].getTraceX()<<
" "<<
584 step[ind].getTraceX2()<<
" "<<
585 step[ind].getPoly()<<
" "<<
586 step[ind].getActualThresh()<<
" "<<
587 step[ind].getDelta()<<
" "<<
588 step[ind].getCorrectOccupation()<<
" "<<
589 step[ind].getXmX2EuclNorm().low()<<
" "<<
590 step[ind].getXmX2EuclNorm().upp()<<
" "<<
591 step[ind].getHomo().midPoint()<<
" "<<
592 step[ind].getLumo().midPoint()<<
" "<<
593 step[ind].getNnzX()<<
" "<<
594 step[ind].getNnzX2()<<
" "<<
595 step[ind].getChosenThresh()<<
" "<<
596 step[ind].getEstimatedStepsLeft()<<
" "<<
601 s<<
"ind = puriMat(:,1);\n";
603 <<
"plot(ind,puriMat(:,12),'x-b',ind,puriMat(:,13),'o-r')\n"
604 <<
"legend('HOMO','LUMO'),xlabel('Iteration')\n"
605 <<
"axis([0 max(ind) 0 1])\n";
607 <<
"plot(ind,100*puriMat(:,2)/(n-nocc),'o-r',...\n"
608 <<
"ind, 100*puriMat(:,3)/nocc,'x-b')\n"
609 <<
"legend('N0','N1'),xlabel('Iteration'), ylabel('Percentage')\n"
610 <<
"axis([0 max(ind) 0 100])\n";
613 <<
"plot(ind,100*puriMat(:,14)/(n*n),'o-r',...\n"
614 <<
"ind, 100*puriMat(:,15)/(n*n),'x-b')\n"
615 <<
"legend('nnz(X)','nnz(X^2)'),xlabel('Iteration') \n"
616 <<
"ylabel('Percentage')\n"
617 <<
"axis([0 max(ind) 0 100])\n";
619 <<
"semilogy(ind,puriMat(:,16),'x-r',ind,puriMat(:,7),'o-b')\n"
620 <<
"xlabel('Iteration'), ylabel('Threshold')\n"
621 <<
"legend('Chosen threshold', 'Actual threshold')\n"
622 <<
"axis([0 max(ind) min(puriMat(:,7))/10 max(puriMat(:,16))*10])\n";
624 <<
"plot(ind,ind(end:-1:1)-1,ind,puriMat(:,17),'x-r')\n"
625 <<
"legend('Steps left', 'Estimated steps left')\n"
626 <<
"xlabel('Iteration')\n";
631 template<
typename Treal,
typename Tvector,
typename TdebugPolicy>
634 s<<
"puriMemUsage = [";
635 for (
int ind = 0; ind < nSteps; ++ind) {
639 memUsageBeforeTrunc.
res <<
" "<<
640 memUsageBeforeTrunc.
virt<<
" "<<
641 memUsageBeforeTrunc.
peak<<
" "<<
642 memUsageInXmX2Diff.
res <<
" "<<
643 memUsageInXmX2Diff.
virt<<
" "<<
644 memUsageInXmX2Diff.
peak<<
" "<<
649 <<
"plot(puriMemUsage(:,1),'x-b')\n"
651 <<
"plot(puriMemUsage(:,2),'o-r')\n"
652 <<
"plot(puriMemUsage(:,3),'^-g')\n"
653 <<
"legend('Resident','Virtual','Peak'),xlabel('Iteration'),ylabel('Mem usage before trunc [GB]')\n"
654 <<
"%force y axis to start at 0\n"
655 <<
"axissaved = axis; axissaved(3) = 0; axis(axissaved);\n"
658 <<
"plot(puriMemUsage(:,4),'x-b')\n"
660 <<
"plot(puriMemUsage(:,5),'o-r')\n"
661 <<
"plot(puriMemUsage(:,6),'^-g')\n"
662 <<
"legend('Resident','Virtual','Peak'),xlabel('Iteration'),ylabel('Mem usage in XmX2Diff [GB]')\n"
663 <<
"%force y axis to start at 0\n"
664 <<
"axissaved = axis; axissaved(3) = 0; axis(axissaved);\n"
668 template<
typename Treal,
typename Tvector,
typename TdebugPolicy>
671 Tvector & eigVecHOMO)
const {
674 for (
int ind = 0; ind < nSteps; ++ind) {
675 if (!haveHOMO && step[ind].getHomoWasComputed()) {
676 eigVecHOMO = *step[ind].getEigVecPtr();
679 if (!haveLUMO && step[ind].getLumoWasComputed()) {
680 eigVecLUMO = *step[ind].getEigVecPtr();
Interval< Treal > getHomoF() const
Returns the best interval containing the homo eigenvalue (not transformed to [0, 1]) ...
Definition: PuriInfo.h:112
std::string getNormTypeString(normType nType)
Definition: matInclude.cc:55
Treal upp() const
Definition: Interval.h:143
Interval< Treal > homoF
Interval containing the HOMO eigenvalue before transformation to [0, 1].
Definition: PuriInfo.h:193
bool ShouldComputeXmX2EuclNormAccurately(Treal &howAccurate) const
Definition: PuriInfo.h:441
Treal length() const
Returns the length of the interval.
Definition: Interval.h:107
void decrease(Treal const value)
Definition: Interval.h:138
float peak
Definition: matInclude.h:174
void improveCorrectOccupation()
Improves the correct occupation indicator Call AFTER convergence and ONLY if it is known that the occ...
Definition: PuriInfo.h:230
PuriStepInfo< Treal, Tvector, TdebugPolicy > * step
Definition: PuriInfo.h:180
void mInfo(std::ostream &file) const
Definition: PuriInfo.h:560
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
static ergo_real distance(const ergo_real *a, const ergo_real *b)
Coomputes distance between two points, they do not need to be of the Vector3D type.
Definition: dft_common.cc:530
void getPolys(std::vector< int > &resultVector)
Definition: PuriInfo.h:120
PuriStepInfo< Treal, Tvector, TdebugPolicy > & operator()(int const ind)
Definition: PuriInfo.h:89
double getAccumulatedTimeXmX2Norm() const
Definition: PuriInfo.h:523
bool empty() const
Definition: Interval.h:49
#define ASSERTALWAYS(x)
Definition: DebugPolicies.h:83
Definition: allocate.cc:30
int nocc
Number of occupied orbitals.
Definition: PuriInfo.h:178
Treal midPoint() const
Definition: Interval.h:113
void mTimings(std::ostream &file) const
Definition: PuriInfo.h:540
Treal template_blas_fabs(Treal x)
Definition: Interval.h:44
void estimateStepsLeft(int &stepsLeft, Treal &thresh) const
Estimates the number of steps (upper bound) needed for convergence based on homo/lumo information and...
Definition: PuriInfo.h:276
bool homoIsAccuratelyKnown() const
HOMO estimation is considered to be accurate if the error is smaller than tolSubspaceError / 100 in s...
Definition: PuriInfo.h:152
int getNSteps() const
Definition: PuriInfo.h:118
void puriStep(int poly)
Definition: Interval.h:228
void mMemUsage(std::ostream &file) const
Definition: PuriInfo.h:633
void forceCorrectOccupation()
Set the correctOccupation flag in the current step to 1.
Definition: PuriInfo.h:221
bool correct_occupation_was_forced() const
Definition: PuriInfo.h:134
Contains information about the matrix before a purification step and about the step.
Definition: PuriStepInfo.h:55
float res
Definition: matInclude.h:172
Contains information about a purification process.
Definition: PuriInfo.h:47
normType getNormForTruncation() const
Definition: PuriInfo.h:169
Treal low() const
Definition: Interval.h:142
bool cover(Treal const value) const
Definition: Interval.h:117
float virt
Definition: matInclude.h:173
Interval< Treal > const eigFInterval
Interval containing all eigenvalues before transformation to the [0, 1] interval. ...
Definition: PuriInfo.h:188
PuriStepInfo< Treal, Tvector, TdebugPolicy > & getNext()
Definition: PuriInfo.h:84
double getAccumulatedTimeTotal() const
Definition: PuriInfo.h:530
Treal getThreshIncreasingGap(Interval< Treal > const &middleGap) const
Definition: PuriInfo.h:399
bool correct_occupation_was_forced_flag
Correct occupation was assumed, not guaranteed.
Definition: PuriInfo.h:186
int getMaxSteps() const
Definition: PuriInfo.h:117
Interval< Treal > lumoF
Interval containing the LUMO eigenvalue before transformation to [0, 1].
Definition: PuriInfo.h:198
Treal subspaceError() const
Returns the subspace error introduced so far.
Definition: PuriInfo.h:97
void getHOMOandLUMOeigVecs(Tvector &eigVecLUMO, Tvector &eigVecHOMO) const
Definition: PuriInfo.h:670
int nSteps
Number of taken steps.
Definition: PuriInfo.h:183
bool converged()
Definition: PuriInfo.h:138
bool lumoIsAccuratelyKnown() const
LUMO estimation is considered to be accurate if the error is smaller than tolSubspaceError / 100 in s...
Definition: PuriInfo.h:162
Treal getOptimalThresh() const
Definition: PuriInfo.h:363
int n
System size.
Definition: PuriInfo.h:177
Treal const tolEigenvalError
Tolerated error in eigenvalues at convergence.
Definition: PuriInfo.h:207
Definition: matInclude.h:171
int maxSteps
Capacity of step array.
Definition: PuriInfo.h:181
void improveHomoLumoF()
Definition: PuriInfo.h:499
double getAccumulatedTimeThresh() const
Definition: PuriInfo.h:516
normType const normForTruncation
Norm used for truncation of small matrix elements.
Definition: PuriInfo.h:210
Treal const tolSubspaceError
Tolerated error in angle between correct and computed subspace.
Definition: PuriInfo.h:203
void improveInfo()
Improve homo / lumo values in each step.
Definition: PuriInfo.h:256
PuriInfo(int nn, int noc, Interval< Treal > eigFInt, Interval< Treal > hoF, Interval< Treal > luF, Treal toleratedEigenvalError, Treal toleratedSubspaceError, normType normForTruncation_, int maxS=100)
Definition: PuriInfo.h:49
void getThreshValues(std::vector< Treal > &resultVector)
Definition: PuriInfo.h:126
Treal template_blas_sqrt(Treal x)
normType
Definition: matInclude.h:135
double getAccumulatedTimeSquare() const
Definition: PuriInfo.h:509
virtual ~PuriInfo()
Definition: PuriInfo.h:67
Interval< Treal > getEigFInterval() const
Definition: PuriInfo.h:81