37 #ifndef SIMPLE_LANCZOS_HEADER
38 #define SIMPLE_LANCZOS_HEADER
53 template<
typename Tmatvecmul>
60 const Tmatvecmul & matvecmul,
70 matvecmul.do_mat_vec_mul(n, tmpVec1, tmpVec2);
72 resultEig_lo = eigenValue;
73 resultEig_hi = eigenValue;
78 int maxIterations = maxIterations_in;
83 for(
int i = 0; i < n; i++)
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);
93 std::vector<ergo_real> bestVector_lo(maxIterations+1);
94 std::vector<ergo_real> bestVector_hi(maxIterations+1);
99 for(
int j = 1; j <= maxIterations; j++) {
101 matvecmul.do_mat_vec_mul(n, q[j], &z[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];
110 for(
int i = 0; i < j*j; i++)
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];
120 for(
int k = 0; k < n; k++) {
122 for(
int i = 1; i <= j; i++)
123 sum += bestVector_lo[i-1] * q[i][k];
124 resultVec_lo[k] = sum;
127 for(
int k = 0; k < n; k++) {
129 for(
int i = 1; i <= j; i++)
130 sum += bestVector_hi[i-1] * q[i][k];
131 resultVec_hi[k] = sum;
133 curr_E_lo = currEig_lo + extraEnergy + shift;
134 curr_E_hi = currEig_hi + extraEnergy + shift;
135 if(beta[j] < 1e-11 && j > 1)
137 if(j == maxIterations)
140 for(
int i = 0; i < n; i++)
141 q[j+1][i] = z[i] / beta[j];
143 resultEig_lo = curr_E_lo;
144 resultEig_hi = curr_E_hi;