Main Page | Class Hierarchy | Class List | File List | Class Members | File Members

fermiqcd_bicgstab_inverter.h

Go to the documentation of this file.
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     // Initial conditions
00066     if(BiCGStabRestart==FALSE) 
00067       psi_out=psi_in; // first guess and a stupid guess!
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       // computation of residue
00106       residue=norm_square(r);
00107       residue=sqrt(residue/r.global_size());
00108 
00109       // computation of rresidue
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 }

Generated on Sun Feb 27 15:12:16 2005 by  doxygen 1.4.1