ergo
Lanczos.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 
00036 #ifndef MAT_LANCZOS
00037 #define MAT_LANCZOS
00038 #include "MatrixTridiagSymmetric.h"
00039 #include "mat_utils.h"
00040 namespace mat { /* Matrix namespace */
00041   namespace arn { /* Arnoldi type methods namespace */
00042 
00056     template<typename Treal, typename Tmatrix, typename Tvector>
00057       class Lanczos {
00058     public:
00059       Lanczos(Tmatrix const & AA, Tvector const & startVec, 
00060               int maxIt = 100, int cap = 100) 
00061         : A(AA), v(new Tvector[cap]), capacity(cap), j(0), maxIter(maxIt),    
00062         alpha(0), beta(0) {
00063         assert(cap > 1);
00064         Treal const ONE = 1.0;    
00065         v[0] = startVec;
00066         if(v[0].eucl() < template_blas_sqrt(getRelPrecision<Treal>())) {
00067           v[0].rand();
00068         }
00069         v[0] *= (ONE / v[0].eucl());
00070         r = v[0];
00071       }
00072         void restart(Tvector const & startVec) {
00073           j = 0;
00074           alpha = 0;
00075           beta = 0;
00076           delete[] v;
00077           v = new Tvector[capacity];
00078           Treal const ONE = 1.0;
00079           v[0] = startVec;
00080           v[0] *= (ONE / startVec.eucl());
00081           r = startVec;   
00082           Tri.clear();
00083         }
00084         
00085         virtual void run() {
00086           do {
00087             step();
00088             update();
00089             if (j > maxIter)
00090               throw AcceptableMaxIter("Lanczos::run() did not converge within maxIter");
00091           }
00092           while (!converged());
00093         }
00094 
00095         inline void copyTridiag(MatrixTridiagSymmetric<Treal> & Tricopy) {
00096           Tricopy = Tri;
00097         }
00098         virtual ~Lanczos() {
00099           delete[] v;
00100         }
00101     protected:
00102       Tmatrix const & A;
00103       Tvector* v; 
00108       Tvector r; 
00109       MatrixTridiagSymmetric<Treal> Tri;
00110       int capacity;
00111       int j;     
00112       int maxIter;
00113       void increaseCapacity(int const newCapacity);
00114       void step();
00115       void getEigVector(Tvector& eigVec, Treal const * const eVecTri) const;
00116       virtual void update() = 0;
00117       virtual bool converged() const = 0;
00118     private:
00119       Treal alpha;
00120       Treal beta;
00121     }; /* end class definition Lanczos */
00122 
00123     template<typename Treal, typename Tmatrix, typename Tvector>
00124       void Lanczos<Treal, Tmatrix, Tvector>::step() {
00125       if (j + 1 >= capacity)
00126         increaseCapacity(capacity * 2);
00127       Treal const ONE = 1.0;
00128       A.matVecProd(r, v[j]);        // r = A * v[j] 
00129       alpha = transpose(v[j]) * r;
00130       r += (-alpha) * v[j];
00131       if (j) {
00132         r += (-beta) * v[j-1];
00133         v[j-1].writeToFile();
00134       }
00135       beta = r.eucl();
00136       v[j+1] = r;
00137       v[j+1] *= ONE / beta;
00138       Tri.increase(alpha, beta);
00139       ++j;
00140     }
00141     
00142     /*  FIXME: If several eigenvectors are needed it is more optimal to
00143      *  compute all of them at once since then the krylov subspace vectors
00144      *  only need to be read from memory once.
00145      */
00146     template<typename Treal, typename Tmatrix, typename Tvector>
00147       void Lanczos<Treal, Tmatrix, Tvector>::
00148       getEigVector(Tvector& eigVec, Treal const * const eVecTri) const {
00149       if (j <= 1) {
00150         eigVec = v[0];
00151       } 
00152       else {
00153         v[0].readFromFile();
00154         eigVec = v[0];
00155         v[0].writeToFile();
00156       }      
00157       eigVec *= eVecTri[0];
00158       for (int ind = 1; ind <= j - 2; ++ind) {
00159         v[ind].readFromFile();
00160         eigVec += eVecTri[ind] * v[ind];
00161         v[ind].writeToFile();
00162       }
00163       eigVec += eVecTri[j-1] * v[j-1];
00164     }
00165 
00166     template<typename Treal, typename Tmatrix, typename Tvector>
00167       void Lanczos<Treal, Tmatrix, Tvector>::
00168       increaseCapacity(int const newCapacity) {
00169       assert(newCapacity > capacity);
00170       assert(j > 0);
00171       capacity = newCapacity;
00172       Tvector* new_v = new Tvector[capacity];
00173       /* FIXME: Fix so that file is copied when operator= is called in Vector
00174        * class
00175        */
00176       for (int ind = 0; ind <= j - 2; ind++){
00177         v[ind].readFromFile();
00178         new_v[ind] = v[ind];
00179         new_v[ind].writeToFile();
00180       }
00181       for (int ind = j - 1; ind <= j; ind++){
00182         new_v[ind] = v[ind];
00183       }
00184       delete[] v;
00185       v = new_v;
00186     }
00187 
00188 
00189   } /* end namespace arn */
00190 } /* end namespace mat */
00191 #endif