Purification.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_PURIFICATION
00037 #define MAT_PURIFICATION
00038 #include <math.h>
00039 #include <iomanip>
00040 #define __ID__ "$Id: Purification.h 4437 2012-07-05 09:01:18Z elias $"
00041 namespace mat {
00042   template<typename Treal, typename Tmatrix, typename TdebugPolicy>
00043     class Purification :public TdebugPolicy {
00044   public: 
00045     typedef typename Tmatrix::VectorType VectorType;
00046     Purification(Tmatrix& M, 
00049                  normType const normXmX2, 
00050                  PuriInfo<Treal, VectorType, TdebugPolicy>& info
00057                  ); 
00058     void step();    
00059     void purify();
00060   protected:
00061     Tmatrix & X;
00062     Tmatrix X2;
00063     normType const normTruncation;
00064     normType const normXmX2;
00065     PuriInfo<Treal, VectorType, TdebugPolicy>& info;
00066     int niter;
00067     
00068     void stepComputeInfo(PuriStepInfo<Treal, VectorType, TdebugPolicy> & currentStep);
00069     
00070   private:
00071   };
00072   
00073   template<typename Treal, typename Tmatrix, typename TdebugPolicy>
00074     Purification<Treal, Tmatrix, TdebugPolicy>::
00075     Purification(Tmatrix& M, 
00076                  normType const normXmX2_inp, 
00077                  PuriInfo<Treal, VectorType, TdebugPolicy>& infop)
00078     :X(M), X2(M), normTruncation( infop.getNormForTruncation() ), normXmX2(normXmX2_inp), info(infop),niter(0) {
00079     Time tStepTotal;
00080     tStepTotal.tic();
00081 
00082     Treal lmin(0);
00083     Treal lmax(0);
00084     PuriStepInfo<Treal, VectorType, TdebugPolicy> & currentStep = info.getNext();
00085     currentStep.improveEigInterval(Interval<Treal>(0.0,1.0));
00086     Interval<Treal> eigFInt = info.getEigFInterval();
00087     lmin = eigFInt.low();
00088     lmax = eigFInt.upp();
00089     X.add_identity(-lmax);      /* Scale to [0, 1] interval and negate */
00090     X *= ((Treal)1.0 / (lmin - lmax));
00091     /* Compute transformed homo and lumo eigenvalues. */
00092     Interval<Treal> homo = info.getHomoF();
00093     Interval<Treal> lumo = info.getLumoF();
00094     homo = (homo - lmax) / (lmin - lmax);
00095     lumo = (lumo - lmax) / (lmin - lmax);
00096 
00097     Treal chosenThresh = info.getOptimalThresh();
00098     currentStep.setChosenThresh(chosenThresh);
00099     /* Consider convergence of eigs in getOptimalThresh */
00100     currentStep.setMemUsageBeforeTrunc();
00101     Time t;
00102     t.tic();
00103     Treal actualThresh = X.thresh(chosenThresh, normTruncation);
00104     currentStep.setTimeThresh(t.toc());
00105     currentStep.setActualThresh(actualThresh);
00106     if (homo.empty() || lumo.empty()) {
00107       throw Failure("Purification<Treal, Tmatrix, TdebugPolicy>::"
00108                     "Purification(Tmatrix&, normType const, " 
00109                     "PuriInfo<Treal, VectorType, TdebugPolicy>&): "
00110                     "HOMO or LUMO empty in purification constructor. ");
00111     }
00112     
00113     homo.increase(actualThresh);
00114     lumo.increase(actualThresh);
00115     stepComputeInfo(currentStep);
00116     currentStep.improveHomoLumo(homo,lumo);
00117     currentStep.setTimeTotal(tStepTotal.toc());
00118   } 
00120     template<typename Treal, typename Tmatrix, typename TdebugPolicy>
00121       void Purification<Treal, Tmatrix, TdebugPolicy>::step() {
00122       Time tStepTotal;
00123       tStepTotal.tic();
00124       PuriStepInfo<Treal, VectorType, TdebugPolicy> & currentStep = info.getNext();
00125       /* It may happen that X2 has many more nonzeros than X, for
00126          example 5 times as many.  Therefore it makes sense to try
00127          having only one "big" matrix in memory at a time. However,
00128          file operations have proved to be quite expensive and should
00129          be avoided if possible. Hence we want to achieve having only
00130          one big matrix in memory without unnecessary file
00131          operations. We are currently hoping that it will be ok to add
00132          a "small" matrix to a "big" one, that the memory usage after
00133          that operation will be like the memory usage for one big
00134          matrix + one small matrix. Therefore we are adding X to X2 (X
00135          is truncated, a "small" matrix) instead of the opposite when
00136          the 2*X-X*X polynomial is evaluated.
00137       */
00138       if (info(niter).getPoly()) {
00139         /*  Polynomial 2 * x - x * x  */
00140         X2 *= (Treal)-1.0;
00141         X2 += (Treal)2.0 * X;
00142         // Now X2 contains 2*X-X*X
00143       }
00144       /* In case of polynomial x * x we only need to transfer the
00145          content of X2 to X.  */
00146       // Transfer content of X2 to X clearing previous content of X if any
00147       // In current implementation this is needed regardless of which
00148       // polynomial is used
00149       X2.transfer(X); 
00150 
00151      /* Obtain homo/lumo information from previous before thresh choice. 
00152        *  Default value of 0.0 for thresh is used. getOptimalThresh can
00153        *  try different thresh values.
00154        */
00155       Treal chosenThresh = info.getOptimalThresh();
00156       currentStep.setChosenThresh(chosenThresh);
00157       /* Consider convergence of eigs in getOptimalThresh? */
00158       currentStep.setMemUsageBeforeTrunc();
00159       Time t;
00160       t.tic();
00161       currentStep.setActualThresh(X.thresh(chosenThresh, normTruncation));
00162       currentStep.setTimeThresh(t.toc());
00163       stepComputeInfo(currentStep);
00164       if (niter > 5 &&
00165           currentStep.getHomo().low() > 0.9 &&
00166           currentStep.getLumo().upp() < 0.1 &&
00167           info(niter-5).getHomo().low() > 0.9 &&
00168           info(niter-5).getLumo().upp() < 0.1 &&
00169           ((1 - currentStep.getHomo().low()) > 
00170            (1 - info(niter-5).getHomo().low()) - 
00171            currentStep.getEigAccLoss()  ||
00172            currentStep.getLumo().upp() >
00173            info(niter-5).getLumo().upp() - 
00174            currentStep.getEigAccLoss())) {      
00175         throw Failure("Purification<Treal, Tmatrix, TdebugPolicy>"
00176                       "::step(): HOMO-LUMO gap has not increased in the "
00177                       "last five iterations. Desired accuracy too high for"
00178                       "current floating point precision?");
00179       }
00180       
00181       ++niter;
00182       currentStep.setTimeTotal(tStepTotal.toc());
00183     }
00184     
00185     template<typename Treal, typename Tmatrix, typename TdebugPolicy>
00186       void Purification<Treal, Tmatrix, TdebugPolicy>::purify() {
00187       
00188       while (niter < info.getMaxSteps() - 1 && !info.converged() ) {
00189         step();
00190       }
00191       if ( info.converged() ) {
00192         // Only if converged because forcing correct occupation can have 
00193         // strange effects otherwise.
00194         info.forceCorrectOccupation();  
00195         info.improveCorrectOccupation();
00196         info.improveInfo();
00197         info.improveHomoLumoF();
00198       }
00199     }
00200 
00201     template<typename Treal, typename Tmatrix, typename TdebugPolicy>
00202       void Purification<Treal, Tmatrix, TdebugPolicy>::
00203       stepComputeInfo(PuriStepInfo<Treal, VectorType, TdebugPolicy> & currentStep) {
00204       Treal const XmX2ENIsSmallValue = 0.207106781186547;
00205       Treal const ONE = 1.0;
00206       
00207       Time t;
00208       t.tic();
00209       X2 = ONE * X * X;
00210       currentStep.setTimeSquare(t.toc());
00211   
00212       currentStep.setNnzX(X.nnz());
00213       currentStep.setNnzX2(X2.nnz());
00214       currentStep.computeEigAccLoss();
00215       currentStep.computeExactValues(X, X2);
00216       currentStep.setTraceX(X.trace());
00217       currentStep.setTraceX2(X2.trace());
00218 
00219       /* Now we are about to compute the Euclidean norm of (X -
00220          X2). Previously this operation needed lots of memory. Now,
00221          however, we hope that the memory usage is smaller because the
00222          difference matrix is never explicitly computed. */
00223       
00224       typename Tmatrix::VectorType * eigVecPtr = new typename Tmatrix::VectorType;
00225       Treal diffAcc;
00226       Interval<Treal> XmX2EN;
00227       t.tic();
00228       if (info.ShouldComputeXmX2EuclNormAccurately(diffAcc)) {
00229         if (normXmX2 == euclNorm)
00230           XmX2EN = Tmatrix::euclDiffIfSmall(X, X2,  
00231                                             diffAcc, 
00232                                             XmX2ENIsSmallValue,
00233                                             eigVecPtr);   
00234         else
00235           XmX2EN = Tmatrix::diffIfSmall(X, X2,  
00236                                         normXmX2, diffAcc, 
00237                                         XmX2ENIsSmallValue);
00238       }
00239       else {
00240         XmX2EN = Tmatrix::diffIfSmall(X, X2,  
00241                                       frobNorm, diffAcc, 
00242                                       XmX2ENIsSmallValue);
00243       }
00244       currentStep.setTimeXmX2Norm(t.toc());
00245 
00246       if (!eigVecPtr->is_empty())
00247         currentStep.setEigVecPtr(eigVecPtr);
00248       else
00249         delete eigVecPtr;
00250 
00251       XmX2EN.increase(currentStep.getEigAccLoss());
00252       Interval<Treal> zeroInt(0.0, XmX2EN.upp());
00253       XmX2EN.intersect(zeroInt);
00254       // FIXME: 
00255       // Consider improving lanczos so that if to good accuracy is requested
00256       // the best possible accuracy is returned
00257       currentStep.setXmX2EuclNorm(XmX2EN);
00258       currentStep.checkIntervals("Purification::stepComputeInfo after setXmX2EuclNorm.");
00259 
00260       Treal tmpVal = template_blas_sqrt(1 + 4 * XmX2EN.upp());
00261       currentStep.improveEigInterval
00262         (Interval<Treal>((1 - tmpVal) / 2, (1 + tmpVal) / 2));
00263 
00264       currentStep.checkIntervals("Purification::stepComputeInfo B");
00265       
00266       info.improveInfo();
00267       if (currentStep.getEigInterval().length() > 1.5)
00268         throw Failure("Purification<Treal, Tmatrix, TdebugPolicy>"
00269                       "::step() : "
00270                       "Eigenvalues moved to far from [0, 1] interval.");
00271       /* Now decide which polynomial to use for first step. */
00272       currentStep.setPoly();
00273       currentStep.checkIntervals("Purification::stepComputeInfo C");
00274     }
00275 
00276 
00277 
00278 
00279 } /* end namespace mat */
00280 #undef __ID__
00281 #endif

Generated on Mon Sep 17 14:30:40 2012 for ergo by  doxygen 1.4.7