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

fermiqcd_minres_inverter.h

Go to the documentation of this file.
00001 
00002 
00003 
00004 
00005 
00006 
00007 
00008 
00009 
00010 
00011 
00012 
00016 class inversion_stats {
00017  public:
00018   mdp_real target_absolute_precision;
00019   mdp_real target_relative_precision;
00020   int      max_steps;
00021   mdp_real residue;
00022   mdp_real absolute_precision;
00023   mdp_real relative_precision;
00024   int      steps;
00025   int      mul_Q_steps;
00026   mdp_real time;
00027 };
00028 
00029 // /////////////////////////////////////////////
00030 // implementation of the minimum residue inversion
00031 // /////////////////////////////////////////////
00032 
00060 class MinRes {
00061  public:
00062   template <class fieldT, class fieldG>
00063   static inversion_stats inverter(fieldT &psi_out, 
00064                                   fieldT &psi_in, 
00065                                   fieldG &U, 
00066                                   coefficients &coeff, 
00067                                   mdp_real absolute_precision=mdp_precision,
00068                                   mdp_real relative_precision=0,
00069                                   int max_steps=2000) {
00070     
00071     mpi.begin_function("MinimumResidueInverter");
00072     
00073     fieldT          q(psi_in);
00074     fieldT          r(psi_in);
00075     double          residue, rresidue=-1, old_rresidue;
00076     mdp_complex         alpha;
00077     long            step=0;
00078     double          time=mpi.time();
00079     inversion_stats stats;
00080     
00081     psi_in.update();
00082     mul_Q(r,psi_in,U,coeff);
00083     r*=-1;
00084     r+=psi_in;
00085     psi_out=psi_in;
00086     
00087     mpi << "\t<target>\n" 
00088         << "\t\t<max_steps>" << max_steps << "</max_steps>\n"
00089         << "\t\t<absolute_precision>" << absolute_precision << "</absolute_precision>\n"
00090         << "\t\t<relative_precision>" << relative_precision << "</relative_precision>\n"
00091         << "\t</target>\n";
00092     
00093     do {
00094       r.update();
00095       mul_Q(q,r,U,coeff);
00096       alpha=q*r;
00097       alpha/=norm_square(q);
00098       mdp_add_scaled_field(psi_out, alpha, r);
00099       mdp_add_scaled_field(r, -alpha, q);
00100       
00101       // computing residue
00102       residue=norm_square(r);
00103       residue=sqrt(residue/r.global_size());
00104       
00105       // computing relative residue
00106       old_rresidue=rresidue;
00107       rresidue=relative_residue(r,psi_out);
00108       
00109       mpi << "\t\t<step>" << step << "</step>\n"
00110           << "\t\t<residue>" << residue << "</residue>\n"
00111           << "\t\t<relative_residue>" << rresidue << "</relative_residue>\n"
00112           << "\t\t<time>" << mpi.time()-time << "</time>\n\n";
00113       
00114       if((step>10) && (rresidue==old_rresidue))
00115         error("not converging"); 
00116       step++;
00117     } while (residue>absolute_precision &&
00118              rresidue>relative_precision &&
00119              step<max_steps);
00120     
00121     psi_out.update();
00122     
00123     stats.target_absolute_precision=absolute_precision;
00124     stats.target_relative_precision=relative_precision;
00125     stats.max_steps=max_steps;
00126     stats.absolute_precision=residue;
00127     stats.relative_precision=rresidue;
00128     stats.residue=residue;
00129     stats.steps=step;
00130     stats.mul_Q_steps=step+1;
00131     stats.time=mpi.time()-time;
00132     
00133     mpi << "\t<stats>\n" 
00134         << "\t\t<max_steps>" << step << "</max_steps>\n"
00135         << "\t\t<absolute_precision>" << residue << "</absolute_precision>\n"
00136         << "\t\t<relative_precision>" << rresidue << "</relative_precision>\n"
00137         << "\t\t<time>" << stats.time << "</time>\n"
00138         << "\t</stats>\n";
00139     
00140     mpi.end_function("MinimumResidueInverter");
00141     return stats;
00142   }
00143 };
00144 
00145 template <class fieldT, class fieldG>
00146 inversion_stats MinimumResidueInverter(fieldT &psi_out, 
00147                                        fieldT &psi_in, 
00148                                        fieldG &U, 
00149                                        coefficients &coeff, 
00150                                        mdp_real absolute_precision=mdp_precision,
00151                                        mdp_real relative_precision=0,
00152                                        int max_steps=2000) {
00153   return MinRes::inverter(psi_out,psi_in,U,coeff,
00154                           absolute_precision,
00155                           relative_precision,
00156                           max_steps);
00157 }

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