ergo
Lanczos.h
Go to the documentation of this file.
1/* Ergo, version 3.8.2, a program for linear scaling electronic structure
2 * calculations.
3 * Copyright (C) 2023 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
41#include "mat_utils.h"
42namespace 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;
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
98 Tricopy = Tri;
99 }
100 virtual ~Lanczos() {
101 delete[] v;
102 }
103 protected:
104 Tmatrix const & A;
105 Tvector* v;
109
110 Tvector r;
113 int j;
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)
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
Class for tridiagonal symmetric matrices.
Definition Failure.h:81
int capacity
Definition Lanczos.h:112
void step()
Definition Lanczos.h:126
virtual void run()
Definition Lanczos.h:87
MatrixTridiagSymmetric< Treal > Tri
Residual vector.
Definition Lanczos.h:111
Treal alpha
Definition Lanczos.h:121
void increaseCapacity(int const newCapacity)
Definition Lanczos.h:170
virtual ~Lanczos()
Definition Lanczos.h:100
Treal beta
Definition Lanczos.h:122
void restart(Tvector const &startVec)
Definition Lanczos.h:74
virtual void update()=0
void copyTridiag(MatrixTridiagSymmetric< Treal > &Tricopy)
Definition Lanczos.h:97
virtual bool converged() const =0
Tvector r
Vectors spanning Krylov subspace.
Definition Lanczos.h:110
int maxIter
Current step.
Definition Lanczos.h:114
void getEigVector(Tvector &eigVec, Treal const *const eVecTri) const
Definition Lanczos.h:150
Tvector * v
Definition Lanczos.h:105
int j
Definition Lanczos.h:113
Lanczos(Tmatrix const &AA, Tvector const &startVec, int maxIt=100, int cap=100)
Definition Lanczos.h:61
Tmatrix const & A
Definition Lanczos.h:104
Tridiagonal symmetric matrix class template.
Definition MatrixTridiagSymmetric.h:47
Utilities used by other parts of the hierarchical matrix library.
Definition Lanczos.h:43
Definition allocate.cc:39
const Treal Matrix< Treal >::ONE
Definition Matrix.h:3448
Xtrans< TX > transpose(TX const &A)
Transposition.
Definition matrix_proxy.h:131
static Treal getRelPrecision()
Definition matInclude.h:152
Treal template_blas_sqrt(Treal x)