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_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