00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00025 template<class fieldT>
00026 class Lanczos {
00027 public:
00028 static mdp_complex step(fieldT &psi,
00029 gauge_field &U,
00030 coefficients &coeff,
00031 bool force=false,
00032 bool output_check=false) {
00033 begin_function("Lanczos__step");
00034 static bool init=true;
00035 static fieldT p(psi);
00036 static fieldT q(psi);
00037 static fieldT r(psi);
00038 static double alpha, beta;
00039 static site x(psi.lattice());
00040
00041 if(init || force) {
00042 mdp << "Initializing Lanczos vectors\n";
00043
00044
00045
00046
00047
00048 psi.update();
00049 q=psi;
00050 double norm = sqrt(norm_square(q));
00051 q /= norm;
00052 p = 0.0;
00053 r = 0.0;
00054 beta = 1.0;
00055 init=false;
00056 }
00057
00058
00059
00060
00061
00062 r=p;
00063 q/=beta;
00064 p=q;
00065 r*=-beta;
00066 q=r;
00067
00068 p.update();
00069 mul_Q(r,p,U,coeff);
00070 multiply_by_gamma5(r,r);
00071 q += r;
00072 alpha = real_scalar_product(p,q);
00073 mdp_add_scaled_field(q,-alpha,p);
00074 beta = sqrt(norm_square(q));
00075
00076 if(output_check) {
00077
00078 static fieldT s(psi);
00079 double pp, qq;
00080 mdp_complex pq, qr, ps;
00081 site x(psi.lattice());
00082
00083 pp=norm_square(p);
00084 qq=norm_square(q);
00085 pq=p*q;
00086
00087 mul_Q(s,q,U,coeff);
00088 multiply_by_gamma5(s,s);
00089 s.update();
00090 qr=q*r;
00091
00092 mul_Q(r,p,U,coeff);
00093 multiply_by_gamma5(r,r);
00094 r.update();
00095 ps=p*s;
00096
00097 mdp << "|p|=" << pp << '\n';
00098 mdp << "|q|=" << qq << '\n';
00099 mdp << "<p|q>=" << pq << '\n';
00100 mdp << "<q|Q|p>=" << qr << '\n';
00101 mdp << "<p|Q|q>^*=" << ps << '\n';
00102 mdp << "alpha=" << alpha << '\n';
00103 mdp << "beta=" << beta << '\n';
00104
00105 }
00106 end_function("Lanczos__step");
00107 return mdp_complex(alpha, beta);
00108 }
00109 };
00110
00111
00112
00113
00114
00115
00116
00117
00118
00119
00120
00121
00122
00123
00124