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