ergo
simple_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
37#ifndef SIMPLE_LANCZOS_HEADER
38#define SIMPLE_LANCZOS_HEADER
39
40#include "realtype.h"
41#include <vector>
42#include <cstdio>
43
44namespace simple_lanczos {
45
48 void simple_lanczos_get_eigs(int n, ergo_real* M,
49 ergo_real & currEig_lo, ergo_real* bestVector_lo,
50 ergo_real & currEig_hi, ergo_real* bestVector_hi,
51 ergo_real* eigValListResult);
52
53 template<typename Tmatvecmul>
55 const ergo_real* guessVector,
56 ergo_real & resultEig_lo,
57 ergo_real* resultVec_lo,
58 ergo_real & resultEig_hi,
59 ergo_real* resultVec_hi,
60 const Tmatvecmul & matvecmul,
61 int maxIterations_in,
62 ergo_real shift,
63 ergo_real extraEnergy) {
64 if(n == 1) {
65 // Special case for n=1, in this case we need only one "matrix-vector" (really scalar) multiplication to get all info we need.
66 ergo_real tmpVec1[1];
67 tmpVec1[0] = 1;
68 ergo_real tmpVec2[1];
69 tmpVec2[0] = 0;
70 matvecmul.do_mat_vec_mul(n, tmpVec1, tmpVec2);
71 ergo_real eigenValue = tmpVec2[0];
72 resultEig_lo = eigenValue;
73 resultEig_hi = eigenValue;
74 resultVec_lo[0] = 1;
75 resultVec_hi[0] = 1;
76 }
77 typedef ergo_real* ergo_real_ptr;
78 int maxIterations = maxIterations_in;
79 if(maxIterations > n)
80 maxIterations = n;
81 ergo_real** q = new ergo_real_ptr[n+1];
82 q[0] = new ergo_real[n];
83 for(int i = 0; i < n; i++)
84 q[0][i] = 0;
85 q[1] = new ergo_real[n];
86 for(int i = 0; i < n; i++)
87 q[1][i] = guessVector[i];
89 std::vector<ergo_real> z(n);
90 std::vector<ergo_real> alpha(n+1);
91 std::vector<ergo_real> beta(n+1);
92 beta[0] = 0;
93 std::vector<ergo_real> bestVector_lo(maxIterations+1);
94 std::vector<ergo_real> bestVector_hi(maxIterations+1);
95 ergo_real currEig_lo = 0;
96 ergo_real currEig_hi = 0;
97 ergo_real curr_E_lo = 0;
98 ergo_real curr_E_hi = 0;
99 for(int j = 1; j <= maxIterations; j++) {
100 // Do matrix-vector multiplication
101 matvecmul.do_mat_vec_mul(n, q[j], &z[0]);
102 // OK, matrix-vector multiplication done
103 alpha[j] = 0;
104 for(int i = 0; i < n; i++)
105 alpha[j] += q[j][i] * z[i];
106 for(int i = 0; i < n; i++)
107 z[i] = z[i] - alpha[j] * q[j][i] - beta[j-1] * q[j-1][i];
108 beta[j] = simple_lanczos_get_vector_norm(n, &z[0]);
109 ergo_real* T = new ergo_real[j*j];
110 for(int i = 0; i < j*j; i++)
111 T[i] = 0;
112 for(int i = 0; i < j; i++)
113 T[i*j+i] = alpha[i+1];
114 for(int i = 0; i < j-1; i++) {
115 T[i*j+(i+1)] = beta[i+1];
116 T[(i+1)*j+i] = beta[i+1];
117 }
118 simple_lanczos_get_eigs(j, T, currEig_lo, &bestVector_lo[0], currEig_hi, &bestVector_hi[0], NULL);
119 // Set resultVec_lo
120 for(int k = 0; k < n; k++) {
121 ergo_real sum = 0;
122 for(int i = 1; i <= j; i++)
123 sum += bestVector_lo[i-1] * q[i][k];
124 resultVec_lo[k] = sum;
125 }
126 // Set resultVec_hi
127 for(int k = 0; k < n; k++) {
128 ergo_real sum = 0;
129 for(int i = 1; i <= j; i++)
130 sum += bestVector_hi[i-1] * q[i][k];
131 resultVec_hi[k] = sum;
132 }
133 curr_E_lo = currEig_lo + extraEnergy + shift;
134 curr_E_hi = currEig_hi + extraEnergy + shift;
135 if(beta[j] < 1e-11 && j > 1)
136 break;
137 if(j == maxIterations)
138 break;
139 q[j+1] = new ergo_real[n];
140 for(int i = 0; i < n; i++)
141 q[j+1][i] = z[i] / beta[j];
142 } // end for j
143 resultEig_lo = curr_E_lo;
144 resultEig_hi = curr_E_hi;
145 simple_lanczos_normalize_vector(n, resultVec_lo);
146 simple_lanczos_normalize_vector(n, resultVec_hi);
147 }
148
149} // end namespace simple_lanczos
150
151#endif
ergo_real * ergo_real_ptr
Definition: integral_matrix_wrappers.cc:519
Definition: simple_lanczos.cc:43
void do_lanczos_method(int n, const ergo_real *guessVector, ergo_real &resultEig_lo, ergo_real *resultVec_lo, ergo_real &resultEig_hi, ergo_real *resultVec_hi, const Tmatvecmul &matvecmul, int maxIterations_in, ergo_real shift, ergo_real extraEnergy)
Definition: simple_lanczos.h:54
ergo_real simple_lanczos_get_vector_norm(int n, const ergo_real *v)
Definition: simple_lanczos.cc:45
void simple_lanczos_normalize_vector(int n, ergo_real *v)
Definition: simple_lanczos.cc:52
void simple_lanczos_get_eigs(int n, ergo_real *M, ergo_real &currEig_lo, ergo_real *bestVector_lo, ergo_real &currEig_hi, ergo_real *bestVector_hi, ergo_real *eigValListResult)
Definition: simple_lanczos.cc:58
Definition of the main floating-point datatype used; the ergo_real type.
double ergo_real
Definition: realtype.h:69