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>
229 template<
typename Treal,
typename Tvector,
typename TdebugPolicy>
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) {
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>
282 Treal initialAltThresh = 1;
293 Treal lastThresh = 0;
299 (
homoF.upp() - lmax) / (lmin - lmax));
305 lastThresh =
step[
nSteps - 2].getChosenThresh();
306 initialAltThresh = initialAltThresh > lastThresh / 10 ?
307 initialAltThresh : lastThresh / 10;
310 if (initialGap.
empty())
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() &&
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>
368 step[
nSteps - 1].setEstimatedStepsLeft(stepsLeft);
375 if (!middleGap.
empty()) {
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>
447 Treal ep = 0.207106781186547;
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) {
475 step[
nSteps - 2].getXmX2EuclNorm().midPoint() < ep)
482 bool computeHomo =
true;
483 bool computeLumo =
true;
485 homoCopy.
upp() < 0.5 ||
489 lumoCopy.
low() > 0.5 ||
492 return computeHomo || computeLumo;
497 template<
typename Treal,
typename Tvector,
typename TdebugPolicy>
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;
568 s<<
"correct_occupation_was_forced_flag = " 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();
Treal upp() const
Definition: Interval.h:143
bool lumoIsAccuratelyKnown() const
LUMO estimation is considered to be accurate if the error is smaller than tolSubspaceError / 100 in s...
Definition: PuriInfo.h:162
double getAccumulatedTimeXmX2Norm() const
Definition: PuriInfo.h:523
void getHOMOandLUMOeigVecs(Tvector &eigVecLUMO, Tvector &eigVecHOMO) const
Definition: PuriInfo.h:670
double getAccumulatedTimeSquare() const
Definition: PuriInfo.h:509
std::string getNormTypeString(normType nType)
Definition: matInclude.cc:55
bool homoIsAccuratelyKnown() const
HOMO estimation is considered to be accurate if the error is smaller than tolSubspaceError / 100 in s...
Definition: PuriInfo.h:152
Interval< Treal > homoF
Interval containing the HOMO eigenvalue before transformation to [0, 1].
Definition: PuriInfo.h:193
double getAccumulatedTimeTotal() const
Definition: PuriInfo.h:530
void decrease(Treal const value)
Definition: Interval.h:138
Treal low() const
Definition: Interval.h:142
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
static Interval intersect(Interval const &A, Interval const &B)
Definition: Interval.h:51
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 getAccumulatedTimeThresh() const
Definition: PuriInfo.h:516
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
#define ASSERTALWAYS(x)
Definition: DebugPolicies.h:83
Definition: allocate.cc:30
int nocc
Number of occupied orbitals.
Definition: PuriInfo.h:178
Treal template_blas_fabs(Treal x)
Definition: Interval.h:44
void mTimings(std::ostream &file) const
Definition: PuriInfo.h:540
void puriStep(int poly)
Definition: Interval.h:228
Treal getOptimalThresh() const
Definition: PuriInfo.h:363
void forceCorrectOccupation()
Set the correctOccupation flag in the current step to 1.
Definition: PuriInfo.h:221
Contains information about the matrix before a purification step and about the step.
Definition: PuriStepInfo.h:55
float res
Definition: matInclude.h:172
Interval< Treal > getHomoF() const
Returns the best interval containing the homo eigenvalue (not transformed to [0, 1]) ...
Definition: PuriInfo.h:112
void mMemUsage(std::ostream &file) const
Definition: PuriInfo.h:633
Contains information about a purification process.
Definition: PuriInfo.h:47
void mInfo(std::ostream &file) const
Definition: PuriInfo.h:560
float virt
Definition: matInclude.h:173
Treal midPoint() const
Definition: Interval.h:113
Treal getThreshIncreasingGap(Interval< Treal > const &middleGap) const
Definition: PuriInfo.h:399
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
Treal subspaceError() const
Returns the subspace error introduced so far.
Definition: PuriInfo.h:97
bool correct_occupation_was_forced_flag
Correct occupation was assumed, not guaranteed.
Definition: PuriInfo.h:186
bool ShouldComputeXmX2EuclNormAccurately(Treal &howAccurate) const
Definition: PuriInfo.h:441
Interval< Treal > lumoF
Interval containing the LUMO eigenvalue before transformation to [0, 1].
Definition: PuriInfo.h:198
Interval< Treal > getLumoF() const
Returns the best interval containing the lumo eigenvalue (not transformed to [0, 1]) ...
Definition: PuriInfo.h:116
bool cover(Treal const value) const
Definition: Interval.h:117
int nSteps
Number of taken steps.
Definition: PuriInfo.h:183
bool converged()
Definition: PuriInfo.h:138
Treal length() const
Returns the length of the interval.
Definition: Interval.h:107
int n
System size.
Definition: PuriInfo.h:177
bool empty() const
Definition: Interval.h:49
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
bool correct_occupation_was_forced() const
Definition: PuriInfo.h:134
normType getNormForTruncation() const
Definition: PuriInfo.h:169
void improveHomoLumoF()
Definition: PuriInfo.h:499
Interval< Treal > getEigFInterval() const
Definition: PuriInfo.h:81
normType const normForTruncation
Norm used for truncation of small matrix elements.
Definition: PuriInfo.h:210
int getNSteps() const
Definition: PuriInfo.h:118
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
int getMaxSteps() const
Definition: PuriInfo.h:117
virtual ~PuriInfo()
Definition: PuriInfo.h:67