ergo
Lanczos.h
Go to the documentation of this file.
1 /* Ergo, version 3.4, a program for linear scaling electronic structure
2  * calculations.
3  * Copyright (C) 2014 Elias Rudberg, Emanuel H. Rubensson, and Pawel Salek.
4  *
5  * This program is free software: you can redistribute it and/or modify
6  * it under the terms of the GNU General Public License as published by
7  * the Free Software Foundation, either version 3 of the License, or
8  * (at your option) any later version.
9  *
10  * This program is distributed in the hope that it will be useful,
11  * but WITHOUT ANY WARRANTY; without even the implied warranty of
12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13  * GNU General Public License for more details.
14  *
15  * You should have received a copy of the GNU General Public License
16  * along with this program. If not, see <http://www.gnu.org/licenses/>.
17  *
18  * Primary academic reference:
19  * Kohn−Sham Density Functional Theory Electronic Structure Calculations
20  * with Linearly Scaling Computational Time and Memory Usage,
21  * Elias Rudberg, Emanuel H. Rubensson, and Pawel Salek,
22  * J. Chem. Theory Comput. 7, 340 (2011),
23  * <http://dx.doi.org/10.1021/ct100611z>
24  *
25  * For further information about Ergo, see <http://www.ergoscf.org>.
26  */
27 
36 #ifndef MAT_LANCZOS
37 #define MAT_LANCZOS
38 #include "MatrixTridiagSymmetric.h"
39 #include "mat_utils.h"
40 namespace mat { /* Matrix namespace */
41  namespace arn { /* Arnoldi type methods namespace */
42 
56  template<typename Treal, typename Tmatrix, typename Tvector>
57  class Lanczos {
58  public:
59  Lanczos(Tmatrix const & AA, Tvector const & startVec,
60  int maxIt = 100, int cap = 100)
61  : A(AA), v(new Tvector[cap]), capacity(cap), j(0), maxIter(maxIt),
62  alpha(0), beta(0) {
63  assert(cap > 1);
64  Treal const ONE = 1.0;
65  v[0] = startVec;
66  if(v[0].eucl() < template_blas_sqrt(getRelPrecision<Treal>())) {
67  v[0].rand();
68  }
69  v[0] *= (ONE / v[0].eucl());
70  r = v[0];
71  }
72  void restart(Tvector const & startVec) {
73  j = 0;
74  alpha = 0;
75  beta = 0;
76  delete[] v;
77  v = new Tvector[capacity];
78  Treal const ONE = 1.0;
79  v[0] = startVec;
80  v[0] *= (ONE / startVec.eucl());
81  r = startVec;
82  Tri.clear();
83  }
84 
85  virtual void run() {
86  do {
87  step();
88  update();
89  if (j > maxIter)
90  throw AcceptableMaxIter("Lanczos::run() did not converge within maxIter");
91  }
92  while (!converged());
93  }
94 
95  inline void copyTridiag(MatrixTridiagSymmetric<Treal> & Tricopy) {
96  Tricopy = Tri;
97  }
98  virtual ~Lanczos() {
99  delete[] v;
100  }
101  protected:
102  Tmatrix const & A;
103  Tvector* v;
108  Tvector r;
110  int capacity;
111  int j;
112  int maxIter;
113  void increaseCapacity(int const newCapacity);
114  void step();
115  void getEigVector(Tvector& eigVec, Treal const * const eVecTri) const;
116  virtual void update() = 0;
117  virtual bool converged() const = 0;
118  private:
119  Treal alpha;
120  Treal beta;
121  }; /* end class definition Lanczos */
122 
123  template<typename Treal, typename Tmatrix, typename Tvector>
125  if (j + 1 >= capacity)
127  Treal const ONE = 1.0;
128  A.matVecProd(r, v[j]); // r = A * v[j]
129  alpha = transpose(v[j]) * r;
130  r += (-alpha) * v[j];
131  if (j) {
132  r += (-beta) * v[j-1];
133  v[j-1].writeToFile();
134  }
135  beta = r.eucl();
136  v[j+1] = r;
137  v[j+1] *= ONE / beta;
138  Tri.increase(alpha, beta);
139  ++j;
140  }
141 
142  /* FIXME: If several eigenvectors are needed it is more optimal to
143  * compute all of them at once since then the krylov subspace vectors
144  * only need to be read from memory once.
145  */
146  template<typename Treal, typename Tmatrix, typename Tvector>
148  getEigVector(Tvector& eigVec, Treal const * const eVecTri) const {
149  if (j <= 1) {
150  eigVec = v[0];
151  }
152  else {
153  v[0].readFromFile();
154  eigVec = v[0];
155  v[0].writeToFile();
156  }
157  eigVec *= eVecTri[0];
158  for (int ind = 1; ind <= j - 2; ++ind) {
159  v[ind].readFromFile();
160  eigVec += eVecTri[ind] * v[ind];
161  v[ind].writeToFile();
162  }
163  eigVec += eVecTri[j-1] * v[j-1];
164  }
165 
166  template<typename Treal, typename Tmatrix, typename Tvector>
168  increaseCapacity(int const newCapacity) {
169  assert(newCapacity > capacity);
170  assert(j > 0);
171  capacity = newCapacity;
172  Tvector* new_v = new Tvector[capacity];
173  /* FIXME: Fix so that file is copied when operator= is called in Vector
174  * class
175  */
176  for (int ind = 0; ind <= j - 2; ind++){
177  v[ind].readFromFile();
178  new_v[ind] = v[ind];
179  new_v[ind].writeToFile();
180  }
181  for (int ind = j - 1; ind <= j; ind++){
182  new_v[ind] = v[ind];
183  }
184  delete[] v;
185  v = new_v;
186  }
187 
188 
189  } /* end namespace arn */
190 } /* end namespace mat */
191 #endif
Tvector * v
Definition: Lanczos.h:103
virtual void run()
Definition: Lanczos.h:85
int j
Definition: Lanczos.h:111
void step()
Definition: Lanczos.h:124
Treal beta
Definition: Lanczos.h:120
Definition: Failure.h:71
virtual bool converged() const =0
void copyTridiag(MatrixTridiagSymmetric< Treal > &Tricopy)
Definition: Lanczos.h:95
virtual void update()=0
Tridiagonal symmetric matrix class template.
Definition: MatrixTridiagSymmetric.h:45
Treal alpha
Definition: Lanczos.h:119
Class for tridiagonal symmetric matrices.
virtual ~Lanczos()
Definition: Lanczos.h:98
Definition: allocate.cc:30
MatrixTridiagSymmetric< Treal > Tri
Residual vector.
Definition: Lanczos.h:109
void getEigVector(Tvector &eigVec, Treal const *const eVecTri) const
Definition: Lanczos.h:148
int capacity
Definition: Lanczos.h:110
int maxIter
Current step.
Definition: Lanczos.h:112
Tmatrix const & A
Definition: Lanczos.h:102
void increaseCapacity(int const newCapacity)
Definition: Lanczos.h:168
Lanczos(Tmatrix const &AA, Tvector const &startVec, int maxIt=100, int cap=100)
Definition: Lanczos.h:59
void restart(Tvector const &startVec)
Definition: Lanczos.h:72
Tvector r
Vectors spanning Krylov subspace.
Definition: Lanczos.h:108
Class template for building Krylov subspaces with Lanczos.
Definition: Lanczos.h:57
Xtrans< TX > transpose(TX const &A)
Transposition.
Definition: matrix_proxy.h:129
Treal template_blas_sqrt(Treal x)