00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00036 #ifndef MAT_PURIINFO
00037 #define MAT_PURIINFO
00038 #include <math.h>
00039 #include <iomanip>
00040 #include "PuriStepInfo.h"
00041 #define __ID__ "$Id: PuriInfo.h 4437 2012-07-05 09:01:18Z elias $"
00042 namespace mat {
00046 template<typename Treal, typename Tvector, typename TdebugPolicy>
00047 class PuriInfo : public TdebugPolicy {
00048 public:
00049 PuriInfo(int nn, int noc,
00050 Interval<Treal> eigFInt,
00051 Interval<Treal> hoF,
00052 Interval<Treal> luF,
00053 Treal toleratedEigenvalError,
00054 Treal toleratedSubspaceError,
00055 normType normForTruncation_,
00056 int maxS = 100)
00057 : n(nn),nocc(noc), step(new PuriStepInfo<Treal, Tvector, TdebugPolicy>[maxS]),
00058 maxSteps(maxS), nSteps(0), correct_occupation_was_forced_flag(false),
00059 eigFInterval(eigFInt), homoF(hoF), lumoF(luF),
00060 tolSubspaceError(toleratedSubspaceError),
00061 tolEigenvalError(toleratedEigenvalError),
00062 normForTruncation(normForTruncation_) {
00063 for (int ind = 0; ind < maxSteps; ++ind)
00064 step[ind] =
00065 PuriStepInfo<Treal, Tvector, TdebugPolicy>(n, nocc, tolEigenvalError);
00066 }
00067 virtual ~PuriInfo() {
00068 delete[] step;
00069 }
00071 void forceCorrectOccupation();
00076 void improveCorrectOccupation();
00080 void improveInfo();
00081 inline Interval<Treal> getEigFInterval() const {
00082 return eigFInterval;
00083 }
00084 inline PuriStepInfo<Treal, Tvector, TdebugPolicy> & getNext() {
00085 nSteps++;
00086 ASSERTALWAYS(nSteps - 1 < maxSteps);
00087 return step[nSteps - 1];
00088 }
00089 inline PuriStepInfo<Treal, Tvector, TdebugPolicy> & operator()(int const ind) {
00090 assert(ind >= 0);
00091 assert(ind < nSteps);
00092 return step[ind];
00093 }
00094
00095
00097 inline Treal subspaceError() const {
00098 return subspaceError(nSteps);
00099 }
00105 void estimateStepsLeft(int& stepsLeft, Treal& thresh) const;
00106 Treal getOptimalThresh() const;
00107 Treal getThreshIncreasingGap(Interval<Treal> const & middleGap) const;
00108 bool ShouldComputeXmX2EuclNormAccurately(Treal & howAccurate) const;
00112 inline Interval<Treal> getHomoF() const {return homoF;}
00116 inline Interval<Treal> getLumoF() const {return lumoF;}
00117 inline int getMaxSteps() const {return maxSteps;}
00118 inline int getNSteps() const {return nSteps;}
00119
00120
00121
00122 bool correct_occupation_was_forced() const
00123 {return correct_occupation_was_forced_flag;}
00124 void improveHomoLumoF();
00125
00126 inline bool converged() {return step[nSteps - 1].converged();}
00127
00128 double getAccumulatedTimeSquare() const;
00129 double getAccumulatedTimeThresh() const;
00130 double getAccumulatedTimeXmX2Norm() const;
00131 double getAccumulatedTimeTotal() const;
00132
00133 void mTimings(std::ostream & file) const;
00134 void mInfo(std::ostream & file) const;
00135 void mMemUsage(std::ostream & file) const;
00140 bool homoIsAccuratelyKnown() const {
00141 bool res = false;
00142 for(int ind = 0; ind < nSteps; ++ind)
00143 res = res || step[ind].homoIsAccuratelyKnown(tolSubspaceError / 100);
00144 return res;
00145 }
00150 bool lumoIsAccuratelyKnown() const {
00151 bool res = false;
00152 for(int ind = 0; ind < nSteps; ++ind)
00153 res = res || step[ind].lumoIsAccuratelyKnown(tolSubspaceError / 100);
00154 return res;
00155 }
00156
00157 normType getNormForTruncation() const {
00158 return normForTruncation;
00159 }
00160
00161 void getHOMOandLUMOeigVecs(Tvector & eigVecLUMO,
00162 Tvector & eigVecHOMO) const;
00163
00164 protected:
00165 int n;
00166 int nocc;
00168 PuriStepInfo<Treal, Tvector, TdebugPolicy>* step;
00169 int maxSteps;
00171 int nSteps;
00174 bool correct_occupation_was_forced_flag;
00176 Interval<Treal> const eigFInterval;
00181 Interval<Treal> homoF;
00186 Interval<Treal> lumoF;
00191 Treal const tolSubspaceError;
00195 Treal const tolEigenvalError;
00198 normType const normForTruncation;
00202 Treal subspaceError(int end) const;
00203
00204 private:
00205 };
00206
00207
00208 template<typename Treal, typename Tvector, typename TdebugPolicy>
00209 void PuriInfo<Treal, Tvector, TdebugPolicy>::forceCorrectOccupation() {
00210 if ( step[nSteps-1].getCorrectOccupation() )
00211 return;
00212 step[nSteps-1].setCorrectOccupation();
00213 correct_occupation_was_forced_flag = true;
00214 return;
00215 }
00216
00217 template<typename Treal, typename Tvector, typename TdebugPolicy>
00218 void PuriInfo<Treal, Tvector, TdebugPolicy>::improveCorrectOccupation() {
00219 if (step[nSteps-1].getCorrectOccupation()) {
00220 Treal distance;
00221 Interval<Treal> middleInt;
00222 Interval<Treal> zeroOneInt(0.0,1.0);
00223 for (int ind = 0; ind < nSteps; ind++) {
00224 Treal XmX2Eucl = step[ind].getXmX2EuclNorm().upp();
00225 if ( XmX2Eucl < 1 / (Treal)4) {
00226 distance = (1 - template_blas_sqrt(1 - 4 * XmX2Eucl)) / 2;
00227 middleInt = Interval<Treal>(distance, 1.0 - distance);
00228 int i = ind;
00229 while (!middleInt.empty() && i < nSteps - 1) {
00230 middleInt.puriStep(step[i].getPoly());
00231 middleInt.decrease(step[i+1].getActualThresh());
00232
00233 middleInt.intersect(zeroOneInt);
00234 ++i;
00235 }
00236 if (middleInt.cover(0.5))
00237 step[ind].setCorrectOccupation();
00238 }
00239 }
00240 }
00241 }
00242
00243 template<typename Treal, typename Tvector, typename TdebugPolicy>
00244 void PuriInfo<Treal, Tvector, TdebugPolicy>::improveInfo() {
00245 for (int ind = nSteps - 2; ind >= 0; ind--) {
00246 step[ind].exchangeInfoWithNext(step[ind + 1]);
00247 }
00248 for (int ind = 0; ind < nSteps - 1; ind++) {
00249 step[ind].exchangeInfoWithNext(step[ind + 1]);
00250 }
00251 }
00252
00253 template<typename Treal, typename Tvector, typename TdebugPolicy>
00254 Treal PuriInfo<Treal, Tvector, TdebugPolicy>::subspaceError(int end) const {
00255 Treal error = 0;
00256 for (int ind = 0; ind < end; ind++) {
00257 error += step[ind].subspaceError();
00258 }
00259 return error;
00260 }
00261
00262 template<typename Treal, typename Tvector, typename TdebugPolicy>
00263 void PuriInfo<Treal, Tvector, TdebugPolicy>::
00264 estimateStepsLeft(int& stepsLeft, Treal& thresh) const {
00265 stepsLeft = -1;
00266 thresh = 0;
00267 Interval<Treal> initialGap;
00268 Interval<Treal> gap;
00269 Treal tolError = tolSubspaceError - subspaceError(nSteps - 1);
00270 Treal initialAltThresh = 1;
00271 if (tolError <= 0.0)
00272 return;
00273
00274
00275
00276
00277
00278
00279
00280
00281 Treal lastThresh = 0;
00282 if (nSteps == 0 || nSteps == 1) {
00283 Treal lmax = eigFInterval.upp();
00284 Treal lmin = eigFInterval.low();
00285
00286 initialGap = Interval<Treal>((lumoF.low() - lmax) / (lmin - lmax),
00287 (homoF.upp() - lmax) / (lmin - lmax));
00288 }
00289 else {
00290 initialGap = Interval<Treal>(step[nSteps - 2].getLumo().upp(),
00291 step[nSteps - 2].getHomo().low());
00292 initialAltThresh = getThreshIncreasingGap(initialGap);
00293 lastThresh = step[nSteps - 2].getChosenThresh();
00294 initialAltThresh = initialAltThresh > lastThresh / 10 ?
00295 initialAltThresh : lastThresh / 10;
00296 initialGap.puriStep(step[nSteps - 2].getPoly());
00297 }
00298 if (initialGap.empty())
00299 return;
00300 #if 0
00301
00302 if (1.0 - initialGap.upp() < tolEigenvalError &&
00303 initialGap.low() < tolEigenvalError) {
00304 stepsLeft = 0;
00305 thresh = 0;
00306 return;
00307 }
00308 #endif
00309 Treal thetaPerStep = 0.0;
00310 Treal altThresh = 0;
00311
00312
00313 Treal currThresh = 0;
00314 int stepsLeftPrev = -2;
00315
00316 while (stepsLeft > stepsLeftPrev) {
00317 gap = initialGap;
00318 altThresh = initialAltThresh;
00319 currThresh = thetaPerStep * gap.length() / (1 + thetaPerStep);
00320 currThresh = currThresh < altThresh ? currThresh : altThresh;
00321 gap.decrease(currThresh);
00322 lastThresh = currThresh;
00323 stepsLeftPrev = stepsLeft;
00324 stepsLeft = 0;
00325 while (!gap.empty() &&
00326 (1.0 - gap.upp() > tolEigenvalError ||
00327 gap.low() > tolEigenvalError)) {
00328 altThresh = getThreshIncreasingGap(gap);
00329 altThresh = altThresh > lastThresh / 10 ? altThresh : lastThresh / 10;
00330
00331 gap.puriStep(gap.upp() + gap.low() < 1);
00332
00333 currThresh = thetaPerStep * gap.length() / (1 + thetaPerStep);
00334 currThresh = currThresh < altThresh ? currThresh : altThresh;
00335 gap.decrease(currThresh);
00336 lastThresh = currThresh;
00337 ++stepsLeft;
00338 }
00339 thetaPerStep = tolError / (stepsLeft + 1);
00340 }
00341
00342
00343 thetaPerStep = tolError / (stepsLeftPrev + 1);
00344 thresh = thetaPerStep * initialGap.length() / (1 + thetaPerStep);
00345 return;
00346 }
00347
00348
00349
00350 template<typename Treal, typename Tvector, typename TdebugPolicy>
00351 Treal PuriInfo<Treal, Tvector, TdebugPolicy>::getOptimalThresh() const {
00352 int stepsLeft = -1;
00353 Treal thresh = 0.0;
00354 estimateStepsLeft(stepsLeft, thresh);
00355 if (nSteps > 0)
00356 step[nSteps - 1].setEstimatedStepsLeft(stepsLeft);
00357 if (stepsLeft < 0)
00358 thresh = tolSubspaceError / 100;
00359 if (nSteps > 1) {
00360 Interval<Treal> middleGap =
00361 Interval<Treal>(step[nSteps - 2].getLumo().upp(),
00362 step[nSteps - 2].getHomo().low());
00363 if (!middleGap.empty()) {
00364
00365 Treal altThresh = getThreshIncreasingGap(middleGap);
00366
00367 Treal lastThresh = step[nSteps - 2].getChosenThresh();
00368 altThresh = altThresh > lastThresh / 10 ? altThresh : lastThresh / 10;
00369 thresh = thresh < altThresh ? thresh : altThresh;
00370 #if 0
00371 Interval<Treal> xmX2Eucl = step[nSteps - 2].getXmX2EuclNorm();
00372 Treal tmp = 1 - xmX2Eucl.upp() * 4.0;
00373 if (tmp > 0) {
00374 Treal altThresh = 0.25 * (1 - template_blas_sqrt(tmp)) / 2;
00375 thresh = thresh < altThresh ? thresh : altThresh;
00376 }
00377 #endif
00378 }
00379 }
00380
00381 ASSERTALWAYS(thresh >= 0.0);
00382 return thresh;
00383 }
00384
00385 template<typename Treal, typename Tvector, typename TdebugPolicy>
00386 Treal PuriInfo<Treal, Tvector, TdebugPolicy>::
00387 getThreshIncreasingGap(Interval<Treal> const & middleGap) const {
00388 Treal x = middleGap.midPoint();
00389 Treal delta = middleGap.upp() - x;
00390 Treal thresh;
00391 #if 0
00392 thresh = delta * template_blas_fabs(2 * x - 1) / 10;
00393 #else
00394
00395
00396
00397
00398
00399
00400
00401 if (delta > 0.4)
00402
00403
00404
00405
00406
00407
00408
00409 thresh = delta * template_blas_fabs(2 * x - 1) * template_blas_fabs(2 * x - 1);
00410 else
00411 thresh = delta * template_blas_fabs(2 * x - 1) / 2;
00412
00413 #endif
00414
00415
00416
00417
00418
00419
00420
00421
00422 ASSERTALWAYS(thresh <= delta * template_blas_fabs(2 * x - 1) + std::numeric_limits<Treal>::epsilon());
00423 ASSERTALWAYS(thresh >= 0.0);
00424 return thresh;
00425 }
00426
00427 template<typename Treal, typename Tvector, typename TdebugPolicy>
00428 bool PuriInfo<Treal, Tvector, TdebugPolicy>::
00429 ShouldComputeXmX2EuclNormAccurately(Treal & howAccurate) const {
00430
00431 if (nSteps == 0 || nSteps == 1) {
00432 howAccurate = 0;
00433 return false;
00434 }
00435 Treal ep = 0.207106781186547;
00436
00437
00438
00439
00440
00441
00442
00443
00444 Interval<Treal> homoCopy = step[nSteps - 2].getHomo();
00445 Interval<Treal> lumoCopy = step[nSteps - 2].getLumo();
00446 homoCopy.puriStep(step[nSteps - 2].getPoly());
00447 lumoCopy.puriStep(step[nSteps - 2].getPoly());
00448
00449
00450
00451
00452 howAccurate = step[nSteps - 1].getChosenThresh() / 100;
00453 Treal highestPossibleAccuracy = 10.0 * step[nSteps - 1].getEigAccLoss();
00454 ASSERTALWAYS(howAccurate >= 0);
00455 ASSERTALWAYS(highestPossibleAccuracy >= 0);
00456 howAccurate = howAccurate > highestPossibleAccuracy ?
00457 howAccurate : highestPossibleAccuracy;
00458
00459 if (homoCopy.length() > 0.2 || lumoCopy.length() > 0.2) {
00460
00461 if (step[nSteps - 2].getN0() / (n - nocc) > 0.5 &&
00462 step[nSteps - 2].getN1() / nocc > 0.5 &&
00463 step[nSteps - 2].getXmX2EuclNorm().midPoint() < ep)
00464 return true;
00465 else
00466 return false;
00467 }
00468 else {
00469
00470 bool computeHomo = true;
00471 bool computeLumo = true;
00472 if (homoIsAccuratelyKnown() ||
00473 homoCopy.upp() < 0.5 ||
00474 template_blas_fabs(lumoCopy.low() - 0.5) < template_blas_fabs(homoCopy.low() - 0.5))
00475 computeHomo = false;
00476 if (lumoIsAccuratelyKnown() ||
00477 lumoCopy.low() > 0.5 ||
00478 template_blas_fabs(homoCopy.upp() - 0.5) < template_blas_fabs(lumoCopy.upp() - 0.5))
00479 computeLumo = false;
00480 return computeHomo || computeLumo;
00481 }
00482 }
00483
00484
00485 template<typename Treal, typename Tvector, typename TdebugPolicy>
00486 void PuriInfo<Treal, Tvector, TdebugPolicy>::
00487 improveHomoLumoF() {
00488 Treal lmax = eigFInterval.upp();
00489 Treal lmin = eigFInterval.low();
00490 Interval<Treal> altHomo(step[0].getHomo() * (lmin - lmax) + lmax);
00491 Interval<Treal> altLumo(step[0].getLumo() * (lmin - lmax) + lmax);
00492 homoF.intersect(altHomo);
00493 lumoF.intersect(altLumo);
00494 }
00495
00496 template<typename Treal, typename Tvector, typename TdebugPolicy>
00497 double PuriInfo<Treal, Tvector, TdebugPolicy>::getAccumulatedTimeSquare() const {
00498 double accTime = 0;
00499 for (int ind = 0; ind < nSteps; ++ind)
00500 accTime += (double)step[ind].getTimeSquare();
00501 return accTime;
00502 }
00503 template<typename Treal, typename Tvector, typename TdebugPolicy>
00504 double PuriInfo<Treal, Tvector, TdebugPolicy>::getAccumulatedTimeThresh() const {
00505 double accTime = 0;
00506 for (int ind = 0; ind < nSteps; ++ind)
00507 accTime += (double)step[ind].getTimeThresh();
00508 return accTime;
00509 }
00510 template<typename Treal, typename Tvector, typename TdebugPolicy>
00511 double PuriInfo<Treal, Tvector, TdebugPolicy>::getAccumulatedTimeXmX2Norm() const {
00512 double accTime = 0;
00513 for (int ind = 0; ind < nSteps; ++ind)
00514 accTime += (double)step[ind].getTimeXmX2Norm();
00515 return accTime;
00516 }
00517 template<typename Treal, typename Tvector, typename TdebugPolicy>
00518 double PuriInfo<Treal, Tvector, TdebugPolicy>::getAccumulatedTimeTotal() const {
00519 double accTime = 0;
00520 for (int ind = 0; ind < nSteps; ++ind)
00521 accTime += (double)step[ind].getTimeTotal();
00522 return accTime;
00523 }
00524
00525
00526 template<typename Treal, typename Tvector, typename TdebugPolicy>
00527 void PuriInfo<Treal, Tvector, TdebugPolicy>::
00528 mTimings(std::ostream & s) const {
00529 s<<"puriTime = [";
00530 for (int ind = 0; ind < nSteps; ++ind) {
00531 s<<
00532 step[ind].getTimeSquare()<<" "<<
00533 step[ind].getTimeThresh()<<" "<<
00534 step[ind].getTimeXmX2Norm()<<" "<<
00535 step[ind].getTimeTotal()<<" "<<
00536 step[ind].getTimeXX2Write()<<" "<<
00537 step[ind].getTimeXX2Read()<<" "<<
00538 std::endl;
00539 }
00540 s<<"];\n"<<"figure; bar(puriTime(:,1:3),'stacked')"<<std::endl<<
00541 "legend('Matrix Square', 'Truncation', '||X-X^2||'),"<<
00542 " xlabel('Iteration'), ylabel('Time (seconds)')"<<std::endl;
00543 s<<"figure; plot(puriTime(:,4),'-x')"<<std::endl;
00544 }
00545
00546 template<typename Treal, typename Tvector, typename TdebugPolicy>
00547 void PuriInfo<Treal, Tvector, TdebugPolicy>::
00548 mInfo(std::ostream & s) const {
00549 s<<"% PURIFICATION INFO IN MATLAB/OCTAVE FILE FORMAT"<<std::endl;
00550 s<<"% Norm for truncation: "<< getNormTypeString(normForTruncation)
00551 << std::endl;
00552 s<<"n = "<<n<<";"<<std::endl;
00553 s<<"nocc = "<<nocc<<";"<<std::endl;
00554 s<<"tolSubspaceError = "<<tolSubspaceError<<";"<<std::endl;
00555 s<<"tolEigenvalError = "<<tolEigenvalError<<";"<<std::endl;
00556 s<<"correct_occupation_was_forced_flag = "
00557 <<correct_occupation_was_forced_flag<<";"<<std::endl;
00558 s<<"% Each row of the following matrix contains purification info \n"
00559 <<"% for one step.\n"
00560 <<"% The columns are arranged as: \n"
00561 <<"% ind, n0, n1, trace(X), trace(X*X), poly, actualThresh, delta, "
00562 <<"correctOcc, XmX2low, XmX2upp, HOMOmid, LUMOmid, "
00563 <<"nnz(X), nnz(X*X), chosenThresh, estimatedStepsLeft "
00564 <<std::endl;
00565 s<<"puriMat = [";
00566 for (int ind = 0; ind < nSteps; ++ind) {
00567 s<<
00568 ind+1<<" "<<
00569 step[ind].getN0()<<" "<<
00570 step[ind].getN1()<<" "<<
00571 step[ind].getTraceX()<<" "<<
00572 step[ind].getTraceX2()<<" "<<
00573 step[ind].getPoly()<<" "<<
00574 step[ind].getActualThresh()<<" "<<
00575 step[ind].getDelta()<<" "<<
00576 step[ind].getCorrectOccupation()<<" "<<
00577 step[ind].getXmX2EuclNorm().low()<<" "<<
00578 step[ind].getXmX2EuclNorm().upp()<<" "<<
00579 step[ind].getHomo().midPoint()<<" "<<
00580 step[ind].getLumo().midPoint()<<" "<<
00581 step[ind].getNnzX()<<" "<<
00582 step[ind].getNnzX2()<<" "<<
00583 step[ind].getChosenThresh()<<" "<<
00584 step[ind].getEstimatedStepsLeft()<<" "<<
00585 std::endl;
00586 }
00587 s<<"];\n";
00588 s<<"if(1)\n";
00589 s<<"ind = puriMat(:,1);\n";
00590 s<<"figure; \n"
00591 <<"plot(ind,puriMat(:,12),'x-b',ind,puriMat(:,13),'o-r')\n"
00592 <<"legend('HOMO','LUMO'),xlabel('Iteration')\n"
00593 <<"axis([0 max(ind) 0 1])\n";
00594 s<<"figure; \n"
00595 <<"plot(ind,100*puriMat(:,2)/(n-nocc),'o-r',...\n"
00596 <<"ind, 100*puriMat(:,3)/nocc,'x-b')\n"
00597 <<"legend('N0','N1'),xlabel('Iteration'), ylabel('Percentage')\n"
00598 <<"axis([0 max(ind) 0 100])\n";
00599 s<<"figure; \n"
00600 <<"subplot(211)\n"
00601 <<"plot(ind,100*puriMat(:,14)/(n*n),'o-r',...\n"
00602 <<"ind, 100*puriMat(:,15)/(n*n),'x-b')\n"
00603 <<"legend('nnz(X)','nnz(X^2)'),xlabel('Iteration') \n"
00604 <<"ylabel('Percentage')\n"
00605 <<"axis([0 max(ind) 0 100])\n";
00606 s<<"subplot(212)\n"
00607 <<"semilogy(ind,puriMat(:,16),'x-r',ind,puriMat(:,7),'o-b')\n"
00608 <<"xlabel('Iteration'), ylabel('Threshold')\n"
00609 <<"legend('Chosen threshold', 'Actual threshold')\n"
00610 <<"axis([0 max(ind) min(puriMat(:,7))/10 max(puriMat(:,16))*10])\n";
00611 s<<"figure; \n"
00612 <<"plot(ind,ind(end:-1:1)-1,ind,puriMat(:,17),'x-r')\n"
00613 <<"legend('Steps left', 'Estimated steps left')\n"
00614 <<"xlabel('Iteration')\n";
00615
00616 s<<"end\n";
00617 }
00618
00619 template<typename Treal, typename Tvector, typename TdebugPolicy>
00620 void PuriInfo<Treal, Tvector, TdebugPolicy>::
00621 mMemUsage(std::ostream & s) const {
00622 s<<"puriMemUsage = [";
00623 for (int ind = 0; ind < nSteps; ++ind) {
00624 MemUsage::Values memUsageBeforeTrunc = step[ind].getMemUsageBeforeTrunc();
00625 MemUsage::Values memUsageInXmX2Diff = step[ind].getMemUsageInXmX2Diff();
00626 s<<
00627 memUsageBeforeTrunc.res <<" "<<
00628 memUsageBeforeTrunc.virt<<" "<<
00629 memUsageBeforeTrunc.peak<<" "<<
00630 memUsageInXmX2Diff.res <<" "<<
00631 memUsageInXmX2Diff.virt<<" "<<
00632 memUsageInXmX2Diff.peak<<" "<<
00633 std::endl;
00634 }
00635 s<<"];\n";
00636 s<<"figure; \n"
00637 <<"plot(puriMemUsage(:,1),'x-b')\n"
00638 << "hold on\n"
00639 <<"plot(puriMemUsage(:,2),'o-r')\n"
00640 <<"plot(puriMemUsage(:,3),'^-g')\n"
00641 <<"legend('Resident','Virtual','Peak'),xlabel('Iteration'),ylabel('Mem usage before trunc [GB]')\n"
00642 <<"%force y axis to start at 0\n"
00643 <<"axissaved = axis; axissaved(3) = 0; axis(axissaved);\n"
00644 << std::endl;
00645 s<<"figure; \n"
00646 <<"plot(puriMemUsage(:,4),'x-b')\n"
00647 << "hold on\n"
00648 <<"plot(puriMemUsage(:,5),'o-r')\n"
00649 <<"plot(puriMemUsage(:,6),'^-g')\n"
00650 <<"legend('Resident','Virtual','Peak'),xlabel('Iteration'),ylabel('Mem usage in XmX2Diff [GB]')\n"
00651 <<"%force y axis to start at 0\n"
00652 <<"axissaved = axis; axissaved(3) = 0; axis(axissaved);\n"
00653 << std::endl;
00654 }
00655
00656 template<typename Treal, typename Tvector, typename TdebugPolicy>
00657 void PuriInfo<Treal, Tvector, TdebugPolicy>::
00658 getHOMOandLUMOeigVecs(Tvector & eigVecLUMO,
00659 Tvector & eigVecHOMO) const {
00660 bool haveHOMO = 0;
00661 bool haveLUMO = 0;
00662 for (int ind = 0; ind < nSteps; ++ind) {
00663 if (!haveHOMO && step[ind].getHomoWasComputed()) {
00664 eigVecHOMO = *step[ind].getEigVecPtr();
00665 haveHOMO = 1;
00666 }
00667 if (!haveLUMO && step[ind].getLumoWasComputed()) {
00668 eigVecLUMO = *step[ind].getEigVecPtr();
00669 haveLUMO = 1;
00670 }
00671 }
00672 }
00673
00674
00675 }
00676 #undef __ID__
00677 #endif