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_PURISTEPINFODEBUG 00037 #define MAT_PURISTEPINFODEBUG 00038 #include "matInclude.h" 00039 #include "Interval.h" 00040 #define __ID__ "$Id: PuriStepInfoDebug.h 4437 2012-07-05 09:01:18Z elias $" 00041 namespace mat { 00042 template<typename Treal, typename TdebugPolicy> 00043 class PuriStepInfoDebug : public TdebugPolicy { 00044 public: 00045 00046 inline void checkIntervals(Interval<Treal> const & eigInterval, 00047 Interval<Treal> const & homo, 00048 Interval<Treal> const & lumo, 00049 Interval<Treal> const & XmX2EuclNorm, 00050 const char* descriptionString) const {} 00051 template<typename Tmatrix> 00052 inline void computeExactValues(Tmatrix const & X, Tmatrix const & X2, 00053 int const n, int const nocc) {} 00054 protected: 00055 PuriStepInfoDebug(){} 00056 private: 00057 }; 00058 00059 00060 00061 /* Specialization for high debug level */ 00062 template<typename Treal> 00063 class PuriStepInfoDebug<Treal, DebugLevelHigh> : public DebugLevelHigh { 00064 public: 00065 void checkIntervals(Interval<Treal> const & eigInterval, 00066 Interval<Treal> const & homo, 00067 Interval<Treal> const & lumo, 00068 Interval<Treal> const & XmX2EuclNorm, 00069 const char* descriptionString) const; 00070 template<typename Tmatrix> 00071 void computeExactValues(Tmatrix const & X, Tmatrix const & X2, 00072 int const n, int const nocc); 00073 protected: 00074 PuriStepInfoDebug() 00075 :homoExact(), lumoExact(), 00076 lmaxExact(), lminExact(), 00077 XmX2EuclNormExact() 00078 {} 00079 Interval<Treal> homoExact; 00080 Interval<Treal> lumoExact; 00081 Interval<Treal> lmaxExact; 00082 Interval<Treal> lminExact; 00083 Interval<Treal> XmX2EuclNormExact; 00084 private: 00085 }; 00086 00087 template<typename Treal> 00088 void PuriStepInfoDebug<Treal, DebugLevelHigh>:: 00089 checkIntervals(Interval<Treal> const & eigInterval, 00090 Interval<Treal> const & homo, 00091 Interval<Treal> const & lumo, 00092 Interval<Treal> const & XmX2EuclNorm, 00093 const char* descriptionString) const { 00094 /* Failure here probably means that checkIntervals was called before 00095 * the exact values were calculated. 00096 */ 00097 if(lminExact.empty()) 00098 std::cout << "PuriStepInfoDebug::checkIntervals failed, descriptionString = '" 00099 << descriptionString << "'" << std::endl; 00100 ASSERTALWAYS(!lminExact.empty()); 00101 Interval<Treal> zeroOneInt(0.0,1.0); 00102 if (Interval<Treal>::intersect(eigInterval, lminExact).empty()) 00103 std::cout<<std::endl<<" eigInterval "<< 00104 std::setprecision(25)<<eigInterval<<std::endl 00105 <<" lminExact "<<std::setprecision(25) 00106 <<lminExact<<std::endl; 00107 ASSERTDEBUG(!Interval<Treal>::intersect(eigInterval, lminExact).empty()); 00108 if (Interval<Treal>::intersect(eigInterval, lmaxExact).empty()) 00109 std::cout<<std::endl<<" eigInterval "<< 00110 std::setprecision(25)<<eigInterval<<std::endl 00111 <<" lmaxExact "<<lmaxExact<<std::endl; 00112 ASSERTDEBUG(!Interval<Treal>::intersect(eigInterval, lmaxExact).empty()); 00113 /* The homo and lumo intervals only provides information if 00114 * the exact values lie in [0, 1] otherwise no check is done. 00115 */ 00116 if (!Interval<Treal>::intersect(zeroOneInt, homoExact).empty()) { 00117 if (Interval<Treal>::intersect(homo, homoExact).empty()) 00118 std::cout<<std::endl<<" homo "<< 00119 std::setprecision(25)<<homo<<std::endl 00120 <<" homoExact "<<homoExact<<std::endl; 00121 if(Interval<Treal>::intersect(homo, homoExact).empty()) 00122 std::cout << "PuriStepInfoDebug::checkIntervals failed due to intersect(homo, homoExact) , descriptionString = '" 00123 << descriptionString << "'" << std::endl; 00124 ASSERTDEBUG(!Interval<Treal>::intersect(homo, homoExact).empty()); 00125 } 00126 if (!Interval<Treal>::intersect(zeroOneInt, lumoExact).empty()) { 00127 if (Interval<Treal>::intersect(lumo, lumoExact).empty()) 00128 std::cout<<std::endl<<" lumo "<< 00129 std::setprecision(25)<<lumo<<std::endl 00130 <<" lumoExact "<<lumoExact<<std::endl; 00131 ASSERTDEBUG(!Interval<Treal>::intersect(lumo, lumoExact).empty()); 00132 } 00133 if (Interval<Treal>:: 00134 intersect(XmX2EuclNorm, XmX2EuclNormExact).empty()) { 00135 std::cout<<std::endl<<" XmX2EuclNorm "<< 00136 std::setprecision(25)<<XmX2EuclNorm<<std::endl 00137 <<" XmX2EuclNormExact "<<XmX2EuclNormExact<<std::endl; 00138 std::cout << "PuriStepInfoDebug::checkIntervals failed, descriptionString = '" 00139 << descriptionString << "'" << std::endl; 00140 } 00141 ASSERTDEBUG(!Interval<Treal>:: 00142 intersect(XmX2EuclNorm, XmX2EuclNormExact).empty()); 00143 } 00144 00145 template<typename Treal> 00146 template<typename Tmatrix> 00147 void PuriStepInfoDebug<Treal, DebugLevelHigh>:: 00148 computeExactValues(Tmatrix const & X, Tmatrix const & X2, 00149 int const n, int const nocc) { 00150 /* Set up work space */ 00151 std::vector<Treal> full(n*n); 00152 // Treal* full = new Treal[n*n]; 00153 Treal* eigvals = new Treal[n]; 00154 int lwork = 3*n-1; 00155 Treal* work = new Treal[lwork]; 00156 int info = 0; 00157 Treal euclNormMatrix; 00158 Treal precSyev; 00159 /* Compute lumoExact, homoExact, lminExact, and lmaxExact. */ 00160 X.fullMatrix(full); 00161 syev("N", "U", &n, &full[0], &n, eigvals, work, &lwork, &info); 00162 ASSERTALWAYS(!info); 00163 euclNormMatrix = template_blas_fabs(eigvals[0]) > template_blas_fabs(eigvals[n-1]) ? 00164 template_blas_fabs(eigvals[0]) : template_blas_fabs(eigvals[n-1]); 00165 precSyev = euclNormMatrix * getRelPrecision<Treal>(); 00166 lumoExact = Interval<Treal>(eigvals[n-nocc-1] - precSyev, 00167 eigvals[n-nocc-1] + precSyev); 00168 homoExact = Interval<Treal>(eigvals[n-nocc] - precSyev, 00169 eigvals[n-nocc] + precSyev); 00170 lminExact = Interval<Treal>(eigvals[0] - precSyev, 00171 eigvals[0] + precSyev); 00172 lmaxExact = Interval<Treal>(eigvals[n-1] - precSyev, 00173 eigvals[n-1] + precSyev); 00174 /* Compute largest magnitude eigenvalue of X-X^2 */ 00175 Tmatrix XmX2(X); 00176 XmX2 += -1.0 * X2; 00177 XmX2.fullMatrix(full); 00178 syev("N", "U", &n, &full[0], &n, eigvals, work, &lwork, &info); 00179 ASSERTALWAYS(!info); 00180 euclNormMatrix = template_blas_fabs(eigvals[0]) > template_blas_fabs(eigvals[n-1]) ? 00181 template_blas_fabs(eigvals[0]) : template_blas_fabs(eigvals[n-1]); 00182 precSyev = euclNormMatrix * getRelPrecision<Treal>(); 00183 00184 #if 0 00185 std::cout<<"EIGS XmX2: "<<std::endl; 00186 for(int ind = 0; ind < n; ++ind) 00187 std::cout<<std::setprecision(20)<<eigvals[ind]<<std::endl; 00188 std::cout<<"end EIGS XmX2: "<<std::endl<<std::endl; 00189 #endif 00190 00191 Treal lowBound = euclNormMatrix - precSyev > 0 ? 00192 euclNormMatrix - precSyev : 0; 00193 XmX2EuclNormExact = Interval<Treal>(lowBound, euclNormMatrix + precSyev); 00194 /* Free memory */ 00195 delete[] eigvals; 00196 delete[] work; 00197 } 00198 00199 00200 } /* end namespace mat */ 00201 #undef __ID__ 00202 #endif