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 00037 #ifndef MAT_LANCZOSLARGESTMAGNITUDEEIG 00038 #define MAT_LANCZOSLARGESTMAGNITUDEEIG 00039 #include <limits> 00040 #include "Lanczos.h" 00041 namespace mat { /* Matrix namespace */ 00042 namespace arn { /* Arnoldi type methods namespace */ 00043 template<typename Treal, typename Tmatrix, typename Tvector> 00044 class LanczosLargestMagnitudeEig 00045 : public Lanczos<Treal, Tmatrix, Tvector> { 00046 public: 00047 LanczosLargestMagnitudeEig(Tmatrix const & AA, Tvector const & startVec, 00048 int maxIter = 100, int cap = 100) 00049 : Lanczos<Treal, Tmatrix, Tvector>(AA, startVec, maxIter, cap), 00050 eVal(0), acc(0), eigVectorTri(0), absTol( std::numeric_limits<Treal>::max() ), 00051 relTol(template_blas_sqrt(template_blas_sqrt(getRelPrecision<Treal>()))), 00052 eValTmp(0), accTmp(0) {} 00053 void setRelTol(Treal const newTol) { relTol = newTol; } 00054 void setAbsTol(Treal const newTol) { absTol = newTol; } 00055 00056 inline void getLargestMagnitudeEig(Treal& ev, Treal& accuracy) { 00057 ev = eVal; 00058 accuracy = acc; 00059 } 00060 void getLargestMagnitudeEigPair(Treal& eValue, 00061 Tvector& eVector, 00062 Treal& accuracy); 00063 00064 virtual void run() { 00065 Lanczos<Treal, Tmatrix, Tvector>::run(); 00066 computeEigVec(); 00067 eVal = eValTmp; 00068 acc = accTmp; 00069 rerun(); 00070 /* FIXME! Elias added these extra Lanczos reruns for small 00071 matrices to make the tests pass. This is bad. 00072 Elias note 2011-01-19: now added one more rerun() 00073 because otherwise the test failed when run on Mac. 00074 This is really bad. */ 00075 if(this->A.get_nrows() == 5) { 00076 rerun(); 00077 rerun(); 00078 rerun(); 00079 } 00080 } 00081 00082 void rerun() { 00083 #if 1 00084 /* Because of the statistical chance of misconvergence: 00085 * Compute new eigenpair with eigenvector in direction orthogonal 00086 * to the first eigenvector. 00087 */ 00088 Tvector newResidual(eVec); 00089 newResidual.rand(); 00090 Treal sP = transpose(eVec) * newResidual; 00091 newResidual += -sP * eVec; 00092 this->restart(newResidual); 00093 Lanczos<Treal, Tmatrix, Tvector>::run(); 00094 /* If the new eigenvalue has larger magnitude: 00095 * Use those instead 00096 */ 00097 if (template_blas_fabs(eValTmp) > template_blas_fabs(eVal)) { 00098 eVal = eValTmp; 00099 acc = accTmp; 00100 computeEigVec(); 00101 } 00102 // Guard against unrealistically small accuracies 00103 acc = acc / template_blas_fabs(eVal) > 2 * std::numeric_limits<Treal>::epsilon() ? acc : 2 * template_blas_fabs(eVal) * std::numeric_limits<Treal>::epsilon(); 00104 #endif 00105 } 00106 00107 virtual ~LanczosLargestMagnitudeEig() { 00108 delete[] eigVectorTri; 00109 } 00110 protected: 00111 Treal eVal; 00112 Tvector eVec; 00113 Treal acc; 00114 Treal* eigVectorTri; 00117 Treal absTol; 00118 Treal relTol; 00119 void computeEigenPairTri(); 00120 void computeEigVec(); 00121 virtual void update() { 00122 computeEigenPairTri(); 00123 } 00124 virtual bool converged() const; 00125 Treal eValTmp; 00126 Treal accTmp; 00127 private: 00128 }; 00129 00130 template<typename Treal, typename Tmatrix, typename Tvector> 00131 void LanczosLargestMagnitudeEig<Treal, Tmatrix, Tvector>:: 00132 getLargestMagnitudeEigPair(Treal& eValue, 00133 Tvector& eVector, 00134 Treal& accuracy) { 00135 eValue = eVal; 00136 accuracy = acc; 00137 eVector = eVec; 00138 } 00139 00140 00141 template<typename Treal, typename Tmatrix, typename Tvector> 00142 void LanczosLargestMagnitudeEig<Treal, Tmatrix, Tvector>:: 00143 computeEigenPairTri() { 00144 delete[] eigVectorTri; 00145 Treal* eigVectorMax = new Treal[this->j]; 00146 Treal* eigVectorMin = new Treal[this->j]; 00147 00148 /* Get largest eigenvalue. */ 00149 int const lMax = this->j - 1; 00150 Treal eValMax; 00151 Treal accMax; 00152 this->Tri.getEigsByIndex(&eValMax, eigVectorMax, &accMax, 00153 lMax, lMax); 00154 /* Get smallest eigenvalue. */ 00155 int const lMin = 0; 00156 Treal eValMin; 00157 Treal accMin; 00158 this->Tri.getEigsByIndex(&eValMin, eigVectorMin, &accMin, 00159 lMin, lMin); 00160 if (template_blas_fabs(eValMin) > template_blas_fabs(eValMax)) { 00161 eValTmp = eValMin; 00162 accTmp = accMin; 00163 delete[] eigVectorMax; 00164 eigVectorTri = eigVectorMin; 00165 } 00166 else { 00167 eValTmp = eValMax; 00168 accTmp = accMax; 00169 delete[] eigVectorMin; 00170 eigVectorTri = eigVectorMax; 00171 } 00172 } 00173 00174 00175 template<typename Treal, typename Tmatrix, typename Tvector> 00176 void LanczosLargestMagnitudeEig<Treal, Tmatrix, Tvector>:: 00177 computeEigVec() { 00178 /* Compute eigenvector as a linear combination of Krylov vectors */ 00179 assert(eigVectorTri); 00180 this->getEigVector(eVec, eigVectorTri); 00181 } 00182 00183 00184 template<typename Treal, typename Tmatrix, typename Tvector> 00185 bool LanczosLargestMagnitudeEig<Treal, Tmatrix, Tvector>:: 00186 converged() const { 00187 bool conv = accTmp < absTol; /* Absolute accuracy */ 00188 if (template_blas_fabs(eValTmp) > 0) { 00189 conv = conv && 00190 accTmp / template_blas_fabs(eValTmp) < relTol; /* Relative acc.*/ 00191 } 00192 return conv; 00193 } 00194 00195 00196 00197 00198 #if 1 00199 template<typename Treal, typename Tmatrix, typename Tvector> 00200 class LanczosLargestMagnitudeEigIfSmall 00201 : public LanczosLargestMagnitudeEig<Treal, Tmatrix, Tvector> { 00202 public: 00203 LanczosLargestMagnitudeEigIfSmall 00204 (Tmatrix const & AA, Tvector const & startVec, 00205 Treal const maxAbsVal, 00206 int maxIter = 100, int cap = 100) 00207 : LanczosLargestMagnitudeEig<Treal, Tmatrix, Tvector> 00208 (AA, startVec, maxIter, cap), maxAbsValue(maxAbsVal) { 00209 } 00210 bool largestMagEigIsSmall() {return eigIsSmall;} 00211 00212 virtual void run() { 00213 Lanczos<Treal, Tmatrix, Tvector>::run(); 00214 this->computeEigVec(); 00215 this->eVal = this->eValTmp; 00216 this->acc = this->accTmp; 00217 if (eigIsSmall) /* No need to rerun if eigenvalue is large. */ 00218 this->rerun(); 00219 } 00220 00221 protected: 00222 Treal const maxAbsValue; 00223 bool eigIsSmall; 00224 virtual void update() { 00225 LanczosLargestMagnitudeEig<Treal, Tmatrix, Tvector>::update(); 00226 eigIsSmall = template_blas_fabs(this->eValTmp) < maxAbsValue; 00227 } 00228 virtual bool converged() const; 00229 private: 00230 }; 00231 00232 template<typename Treal, typename Tmatrix, typename Tvector> 00233 bool LanczosLargestMagnitudeEigIfSmall<Treal, Tmatrix, Tvector>:: 00234 converged() const { 00235 bool convAccuracy = 00236 LanczosLargestMagnitudeEig<Treal, Tmatrix, Tvector>::converged(); 00237 return convAccuracy || (!eigIsSmall); 00238 } 00239 00240 #endif 00241 00242 } /* end namespace arn */ 00243 00245 00246 // EMANUEL COMMENT: 00247 // A function similar to euclIfSmall below lies/used to lie inside the MatrixSymmetric class. 00248 // There, the matrix was copied and truncated before the calculation. 00249 // It is unclear if that had a significant positive impact on the execution time - it definitely required more memory. 00256 template<typename Tmatrix, typename Treal> 00257 Interval<Treal> euclIfSmall(Tmatrix const & M, 00258 Treal const requestedAbsAccuracy, 00259 Treal const requestedRelAccuracy, 00260 Treal const maxAbsVal, 00261 typename Tmatrix::VectorType * const eVecPtr = 0) { 00262 assert(requestedAbsAccuracy >= 0); 00263 // Treal frobTmp; 00264 /* Check if norm is really small, or can easily be determined to be 'large', in those cases quick return */ 00265 Treal euclLowerBound; 00266 Treal euclUpperBound; 00267 M.quickEuclBounds(euclLowerBound, euclUpperBound); 00268 if ( euclUpperBound < requestedAbsAccuracy ) 00269 // Norm is really small, quick return 00270 return Interval<Treal>( euclLowerBound, euclUpperBound ); 00271 if ( euclLowerBound > maxAbsVal ) 00272 // Norm is not small, quick return 00273 return Interval<Treal>( euclLowerBound, euclUpperBound ); 00274 int maxIter = M.get_nrows() * 100; 00275 typename Tmatrix::VectorType guess; 00276 SizesAndBlocks cols; 00277 M.getCols(cols); 00278 guess.resetSizesAndBlocks(cols); 00279 guess.rand(); 00280 mat::arn::LanczosLargestMagnitudeEigIfSmall<Treal, Tmatrix, typename Tmatrix::VectorType> 00281 lan(M, guess, maxAbsVal, maxIter); 00282 lan.setAbsTol( requestedAbsAccuracy ); 00283 lan.setRelTol( requestedRelAccuracy ); 00284 lan.run(); 00285 Treal eVal = 0; 00286 Treal acc = 0; 00287 Treal eValMin = 0; 00288 if (lan.largestMagEigIsSmall()) { 00289 if (eVecPtr) 00290 lan.getLargestMagnitudeEigPair(eVal, *eVecPtr, acc); 00291 else 00292 lan.getLargestMagnitudeEig(eVal, acc); 00293 eVal = template_blas_fabs(eVal); 00294 eValMin = eVal - acc; 00295 eValMin = eValMin >= 0 ? eValMin : (Treal)0; 00296 return Interval<Treal>(eValMin, eVal + acc); 00297 } 00298 else { 00299 eValMin = euclLowerBound; 00300 eValMin = eValMin >= maxAbsVal ? eValMin : maxAbsVal; 00301 return Interval<Treal>(eValMin, euclUpperBound); 00302 } 00303 } 00304 00305 } /* end namespace mat */ 00306 #endif