ergo
|
00001 /* Ergo, version 3.2, a program for linear scaling electronic structure 00002 * calculations. 00003 * Copyright (C) 2012 Elias Rudberg, Emanuel H. Rubensson, and Pawel Salek. 00004 * 00005 * This program is free software: you can redistribute it and/or modify 00006 * it under the terms of the GNU General Public License as published by 00007 * the Free Software Foundation, either version 3 of the License, or 00008 * (at your option) any later version. 00009 * 00010 * This program is distributed in the hope that it will be useful, 00011 * but WITHOUT ANY WARRANTY; without even the implied warranty of 00012 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 00013 * GNU General Public License for more details. 00014 * 00015 * You should have received a copy of the GNU General Public License 00016 * along with this program. If not, see <http://www.gnu.org/licenses/>. 00017 * 00018 * Primary academic reference: 00019 * KohnâSham Density Functional Theory Electronic Structure Calculations 00020 * with Linearly Scaling Computational Time and Memory Usage, 00021 * Elias Rudberg, Emanuel H. Rubensson, and Pawel Salek, 00022 * J. Chem. Theory Comput. 7, 340 (2011), 00023 * <http://dx.doi.org/10.1021/ct100611z> 00024 * 00025 * For further information about Ergo, see <http://www.ergoscf.org>. 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 /* Tries to improve the homoF and lumoF eigenvalues. 00120 * Called after the purification process. 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 /* Make sure we stay in [0, 1] */ 00233 middleInt.intersect(zeroOneInt); 00234 ++i; 00235 } /* end while */ 00236 if (middleInt.cover(0.5)) 00237 step[ind].setCorrectOccupation(); 00238 } /* end if */ 00239 } /* end for */ 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 /* Compute number of steps needed to converge. */ 00275 /* Compute initial interval */ 00276 /* nSteps == 0 means that a Purification object has not been 00277 * associated to this instance yet. 00278 * nSteps == 1 means that a Purification object has been associated 00279 * to this instance but no steps have been taken yet. 00280 */ 00281 Treal lastThresh = 0; 00282 if (nSteps == 0 || nSteps == 1) { 00283 Treal lmax = eigFInterval.upp(); 00284 Treal lmin = eigFInterval.low(); 00285 /* Compute transformed homo and lumo eigenvalues. */ 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 /* Already converged? */ 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; /* Tolerated subspace error per step. */ 00310 Treal altThresh = 0; /* Maximum threshold that guarantees increased 00311 * gap. 00312 */ 00313 Treal currThresh = 0; /* Chosen threshold. */ 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 /* Compute optimal threshold for convergence of subspace. */ 00343 thetaPerStep = tolError / (stepsLeftPrev + 1); 00344 thresh = thetaPerStep * initialGap.length() / (1 + thetaPerStep); 00345 return; 00346 } 00347 00348 00349 /* FIXME !! */ 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; /* No reason */ 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 /* Get thresh that definitely gives increasing gap. */ 00365 Treal altThresh = getThreshIncreasingGap(middleGap); 00366 /* Allow thresh to decrease only to a ten times smaller threshold */ 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 //std::cout<<"nsteps left, thresh: "<<stepsLeft<<" , "<<thresh<<std::endl; 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 /* After a purification step, truncation is done. 00395 * The threshold t (measured by the Euclidean norm of the error matrix) 00396 * has to satisfy t <= | 1 - 2x | * d / 2 where x and d are the midpoint 00397 * and the length of the HOMO-LUMO gap respectively. In this way the 00398 * gap is guaranteed to increase to the next step. Note that x = 0.5 and 00399 * d = 0 gives zero threshold. x -> 0.5 as the process converges. 00400 */ 00401 if (delta > 0.4) 00402 /* This choice gives much better estimation of the remaining 00403 * number of iterations. 00404 */ 00405 /* Close to convergence we chose quadratical dependence on the 00406 * midpoint since we expect the midpoint to go to 0.5 quadratically 00407 * at convergence. 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 // thresh = thresh > 1e-7 ? thresh : 1e-7; 00413 #endif 00414 /************************************************************** 00415 ELIAS NOTE 2010-11-18: I got assertion failure when using gcc 00416 in Fedora 14 [g++ (GCC) 4.5.1 20100924 (Red Hat 4.5.1-4)] and 00417 it turned out to be because the fabs call returned zero while 00418 thresh was something like 1e-30 for single precision. Therefore 00419 I added std::numeric_limits<Treal>::epsilon() in the assertion, 00420 which seemed to solve the problem. 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 /* = (sqrt(2) - 1) / 2 00437 * This value is obtained by noting that x = 1 / sqrt(2) -> x^2 = 1 / 2. 00438 * If an interval ]1 - 1 / sqrt(2), 1 / sqrt(2)[ is empty from eigenvalues, 00439 * and if the occ. count is correct, then the occupation 00440 * count is guaranteed to be correct after one step as well. 00441 * This transforms to the value (sqrt(2) - 1) / 2 for the 00442 * ||X - X^2||_2 norm via theorem 1. 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 /* Note: we changed this from getActualThresh() to 00449 getChosenThresh() to avoid ridiculously small values for 00450 howAccurate, as happened earlier for cases when no truncation 00451 occurred, i.e. when the matrices were small. */ 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 /* Base decision on n0 and n1 and XmX2Norm from previous step */ 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 /* Decision can probably be made from homo and lumo estimates */ 00470 bool computeHomo = true; /* Do we want to try to compute homo */ 00471 bool computeLumo = true; /* Do we want to try to compute lumo */ 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 } /* end namespace mat */ 00676 #undef __ID__ 00677 #endif