PuriStepInfo.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_PURISTEPINFO
00037 #define MAT_PURISTEPINFO
00038 #include <math.h>
00039 #include <iomanip>
00040 #include "PuriStepInfoDebug.h"
00041 #define __ID__ "$Id: PuriStepInfo.h 4437 2012-07-05 09:01:18Z elias $"
00042 namespace mat {
00043   /* IDEA: Adjust threshold (looser) in later iteration if threshold in early
00044    * iteration was tighter than expected.
00045    */
00046 
00054   template<typename Treal, typename Tvector, typename TdebugPolicy>
00055     class PuriStepInfo : public PuriStepInfoDebug<Treal, TdebugPolicy> {
00056   public:
00057     explicit PuriStepInfo(int nn = -1, int noc = -1, Treal eigvalConvCrit = 0.0)
00058       : n(nn),nocc(noc), traceX(-1.0), traceX2(-1.0), poly(-1), 
00059       chosenThresh(0.0), actualThresh(0.0), 
00060       estimatedStepsLeft(-1),
00061       eigInterval(-1.0,2.0), delta(-1.0), correctOccupation(0),
00062       XmX2EuclNorm(0.0, 10.0), eigVecPtr(0), 
00063       lumoWasComputed(0), homoWasComputed(0),
00064       n0(0.0), n1(0.0), 
00065       homo(0.0, 1.0), lumo(0.0, 1.0), eigConvCrit(eigvalConvCrit), 
00066       nnzX(0), nnzX2(0), eigAccLoss(0),
00067       timeThresh(0),  timeSquare(0),  timeXmX2Norm(0), timeTotal(0),
00068       timeXX2Write(0), timeXX2Read(0) {}
00069     
00070     ~PuriStepInfo() {
00071       delete eigVecPtr;
00072     } 
00073       bool converged() const {
00074         bool homoLumo_converged = (1 - homo.low() < eigConvCrit && 
00075                                    lumo.upp() < eigConvCrit);
00076         bool XmX2norm_converged = getXmX2EuclNorm().upp() < eigConvCrit;
00077         // FIXME: Possibly, propagating XmX2norm between the
00078         // iterations can give more accurate values since this norm is not
00079         // computed accurately in each iteration.
00080         return homoLumo_converged || XmX2norm_converged;
00081       }
00082       
00083       inline void setChosenThresh(Treal const thr) {chosenThresh = thr;}
00084       inline Treal getChosenThresh() const { return chosenThresh;}
00085       inline void setActualThresh(Treal const thr) {actualThresh = thr;}
00086       inline Treal getActualThresh() const { return actualThresh;}
00087       inline void setEstimatedStepsLeft(int const stepsleft) {
00088         estimatedStepsLeft = stepsleft;
00089       }
00090       inline int getEstimatedStepsLeft() const { return estimatedStepsLeft;}
00091       inline void setTraceX(Treal const trX) {traceX = trX;}
00092       inline void setTraceX2(Treal const trX2) {traceX2 = trX2;}
00093       inline Treal getTraceX() const {return traceX;}
00094       inline Treal getTraceX2() const {return traceX2;}
00095       void setPoly();
00096       inline int getPoly() const {return poly;}
00098       inline void setXmX2EuclNorm(Interval<Treal> const XmX2eucl) {
00099         XmX2EuclNorm = XmX2eucl;
00100       }
00101       
00102       inline Interval<Treal> getXmX2EuclNorm() const {return XmX2EuclNorm;}
00107       void setEigVecPtr(Tvector * eigVecPtr_) {
00108         delete eigVecPtr;
00109         eigVecPtr = eigVecPtr_;
00110       }
00111       Tvector const * const getEigVecPtr() const {return eigVecPtr;}
00112 
00113       void improveHomoLumo(Interval<Treal> const homoInt, 
00114                            Interval<Treal> const lumoInt);
00115       inline Interval<Treal> const & getHomo() const {
00116         return homo;
00117       }
00118       inline Interval<Treal> const & getLumo() const {
00119         return lumo;
00120       }
00121       inline Interval<Treal> const & getEigInterval() const {
00122         return eigInterval;
00123       }
00124       void exchangeInfoWithNext(PuriStepInfo<Treal, Tvector, TdebugPolicy> & next);
00127       inline void setCorrectOccupation() {correctOccupation = 1;}
00128       inline int getCorrectOccupation() const {return correctOccupation;}
00129       Treal subspaceError() const;
00133       void improveEigInterval(Interval<Treal> const eInt);
00134       
00135       inline void setNnzX(size_t const nzX) {nnzX = nzX;}
00136       inline size_t getNnzX() const {return nnzX;}
00137       inline void setNnzX2(size_t const nzX2) {nnzX2 = nzX2;}
00138       inline size_t getNnzX2() const {return nnzX2;}
00139       void computeEigAccLoss();
00140       inline Treal getEigAccLoss() const {return eigAccLoss;}
00141       
00142       inline Treal getN0() const {return n0;}
00143       inline Treal getN1() const {return n1;}
00144       inline Treal getDelta() const {return delta;}
00145       
00146 
00147       inline void checkIntervals(const char* descriptionString) const {
00148         PuriStepInfoDebug<Treal, TdebugPolicy>::
00149 	  checkIntervals(eigInterval, homo, lumo, XmX2EuclNorm, descriptionString);    
00150       }
00151       template<typename Tmatrix>
00152         inline void computeExactValues(Tmatrix const & X, Tmatrix const & X2) {
00153         PuriStepInfoDebug<Treal, TdebugPolicy>::
00154 	  computeExactValues(X, X2, n, nocc);
00155       }
00156   
00157       /* Functions to set and get timings and mem usage info for the different steps: */
00158       void setMemUsageBeforeTrunc() {
00159         MemUsage::getMemUsage(memUsageBeforeTrunc);
00160       }
00161       void setMemUsageInXmX2Diff(MemUsage::Values & memUsage) {
00162         memUsageInXmX2Diff = memUsage;
00163       }
00164       void setTimeThresh(float t) {timeThresh = t;}
00165       void setTimeSquare(float t) {timeSquare = t;}
00166       void setTimeXmX2Norm(float t) {timeXmX2Norm = t;}
00167       void setTimeTotal(float t) {timeTotal = t;}
00168       void setTimeXX2Write(float t) {timeXX2Write = t;}
00169       void setTimeXX2Read(float t) {timeXX2Read = t;}
00170 
00171       MemUsage::Values getMemUsageBeforeTrunc() { return memUsageBeforeTrunc; }
00172       MemUsage::Values getMemUsageInXmX2Diff()  { return memUsageInXmX2Diff;  }
00173       float getTimeThresh() {return timeThresh;}
00174       float getTimeSquare() {return timeSquare;}
00175       float getTimeXmX2Norm() {return timeXmX2Norm;}
00176       float getTimeTotal() {return timeTotal;}
00177       float getTimeXX2Write() {return timeXX2Write;}
00178       float getTimeXX2Read() {return timeXX2Read;}
00179       
00180       inline bool homoIsAccuratelyKnown
00181         (Treal accuracyLimit 
00184          ) const {
00185         return homo.length() < accuracyLimit;
00186       }
00187       inline bool lumoIsAccuratelyKnown
00188         (Treal accuracyLimit 
00191          ) const {
00192         return lumo.length() < accuracyLimit;
00193       }
00194 
00195       inline bool getLumoWasComputed() {return lumoWasComputed;}
00196       inline bool getHomoWasComputed() {return homoWasComputed;}
00197       
00198   protected:
00202       void computen0n1();
00203       
00204       int n; 
00205       int nocc; 
00206       /* FIXME?: Upper and lower bound on traces? */
00207       Treal traceX; 
00208       Treal traceX2; 
00209       int poly; 
00210       Treal chosenThresh; 
00211       Treal actualThresh; 
00212       int estimatedStepsLeft; 
00216       Interval<Treal> eigInterval; 
00218       Treal delta; 
00219       int correctOccupation; 
00226       Interval<Treal> XmX2EuclNorm; 
00230       Tvector * eigVecPtr; 
00233       bool lumoWasComputed; 
00234       bool homoWasComputed; 
00235       Treal n0; 
00238       Treal n1; 
00242       Interval<Treal> homo; 
00243       Interval<Treal> lumo; 
00244       Treal eigConvCrit; 
00247       size_t nnzX;  
00248       size_t nnzX2; 
00254       Treal eigAccLoss;
00255       
00256       /* Variables for the time and mem usage of various operations */
00257       MemUsage::Values memUsageBeforeTrunc;
00258       MemUsage::Values memUsageInXmX2Diff;
00259       float timeThresh;
00260       float timeSquare;
00261       float timeXmX2Norm;
00262       float timeTotal;
00263       float timeXX2Write;
00264       float timeXX2Read;
00265   private:
00266   };
00267 
00268 #if 1
00269   template<typename Treal, typename Tvector, typename TdebugPolicy>
00270     std::ostream& operator<<(std::ostream& s, 
00271                              PuriStepInfo<Treal, Tvector, TdebugPolicy> const & psi) {
00272     s<<" Trace X: "<<psi.getTraceX()
00273      <<" Trace X2: "<<psi.getTraceX2()
00274      <<" poly: "<<psi.getPoly()
00275      <<" ||E||_2: "<<psi.getActualThresh()
00276      <<" Eiginterval: "<<psi.getEigInterval()
00277      <<" correctOccupation: "<<psi.getCorrectOccupation()
00278      <<std::endl
00279      <<" XmX2EuclNorm: "<<psi.getXmX2EuclNorm()
00280      <<" n0: "<<psi.getN0()
00281      <<" n1: "<<psi.getN1()
00282      <<" homo: "<<psi.getHomo()
00283      <<" lumo: "<<psi.getLumo()
00284      <<" nnzX: "<<psi.getNnzX()
00285      <<" nnzX2: "<<psi.getNnzX2();
00286     return s;
00287   }
00288 
00289 #endif
00290 
00291   template<typename Treal, typename Tvector, typename TdebugPolicy>
00292     void PuriStepInfo<Treal, Tvector, TdebugPolicy>::setPoly() {
00293     if (Interval<Treal>::intersect(homo,lumo).empty()) {
00294       ASSERTALWAYS(homo.low() > lumo.upp());
00295       if (homo.low() + lumo.upp() > 1)
00296         poly = 0; /* x*x */
00297       else
00298         poly = 1; /* 2*x - x*x */
00299     }
00300     else {
00301       if (template_blas_fabs(traceX2 - nocc) < template_blas_fabs(2 * traceX - traceX2 - nocc))
00302         poly = 0; /* x*x */
00303       else
00304         poly = 1; /* 2*x - x*x */
00305     }
00306   }
00307   
00308   template<typename Treal, typename Tvector, typename TdebugPolicy>
00309     void PuriStepInfo<Treal, Tvector, TdebugPolicy>::
00310     improveHomoLumo(Interval<Treal> const homoInt, 
00311                     Interval<Treal> const lumoInt) {
00312     checkIntervals("PuriStepInfo::improveHomoLumo A0");
00313     homo.intersect(homoInt);
00314     checkIntervals("PuriStepInfo::improveHomoLumo A1");
00315     lumo.intersect(lumoInt);
00316     checkIntervals("PuriStepInfo::improveHomoLumo A2");
00317     ASSERTALWAYS(!homo.empty());
00318     ASSERTALWAYS(!lumo.empty());
00319     if (homo.low() > 0.5 && lumo.upp() < 0.5)
00320       this->setCorrectOccupation();
00321     
00322     if (correctOccupation && 1.0 - XmX2EuclNorm.upp() * 4.0 > 0) {
00323       Interval<Treal> tmp = 
00324         sqrtInt((XmX2EuclNorm * (Treal)(-4.0)) + (Treal)1.0);
00325       ASSERTALWAYS(tmp.length() > 0);
00326 
00327       ASSERTALWAYS(!homo.empty());
00328       ASSERTALWAYS(!lumo.empty());
00329 
00330       homo.intersect(Interval<Treal>((1.0 + tmp.low()) / 2.0, homo.upp()));
00331 
00332       ASSERTALWAYS(!homo.empty());
00333       ASSERTALWAYS(!lumo.empty());
00334 
00335       lumo.intersect(Interval<Treal>(lumo.low(), (1.0 - tmp.low()) / 2.0));
00336       checkIntervals("PuriStepInfo::improveHomoLumo B");
00337       ASSERTALWAYS(!homo.empty());
00338       ASSERTALWAYS(!lumo.empty());
00339       if (!eigInterval.cover((1 + template_blas_sqrt(1.0 + 4.0 * XmX2EuclNorm.low())) / 2) &&
00340           !eigInterval.cover((1 - template_blas_sqrt(1.0 + 4.0 * XmX2EuclNorm.low())) / 2)){
00341         /* Either homo lies in homoTmp or lumo lies in lumoTmp. */
00342         Interval<Treal> homoTmp = 
00343           (tmp + (Treal)1.0) / (Treal)2.0;
00344         Interval<Treal> lumoTmp = 
00345           ((tmp * (Treal)(-1.0))  + (Treal)1.0) / (Treal)2.0;
00346         ASSERTALWAYS(!(Interval<Treal>::intersect(homo, homoTmp).empty() && 
00347                        Interval<Treal>::intersect(lumo, lumoTmp).empty()));
00348         if (Interval<Treal>::intersect(homo, homoTmp).empty()) {
00349           // ok, lumo was computed in this iteration
00350           if (eigVecPtr)
00351             lumoWasComputed = 1;
00352           lumo.intersect(lumoTmp);
00353         }
00354         if (Interval<Treal>::intersect(lumo, lumoTmp).empty()) {
00355           // ok, homo was computed in this iteration
00356           if (eigVecPtr)
00357             homoWasComputed = 1;
00358           homo.intersect(homoTmp);
00359         }
00360         checkIntervals("PuriStepInfo::improveHomoLumo C");
00361       }
00362     }
00363   }
00364   
00365   template<typename Treal, typename Tvector, typename TdebugPolicy>
00366     void PuriStepInfo<Treal, Tvector, TdebugPolicy>::
00367     exchangeInfoWithNext(PuriStepInfo<Treal, Tvector, TdebugPolicy> & next) {
00368     Interval<Treal> zeroOneInt(0.0,1.0);
00369 
00370     next.checkIntervals("PuriStepInfo::exchangeInfoWithNext A");
00371 
00372     /* Improve homo/lumo and eig-bounds for next */
00373     Interval<Treal> homoForNext(homo);
00374     Interval<Treal> lumoForNext(lumo);
00375     Interval<Treal> eigIntForNext(eigInterval);
00376     homoForNext.puriStep(poly);
00377     lumoForNext.puriStep(poly);
00378     eigIntForNext.puriStep(poly);
00379     ASSERTALWAYS(!homoForNext.empty());
00380     ASSERTALWAYS(!lumoForNext.empty());
00381     ASSERTALWAYS(!eigIntForNext.empty());
00382     /* Increase intervals because relative precision in matrix-matrix 
00383      * multiplication can result in loss of accuracy in eigenvalues
00384      */
00385     homoForNext.increase(eigAccLoss);
00386     lumoForNext.increase(eigAccLoss);
00387     eigIntForNext.increase(eigAccLoss);
00388     /* Increase intervals because of thresholding */
00389     homoForNext.increase(next.actualThresh);
00390     lumoForNext.increase(next.actualThresh);
00391     homoForNext.intersect(zeroOneInt);
00392     lumoForNext.intersect(zeroOneInt);
00393     eigIntForNext.increase(next.actualThresh);
00394     next.improveEigInterval(eigIntForNext);
00395     next.improveHomoLumo(homoForNext, lumoForNext);
00396     
00397     /* Improve homo/lumo for this */
00398     /* FIXME: Consider improving also eigInterval from next
00399      * This could possibly only be done in one end of the interval since
00400      * for example information about negative eigenvalues is lost in case
00401      * of an x*x step. 
00402      */
00403     Interval<Treal> homoTmp(next.homo);
00404     Interval<Treal> lumoTmp(next.lumo);
00405     ASSERTALWAYS(!homoTmp.empty());
00406     ASSERTALWAYS(!lumoTmp.empty());
00407     /* Increase intervals because of thresholding. */
00408     homoTmp.increase(next.actualThresh);
00409     lumoTmp.increase(next.actualThresh);
00410     homoTmp.intersect(zeroOneInt);
00411     lumoTmp.intersect(zeroOneInt);
00412     homoTmp.invPuriStep(poly);
00413     lumoTmp.invPuriStep(poly);
00414     /* Increase intervals because relative precision in matrix-matrix 
00415      * multiplication can result in loss of accuracy in eigenvalues
00416      */
00417     homoTmp.increase(eigAccLoss);
00418     lumoTmp.increase(eigAccLoss);
00419 
00420     this->improveHomoLumo(homoTmp, lumoTmp);
00421   }
00422 
00423   template<typename Treal, typename Tvector, typename TdebugPolicy>
00424     Treal PuriStepInfo<Treal, Tvector, TdebugPolicy>::subspaceError() const {
00425     Interval<Treal> gap = Interval<Treal>(lumo.upp(),homo.low()); 
00426     if (actualThresh >= gap.length())
00427       return 1.0; /* 1.0 means no accuracy. */
00428     else {
00429       Treal error = actualThresh / (gap.length() - actualThresh); 
00430       return error < 1.0 ? error : (Treal)1.0;
00431     }
00432   }
00433 
00434   template<typename Treal, typename Tvector, typename TdebugPolicy>
00435     void PuriStepInfo<Treal, Tvector, TdebugPolicy>::
00436     improveEigInterval(Interval<Treal> const eInt) {
00437     eigInterval.intersect(eInt);
00438     Treal delta1 = eigInterval.upp() - 1;
00439     Treal delta0 = -eigInterval.low();
00440     delta = delta1 > delta0 ? delta1 : delta0;
00441     this->computen0n1();
00442   }
00443   
00444 
00445   template<typename Treal, typename Tvector, typename TdebugPolicy>
00446     void PuriStepInfo<Treal, Tvector, TdebugPolicy>::computen0n1() {
00447     Treal beta = 0.5; 
00448     /* Increase beta if possible */
00449     if (XmX2EuclNorm.upp() < 1 / (Treal)4)
00450       beta = (1 + template_blas_sqrt(1 - 4 * XmX2EuclNorm.upp())) / 2;
00451     n1 = (traceX2 - delta * (1 - beta) * n - 
00452           (1 - delta - beta) * traceX) / 
00453       ((1 + 2 * delta) * (delta + beta));
00454     n0 = (traceX2 + beta * (1 + delta) * n - 
00455           (1 + delta + beta) * traceX) / 
00456       ((1 + 2 * delta) * (delta + beta));
00457     if (n1 > nocc -1 && 
00458         n0 > n - nocc - 1)
00459       correctOccupation = 1;
00460   }
00461 
00466   template<typename Treal, typename Tvector, typename TdebugPolicy>
00467     void PuriStepInfo<Treal, Tvector, TdebugPolicy>::computeEigAccLoss() {
00468     Treal nnzPerRowX = nnzX / (Treal)n;
00469     Treal maxAbsErrPerElX2 = getRelPrecision<Treal>() * nnzPerRowX;
00470     /*  mah = max(abs(h_ij)) \approx relPrec * (nnz(X) / n)
00471      *  e is the exact eigenvalue of X^2, e' is the eigenvalue of the 
00472      *  computed X^2.
00473      *  | e - e' | <= || H ||_2 <=  || H ||_F = sqrt( sum h_ij^2 ) <=
00474      *  sqrt( mah^2 * nnz(X2))
00475      */
00476     eigAccLoss = maxAbsErrPerElX2 * template_blas_sqrt((Treal)nnzX2);
00477     ASSERTALWAYS(eigAccLoss >= 0);
00478   }
00479   
00480 
00481 
00482 } /* end namespace mat */
00483 #undef __ID__
00484 #endif

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