00001
00002
00003
00004
00005
00006
00007
00033
00034 class CG2 {
00035 public:
00036 template <class fieldT, class fieldG>
00037 static inversion_stats inverter(fieldT &psi_out,
00038 fieldT &psi_in,
00039 fieldG &U,
00040 coefficients &coeff,
00041 mdp_real absolute_precision=mdp_precision,
00042 mdp_real relative_precision=0,
00043 int max_steps=2000,
00044 bool qdaggerq=false) {
00045 mpi.begin_function("ConjugateGradientInverter");
00046 int step=0;
00047 double residue, rresidue=-1, old_rresidue;
00048 double time=mpi.time();
00049 inversion_stats stats;
00050
00051 fieldT r(psi_in);
00052 fieldT rnew(psi_in);
00053 fieldT p(psi_in);
00054 fieldT pnew(psi_in);
00055 fieldT t(psi_in);
00056 fieldT s(psi_in);
00057 fieldT solnew(psi_in);
00058 fieldT b(psi_in);
00059
00060 double alpha, beta, rrtmp, rrtmp2, psdot;
00061
00062
00063
00064
00065 int mulQ=1, guesszero=1,debug=0;
00066
00067 psi_out=0;
00068 r=psi_in;
00069 p=r;
00070
00071
00072
00073
00074
00075
00076
00077 do {
00078 p.update();
00079 mul_Q(t,p,U,coeff); t/=2.0*coeff["kappa"]; t.update();
00080 coeff["sign"]=-1;
00081 mul_Q(s,t,U,coeff);
00082 coeff["sign"]=1;
00083 s/=2.0*coeff["kappa"];
00084 rrtmp=real_scalar_product(r,r);
00085 psdot=real_scalar_product(s,p);
00086 alpha=rrtmp/psdot;
00087 rnew=r;
00088 mdp_add_scaled_field(rnew, -alpha, s);
00089 solnew=psi_out;
00090 mdp_add_scaled_field(solnew, alpha, p);
00091 rrtmp2=real_scalar_product(rnew,rnew);
00092 beta=rrtmp2/rrtmp;
00093 pnew=rnew;
00094 mdp_add_scaled_field(pnew, beta, p);
00095 old_rresidue=rresidue;
00096 residue=sqrt(rrtmp2/(psi_in.global_size()));
00097 rresidue=relative_residue(rnew,psi_out);
00098 p=pnew;
00099 r=rnew;
00100 psi_out=solnew;
00101 step++;
00102
00103
00104
00105
00106
00107
00108
00109
00110
00111
00112 step++;
00113
00114 } while (residue>absolute_precision &&
00115 rresidue>relative_precision &&
00116 step<max_steps);
00117
00118 if(!qdaggerq) {
00119 solnew.update();
00120 mul_Q(psi_out,solnew,U,coeff);
00121 }
00122
00123 psi_out.update();
00124
00125 stats.target_absolute_precision=absolute_precision;
00126 stats.target_relative_precision=relative_precision;
00127 stats.max_steps=max_steps;
00128 stats.absolute_precision=residue;
00129 stats.relative_precision=rresidue;
00130 stats.residue=residue;
00131 stats.steps=step;
00132 stats.mul_Q_steps=2*step+1;
00133 stats.time=mpi.time()-time;
00134
00135 mpi << "\t<stats>\n"
00136 << "\t\t<max_steps>" << step << "</max_steps>\n"
00137 << "\t\t<absolute_precision>" << residue << "</absolute_precision>\n"
00138 << "\t\t<relative_precision>" << rresidue << "</relative_precision>\n"
00139 << "\t\t<time>" << stats.time << "</time>\n"
00140 << "\t</stats>\n";
00141
00142 mpi.end_function("ConjugateGradientInverter");
00143 return stats;
00144 }
00145 };