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