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 
44 namespace 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>
54  void do_lanczos_method(int n,
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
template_blas_sqrt
Treal template_blas_sqrt(Treal x)
realtype.h
Definition of the main floating-point datatype used; the ergo_real type.
LOG_CAT_ERROR
#define LOG_CAT_ERROR
Definition: output.h:47
ergo_real
double ergo_real
Definition: realtype.h:69
simple_lanczos.h
Simple implementation of the Lanczos method.
simple_lanczos
Definition: simple_lanczos.cc:43
simple_lanczos::do_lanczos_method
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_ptr
ergo_real * ergo_real_ptr
Definition: integral_matrix_wrappers.cc:519
A
#define A
simple_lanczos::simple_lanczos_normalize_vector
void simple_lanczos_normalize_vector(int n, ergo_real *v)
Definition: simple_lanczos.cc:52
mat::syev
static void syev(const char *jobz, const char *uplo, const int *n, T *a, const int *lda, T *w, T *work, const int *lwork, int *info)
Definition: mat_gblas.h:382
do_output
void do_output(int logCategory, int logArea, const char *format,...)
Definition: output.cc:53
simple_lanczos::simple_lanczos_get_eigs
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
simple_lanczos::simple_lanczos_get_vector_norm
ergo_real simple_lanczos_get_vector_norm(int n, const ergo_real *v)
Definition: simple_lanczos.cc:45
output.h
Functionality for writing output messages to a text file.
LOG_AREA_CI
#define LOG_AREA_CI
Definition: output.h:64