PuriInfo.h

Go to the documentation of this file.
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

Generated on Wed Nov 21 09:32:11 2012 for ergo by  doxygen 1.4.7