41 #define __ID__ "$Id: PuriInfo.h 4437 2012-07-05 09:01:18Z elias $"
46 template<
typename Treal,
typename Tvector,
typename TdebugPolicy>
53 Treal toleratedEigenvalError,
54 Treal toleratedSubspaceError,
63 for (
int ind = 0; ind <
maxSteps; ++ind)
133 void mTimings(std::ostream & file)
const;
134 void mInfo(std::ostream & file)
const;
135 void mMemUsage(std::ostream & file)
const;
142 for(
int ind = 0; ind <
nSteps; ++ind)
152 for(
int ind = 0; ind <
nSteps; ++ind)
162 Tvector & eigVecHOMO)
const;
208 template<
typename Treal,
typename Tvector,
typename TdebugPolicy>
210 if ( step[nSteps-1].getCorrectOccupation() )
212 step[nSteps-1].setCorrectOccupation();
213 correct_occupation_was_forced_flag =
true;
217 template<
typename Treal,
typename Tvector,
typename TdebugPolicy>
219 if (step[nSteps-1].getCorrectOccupation()) {
223 for (
int ind = 0; ind < nSteps; ind++) {
224 Treal XmX2Eucl = step[ind].getXmX2EuclNorm().upp();
225 if ( XmX2Eucl < 1 / (Treal)4) {
229 while (!middleInt.
empty() && i < nSteps - 1) {
230 middleInt.
puriStep(step[i].getPoly());
231 middleInt.
decrease(step[i+1].getActualThresh());
236 if (middleInt.
cover(0.5))
237 step[ind].setCorrectOccupation();
243 template<
typename Treal,
typename Tvector,
typename TdebugPolicy>
245 for (
int ind = nSteps - 2; ind >= 0; ind--) {
246 step[ind].exchangeInfoWithNext(step[ind + 1]);
248 for (
int ind = 0; ind < nSteps - 1; ind++) {
249 step[ind].exchangeInfoWithNext(step[ind + 1]);
253 template<
typename Treal,
typename Tvector,
typename TdebugPolicy>
256 for (
int ind = 0; ind < end; ind++) {
257 error += step[ind].subspaceError();
262 template<
typename Treal,
typename Tvector,
typename TdebugPolicy>
269 Treal tolError = tolSubspaceError - subspaceError(nSteps - 1);
270 Treal initialAltThresh = 1;
281 Treal lastThresh = 0;
282 if (nSteps == 0 || nSteps == 1) {
283 Treal lmax = eigFInterval.upp();
284 Treal lmin = eigFInterval.low();
287 (homoF.upp() - lmax) / (lmin - lmax));
291 step[nSteps - 2].getHomo().low());
292 initialAltThresh = getThreshIncreasingGap(initialGap);
293 lastThresh = step[nSteps - 2].getChosenThresh();
294 initialAltThresh = initialAltThresh > lastThresh / 10 ?
295 initialAltThresh : lastThresh / 10;
296 initialGap.
puriStep(step[nSteps - 2].getPoly());
298 if (initialGap.
empty())
302 if (1.0 - initialGap.
upp() < tolEigenvalError &&
303 initialGap.
low() < tolEigenvalError) {
309 Treal thetaPerStep = 0.0;
313 Treal currThresh = 0;
314 int stepsLeftPrev = -2;
316 while (stepsLeft > stepsLeftPrev) {
318 altThresh = initialAltThresh;
319 currThresh = thetaPerStep * gap.
length() / (1 + thetaPerStep);
320 currThresh = currThresh < altThresh ? currThresh : altThresh;
322 lastThresh = currThresh;
323 stepsLeftPrev = stepsLeft;
325 while (!gap.
empty() &&
326 (1.0 - gap.
upp() > tolEigenvalError ||
327 gap.
low() > tolEigenvalError)) {
328 altThresh = getThreshIncreasingGap(gap);
329 altThresh = altThresh > lastThresh / 10 ? altThresh : lastThresh / 10;
333 currThresh = thetaPerStep * gap.
length() / (1 + thetaPerStep);
334 currThresh = currThresh < altThresh ? currThresh : altThresh;
336 lastThresh = currThresh;
339 thetaPerStep = tolError / (stepsLeft + 1);
343 thetaPerStep = tolError / (stepsLeftPrev + 1);
344 thresh = thetaPerStep * initialGap.
length() / (1 + thetaPerStep);
350 template<
typename Treal,
typename Tvector,
typename TdebugPolicy>
354 estimateStepsLeft(stepsLeft, thresh);
356 step[nSteps - 1].setEstimatedStepsLeft(stepsLeft);
358 thresh = tolSubspaceError / 100;
362 step[nSteps - 2].getHomo().low());
363 if (!middleGap.
empty()) {
365 Treal altThresh = getThreshIncreasingGap(middleGap);
367 Treal lastThresh = step[nSteps - 2].getChosenThresh();
368 altThresh = altThresh > lastThresh / 10 ? altThresh : lastThresh / 10;
369 thresh = thresh < altThresh ? thresh : altThresh;
372 Treal tmp = 1 - xmX2Eucl.
upp() * 4.0;
375 thresh = thresh < altThresh ? thresh : altThresh;
385 template<
typename Treal,
typename Tvector,
typename TdebugPolicy>
389 Treal delta = middleGap.
upp() - x;
427 template<
typename Treal,
typename Tvector,
typename TdebugPolicy>
431 if (nSteps == 0 || nSteps == 1) {
435 Treal ep = 0.207106781186547;
446 homoCopy.
puriStep(step[nSteps - 2].getPoly());
447 lumoCopy.
puriStep(step[nSteps - 2].getPoly());
452 howAccurate = step[nSteps - 1].getChosenThresh() / 100;
453 Treal highestPossibleAccuracy = 10.0 * step[nSteps - 1].getEigAccLoss();
456 howAccurate = howAccurate > highestPossibleAccuracy ?
457 howAccurate : highestPossibleAccuracy;
459 if (homoCopy.
length() > 0.2 || lumoCopy.
length() > 0.2) {
461 if (step[nSteps - 2].getN0() / (n - nocc) > 0.5 &&
462 step[nSteps - 2].getN1() / nocc > 0.5 &&
463 step[nSteps - 2].getXmX2EuclNorm().midPoint() < ep)
470 bool computeHomo =
true;
471 bool computeLumo =
true;
472 if (homoIsAccuratelyKnown() ||
473 homoCopy.
upp() < 0.5 ||
476 if (lumoIsAccuratelyKnown() ||
477 lumoCopy.
low() > 0.5 ||
480 return computeHomo || computeLumo;
485 template<
typename Treal,
typename Tvector,
typename TdebugPolicy>
488 Treal lmax = eigFInterval.upp();
489 Treal lmin = eigFInterval.low();
492 homoF.intersect(altHomo);
493 lumoF.intersect(altLumo);
496 template<
typename Treal,
typename Tvector,
typename TdebugPolicy>
499 for (
int ind = 0; ind < nSteps; ++ind)
500 accTime += (
double)step[ind].getTimeSquare();
503 template<
typename Treal,
typename Tvector,
typename TdebugPolicy>
506 for (
int ind = 0; ind < nSteps; ++ind)
507 accTime += (
double)step[ind].getTimeThresh();
510 template<
typename Treal,
typename Tvector,
typename TdebugPolicy>
513 for (
int ind = 0; ind < nSteps; ++ind)
514 accTime += (
double)step[ind].getTimeXmX2Norm();
517 template<
typename Treal,
typename Tvector,
typename TdebugPolicy>
520 for (
int ind = 0; ind < nSteps; ++ind)
521 accTime += (
double)step[ind].getTimeTotal();
526 template<
typename Treal,
typename Tvector,
typename TdebugPolicy>
530 for (
int ind = 0; ind < nSteps; ++ind) {
532 step[ind].getTimeSquare()<<
" "<<
533 step[ind].getTimeThresh()<<
" "<<
534 step[ind].getTimeXmX2Norm()<<
" "<<
535 step[ind].getTimeTotal()<<
" "<<
536 step[ind].getTimeXX2Write()<<
" "<<
537 step[ind].getTimeXX2Read()<<
" "<<
540 s<<
"];\n"<<
"figure; bar(puriTime(:,1:3),'stacked')"<<std::endl<<
541 "legend('Matrix Square', 'Truncation', '||X-X^2||'),"<<
542 " xlabel('Iteration'), ylabel('Time (seconds)')"<<std::endl;
543 s<<
"figure; plot(puriTime(:,4),'-x')"<<std::endl;
546 template<
typename Treal,
typename Tvector,
typename TdebugPolicy>
549 s<<
"% PURIFICATION INFO IN MATLAB/OCTAVE FILE FORMAT"<<std::endl;
552 s<<
"n = "<<n<<
";"<<std::endl;
553 s<<
"nocc = "<<nocc<<
";"<<std::endl;
554 s<<
"tolSubspaceError = "<<tolSubspaceError<<
";"<<std::endl;
555 s<<
"tolEigenvalError = "<<tolEigenvalError<<
";"<<std::endl;
556 s<<
"correct_occupation_was_forced_flag = "
557 <<correct_occupation_was_forced_flag<<
";"<<std::endl;
558 s<<
"% Each row of the following matrix contains purification info \n"
559 <<
"% for one step.\n"
560 <<
"% The columns are arranged as: \n"
561 <<
"% ind, n0, n1, trace(X), trace(X*X), poly, actualThresh, delta, "
562 <<
"correctOcc, XmX2low, XmX2upp, HOMOmid, LUMOmid, "
563 <<
"nnz(X), nnz(X*X), chosenThresh, estimatedStepsLeft "
566 for (
int ind = 0; ind < nSteps; ++ind) {
569 step[ind].getN0()<<
" "<<
570 step[ind].getN1()<<
" "<<
571 step[ind].getTraceX()<<
" "<<
572 step[ind].getTraceX2()<<
" "<<
573 step[ind].getPoly()<<
" "<<
574 step[ind].getActualThresh()<<
" "<<
575 step[ind].getDelta()<<
" "<<
576 step[ind].getCorrectOccupation()<<
" "<<
577 step[ind].getXmX2EuclNorm().low()<<
" "<<
578 step[ind].getXmX2EuclNorm().upp()<<
" "<<
579 step[ind].getHomo().midPoint()<<
" "<<
580 step[ind].getLumo().midPoint()<<
" "<<
581 step[ind].getNnzX()<<
" "<<
582 step[ind].getNnzX2()<<
" "<<
583 step[ind].getChosenThresh()<<
" "<<
584 step[ind].getEstimatedStepsLeft()<<
" "<<
589 s<<
"ind = puriMat(:,1);\n";
591 <<
"plot(ind,puriMat(:,12),'x-b',ind,puriMat(:,13),'o-r')\n"
592 <<
"legend('HOMO','LUMO'),xlabel('Iteration')\n"
593 <<
"axis([0 max(ind) 0 1])\n";
595 <<
"plot(ind,100*puriMat(:,2)/(n-nocc),'o-r',...\n"
596 <<
"ind, 100*puriMat(:,3)/nocc,'x-b')\n"
597 <<
"legend('N0','N1'),xlabel('Iteration'), ylabel('Percentage')\n"
598 <<
"axis([0 max(ind) 0 100])\n";
601 <<
"plot(ind,100*puriMat(:,14)/(n*n),'o-r',...\n"
602 <<
"ind, 100*puriMat(:,15)/(n*n),'x-b')\n"
603 <<
"legend('nnz(X)','nnz(X^2)'),xlabel('Iteration') \n"
604 <<
"ylabel('Percentage')\n"
605 <<
"axis([0 max(ind) 0 100])\n";
607 <<
"semilogy(ind,puriMat(:,16),'x-r',ind,puriMat(:,7),'o-b')\n"
608 <<
"xlabel('Iteration'), ylabel('Threshold')\n"
609 <<
"legend('Chosen threshold', 'Actual threshold')\n"
610 <<
"axis([0 max(ind) min(puriMat(:,7))/10 max(puriMat(:,16))*10])\n";
612 <<
"plot(ind,ind(end:-1:1)-1,ind,puriMat(:,17),'x-r')\n"
613 <<
"legend('Steps left', 'Estimated steps left')\n"
614 <<
"xlabel('Iteration')\n";
619 template<
typename Treal,
typename Tvector,
typename TdebugPolicy>
622 s<<
"puriMemUsage = [";
623 for (
int ind = 0; ind < nSteps; ++ind) {
627 memUsageBeforeTrunc.
res <<
" "<<
628 memUsageBeforeTrunc.
virt<<
" "<<
629 memUsageBeforeTrunc.
peak<<
" "<<
630 memUsageInXmX2Diff.
res <<
" "<<
631 memUsageInXmX2Diff.
virt<<
" "<<
632 memUsageInXmX2Diff.
peak<<
" "<<
637 <<
"plot(puriMemUsage(:,1),'x-b')\n"
639 <<
"plot(puriMemUsage(:,2),'o-r')\n"
640 <<
"plot(puriMemUsage(:,3),'^-g')\n"
641 <<
"legend('Resident','Virtual','Peak'),xlabel('Iteration'),ylabel('Mem usage before trunc [GB]')\n"
642 <<
"%force y axis to start at 0\n"
643 <<
"axissaved = axis; axissaved(3) = 0; axis(axissaved);\n"
646 <<
"plot(puriMemUsage(:,4),'x-b')\n"
648 <<
"plot(puriMemUsage(:,5),'o-r')\n"
649 <<
"plot(puriMemUsage(:,6),'^-g')\n"
650 <<
"legend('Resident','Virtual','Peak'),xlabel('Iteration'),ylabel('Mem usage in XmX2Diff [GB]')\n"
651 <<
"%force y axis to start at 0\n"
652 <<
"axissaved = axis; axissaved(3) = 0; axis(axissaved);\n"
656 template<
typename Treal,
typename Tvector,
typename TdebugPolicy>
659 Tvector & eigVecHOMO)
const {
662 for (
int ind = 0; ind < nSteps; ++ind) {
663 if (!haveHOMO && step[ind].getHomoWasComputed()) {
664 eigVecHOMO = *step[ind].getEigVecPtr();
667 if (!haveLUMO && step[ind].getLumoWasComputed()) {
668 eigVecLUMO = *step[ind].getEigVecPtr();