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