00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00037 #ifndef SIMPLE_LANCZOS_HEADER
00038 #define SIMPLE_LANCZOS_HEADER
00039
00040 #include "realtype.h"
00041 #include <vector>
00042 #include <cstdio>
00043
00044 namespace simple_lanczos {
00045
00046 ergo_real simple_lanczos_get_vector_norm(int n, const ergo_real* v);
00047 void simple_lanczos_normalize_vector(int n, ergo_real* v);
00048 void simple_lanczos_get_eigs(int n, ergo_real* M,
00049 ergo_real & currEig_lo, ergo_real* bestVector_lo,
00050 ergo_real & currEig_hi, ergo_real* bestVector_hi,
00051 ergo_real* eigValListResult);
00052
00053 template<typename Tmatvecmul>
00054 void do_lanczos_method(int n,
00055 const ergo_real* guessVector,
00056 ergo_real & resultEig_lo,
00057 ergo_real* resultVec_lo,
00058 ergo_real & resultEig_hi,
00059 ergo_real* resultVec_hi,
00060 const Tmatvecmul & matvecmul,
00061 int maxIterations_in,
00062 ergo_real shift,
00063 ergo_real extraEnergy) {
00064 if(n == 1) {
00065
00066 ergo_real tmpVec1[1];
00067 tmpVec1[0] = 1;
00068 ergo_real tmpVec2[1];
00069 tmpVec2[0] = 0;
00070 matvecmul.do_mat_vec_mul(n, tmpVec1, tmpVec2);
00071 ergo_real eigenValue = tmpVec2[0];
00072 resultEig_lo = eigenValue;
00073 resultEig_hi = eigenValue;
00074 resultVec_lo[0] = 1;
00075 resultVec_hi[0] = 1;
00076 }
00077 typedef ergo_real* ergo_real_ptr;
00078 int maxIterations = maxIterations_in;
00079 if(maxIterations > n)
00080 maxIterations = n;
00081 ergo_real** q = new ergo_real_ptr[n+1];
00082 q[0] = new ergo_real[n];
00083 for(int i = 0; i < n; i++)
00084 q[0][i] = 0;
00085 q[1] = new ergo_real[n];
00086 for(int i = 0; i < n; i++)
00087 q[1][i] = guessVector[i];
00088 simple_lanczos_normalize_vector(n, q[1]);
00089 std::vector<ergo_real> z(n);
00090 std::vector<ergo_real> alpha(n+1);
00091 std::vector<ergo_real> beta(n+1);
00092 beta[0] = 0;
00093 std::vector<ergo_real> bestVector_lo(maxIterations+1);
00094 std::vector<ergo_real> bestVector_hi(maxIterations+1);
00095 ergo_real currEig_lo = 0;
00096 ergo_real currEig_hi = 0;
00097 ergo_real curr_E_lo = 0;
00098 ergo_real curr_E_hi = 0;
00099 for(int j = 1; j <= maxIterations; j++) {
00100
00101 matvecmul.do_mat_vec_mul(n, q[j], &z[0]);
00102
00103 alpha[j] = 0;
00104 for(int i = 0; i < n; i++)
00105 alpha[j] += q[j][i] * z[i];
00106 for(int i = 0; i < n; i++)
00107 z[i] = z[i] - alpha[j] * q[j][i] - beta[j-1] * q[j-1][i];
00108 beta[j] = simple_lanczos_get_vector_norm(n, &z[0]);
00109 ergo_real* T = new ergo_real[j*j];
00110 for(int i = 0; i < j*j; i++)
00111 T[i] = 0;
00112 for(int i = 0; i < j; i++)
00113 T[i*j+i] = alpha[i+1];
00114 for(int i = 0; i < j-1; i++) {
00115 T[i*j+(i+1)] = beta[i+1];
00116 T[(i+1)*j+i] = beta[i+1];
00117 }
00118 simple_lanczos_get_eigs(j, T, currEig_lo, &bestVector_lo[0], currEig_hi, &bestVector_hi[0], NULL);
00119
00120 for(int k = 0; k < n; k++) {
00121 ergo_real sum = 0;
00122 for(int i = 1; i <= j; i++)
00123 sum += bestVector_lo[i-1] * q[i][k];
00124 resultVec_lo[k] = sum;
00125 }
00126
00127 for(int k = 0; k < n; k++) {
00128 ergo_real sum = 0;
00129 for(int i = 1; i <= j; i++)
00130 sum += bestVector_hi[i-1] * q[i][k];
00131 resultVec_hi[k] = sum;
00132 }
00133 curr_E_lo = currEig_lo + extraEnergy + shift;
00134 curr_E_hi = currEig_hi + extraEnergy + shift;
00135 if(beta[j] < 1e-11 && j > 1)
00136 break;
00137 if(j == maxIterations)
00138 break;
00139 q[j+1] = new ergo_real[n];
00140 for(int i = 0; i < n; i++)
00141 q[j+1][i] = z[i] / beta[j];
00142 }
00143 resultEig_lo = curr_E_lo;
00144 resultEig_hi = curr_E_hi;
00145 simple_lanczos_normalize_vector(n, resultVec_lo);
00146 simple_lanczos_normalize_vector(n, resultVec_hi);
00147 }
00148
00149 }
00150
00151 #endif