00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00041 class BiCGStab {
00042 public:
00043 template <class fieldT, class fieldG>
00044 static inversion_stats inverter(fieldT &psi_out,
00045 fieldT &psi_in,
00046 fieldG &U,
00047 coefficients &coeff,
00048 mdp_real absolute_precision=mdp_precision,
00049 mdp_real relative_precision=0,
00050 int max_steps=2000) {
00051 mpi.begin_function("BiConugateGradientStabilizedInverter");
00052 int step=0;
00053 fieldT p(psi_in);
00054 fieldT q(psi_in);
00055 fieldT r(psi_in);
00056 fieldT s(psi_in);
00057 fieldT t(psi_in);
00058 double residue, rresidue=-1, old_rresidue;
00059 mdp_complex alpha, beta, rho, rho_old, omega;
00060 double time=mpi.time();
00061 inversion_stats stats;
00062
00063 mpi << "\tstep\tresidue\t\ttime (sec)\n";
00064
00065
00066 if(BiCGStabRestart==false)
00067 psi_out=psi_in;
00068
00069 psi_out.update();
00070 mul_Q(r,psi_out,U,coeff);
00071 r*=-1;
00072 r+=psi_in;
00073 q=r;
00074
00075 p=0.0;
00076 s=0.0;
00077
00078 rho_old=alpha=omega=1;
00079
00080 mpi << "\t<target>\n"
00081 << "\t\t<max_steps>" << max_steps << "</max_steps>\n"
00082 << "\t\t<absolute_precision>" << absolute_precision << "</absolute_precision>\n"
00083 << "\t\t<relative_precision>" << relative_precision << "</relative_precision>\n"
00084 << "\t</target>\n";
00085
00086 do {
00087
00088 rho=q*r;
00089 beta=(rho/rho_old)*(alpha/omega);
00090 rho_old=rho;
00091 p*=beta;
00092 p+=r;
00093 mdp_add_scaled_field(p, -beta*omega, s);
00094 p.update();
00095 mul_Q(s,p,U,coeff);
00096 alpha=rho/(q*s);
00097 mdp_add_scaled_field(r, -alpha, s);
00098 r.update();
00099 mul_Q(t,r,U,coeff);
00100 omega=t*r;
00101 omega/=norm_square(t);
00102 mdp_add_scaled_field(psi_out, omega, r);
00103 mdp_add_scaled_field(psi_out, alpha, p);
00104
00105
00106 residue=norm_square(r);
00107 residue=sqrt(residue/r.global_size());
00108
00109
00110 old_rresidue=rresidue;
00111 rresidue=relative_residue(r,psi_out);
00112
00113 mdp_add_scaled_field(r, -omega, t);
00114
00115 mpi << "\t\t<step>" << step << "</step>\n"
00116 << "\t\t<residue>" << residue << "</residue>\n"
00117 << "\t\t<relative_residue>" << rresidue << "</relative_residue>\n"
00118 << "\t\t<time>" << mpi.time()-time << "</time>\n\n";
00119
00120 if((step>10) && (rresidue==old_rresidue))
00121 error("not converging");
00122 step++;
00123
00124 } while (residue>absolute_precision &&
00125 rresidue>relative_precision &&
00126 step<max_steps);
00127
00128 psi_out.update();
00129
00130 stats.target_absolute_precision=absolute_precision;
00131 stats.target_relative_precision=relative_precision;
00132 stats.max_steps=max_steps;
00133 stats.absolute_precision=residue;
00134 stats.relative_precision=rresidue;
00135 stats.residue=residue;
00136 stats.steps=step;
00137 stats.mul_Q_steps=2*step+1;
00138 stats.time=mpi.time()-time;
00139
00140 mpi << "\t<stats>\n"
00141 << "\t\t<max_steps>" << step << "</max_steps>\n"
00142 << "\t\t<absolute_precision>" << residue << "</absolute_precision>\n"
00143 << "\t\t<relative_precision>" << rresidue << "</relative_precision>\n"
00144 << "\t\t<time>" << stats.time << "</time>\n"
00145 << "\t</stats>\n";
00146
00147 mpi.end_function("BiConugateGradientStabilizedInverter");
00148 return stats;
00149 }
00150 };
00151
00152 template <class fieldT, class fieldG>
00153 inversion_stats BiConjugateGradientStabilizedInverter(fieldT &psi_out,
00154 fieldT &psi_in,
00155 fieldG &U,
00156 coefficients &coeff,
00157 mdp_real absolute_precision=mdp_precision,
00158 mdp_real relative_precision=0,
00159 int max_steps=2000) {
00160 return BiCGStab::inverter(psi_out,psi_in,U,coeff,
00161 absolute_precision,
00162 relative_precision,
00163 max_steps);
00164 }
