LanczosLargestMagnitudeEig.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 
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

Generated on 21 Nov 2012 for ergo by  doxygen 1.6.1