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

fermiqcd_staggered_uml_inverter.h

Go to the documentation of this file.
00001 
00002 
00003 
00004 
00005 
00006 
00007 
00008 
00009 
00010 
00011 
00012 
00013 // ////////////////////////////////////////////////
00014 // the MILC UML inverter
00015 // ////////////////////////////////////////////////
00016 
00045 class StaggeredBiCGUML {
00046 private:
00047   static inversion_stats staggered_BiCG_QQh(staggered_field &psi_out, 
00048                                             staggered_field &psi_in, 
00049                                             gauge_field &U,
00050                                             coefficients &coeff, 
00051                                             mdp_real mass,
00052                                             int parity=EVEN,
00053                                             mdp_real absolute_precision=staggered_inversion_precision,
00054                                             mdp_real relative_precision=0,
00055                                             int max_steps=2000) {
00056     
00057     int i, opposite_parity=EVENODD, step=0;
00058     int nc=psi_in.nc;
00059     double beta, norm, four_mass_sq;
00060     double pMMp, alpha, residue, rresidue, target_residue, old_residue;
00061     inversion_stats stats;
00062     
00063     staggered_field r(psi_in.lattice(),nc);
00064     staggered_field p(psi_in.lattice(),nc);
00065     staggered_field t(psi_in.lattice(),nc);
00066     staggered_field q(psi_in.lattice(),nc);
00067     site x(psi_in.lattice());
00068     
00069     double time=mpi.time();
00070   
00071     if(!coeff.has_key("mass") || coeff["mass"]!=0) 
00072       error("coefficient mass undeclared or different from zero");
00073     four_mass_sq=pow(2.0*mass,2);
00074     
00075     if(parity==EVEN)    opposite_parity=ODD;
00076     if(parity==ODD)     opposite_parity=EVEN;
00077     if(parity==EVENODD) opposite_parity=EVENODD;
00078 
00079     // psi_out=psi_in;
00080     // t = M^Dagger M psi_out
00081     
00082     mdp << "BEGIN BiConjugate Gradient Inversion of Q Q^DAGGER (?) ...\n";    
00083     mdp << "\tstep\tresidue\t\ttime (sec)\n";
00084     mdp << "\t====\t=======\t\t==========\n";
00085     
00086     forallsitesofparity(x,parity)
00087       for(i=0; i<nc; i++) 
00088         p(x,i)=psi_out(x,i);
00089     
00090     p.update(parity);
00091     mul_Q(q, p, U, coeff, opposite_parity); 
00092     dagger(coeff);
00093     
00094     q.update(opposite_parity);
00095     mul_Q(t, q, U, coeff, parity);
00096     dagger(coeff);
00097 
00098     forallsitesofparity(x,parity)
00099       for(i=0; i<nc; i++) {
00100         t(x,i)+=four_mass_sq*p(x,i);
00101         r(x,i)=p(x,i)=psi_in(x,i)-t(x,i);
00102       }
00103     
00104     residue=sqrt(real_scalar_product(r,r,parity));
00105     norm   =sqrt(real_scalar_product(psi_in,psi_in,parity));
00106     
00107     target_residue=absolute_precision*norm;
00108     
00109     do {
00110       p.update(parity);
00111       mul_Q(q, p, U, coeff, opposite_parity);
00112       dagger(coeff);
00113       
00114       q.update(opposite_parity);
00115       mul_Q(t, q, U, coeff, parity);
00116       dagger(coeff);
00117       
00118       if(four_mass_sq!=0)
00119         forallsitesofparity(x,parity)
00120           for(i=0; i<nc; i++)
00121             t(x,i)+=four_mass_sq*p(x,i);
00122       
00123       pMMp=real_scalar_product(p,t,parity);
00124       
00125       alpha=pow(residue,2)/pMMp; 
00126       
00127       forallsitesofparity(x,parity)
00128         for(i=0; i<nc; i++) {
00129           r(x,i)-=alpha*t(x,i);
00130           psi_out(x,i)+=alpha*p(x,i);
00131       }
00132       
00133       old_residue=residue;
00134       residue=sqrt(real_scalar_product(r,r,parity));
00135       if(residue<target_residue) break;    
00136       rresidue=relative_residue(r,psi_out,parity);
00137       if(rresidue<relative_precision) break;                                    
00138       beta=pow(residue/old_residue,2);
00139       
00140       forallsitesofparity(x,parity)
00141         for(i=0; i<nc; i++)
00142           p(x,i)=r(x,i)+beta*p(x,i);
00143       
00144       mdp << step << "\t" << residue << "\t" << mpi.time()-time << '\n';
00145 
00146       if((ME==0) && (step>100) && (residue==old_residue))
00147         error("fermiqcd_staggered_uml_inverter/staggered_BiCG_QQh: not converging"); 
00148       step++;
00149   } while(step<max_steps);
00150     mdp << "UML BiCG converged in " << step 
00151         << " iterations and " << mpi.time()-time << " seconds\n",
00152     
00153     stats.target_absolute_precision=absolute_precision;
00154     stats.target_relative_precision=relative_precision;
00155     stats.max_steps=max_steps;
00156     stats.absolute_precision=residue;
00157     stats.relative_precision=rresidue;
00158     stats.residue=residue;
00159     stats.steps=step;
00160     stats.mul_Q_steps=step+1;
00161     return stats;
00162   }
00163  public:  
00164   static inversion_stats inverter(staggered_field &psi_out, 
00165                                   staggered_field &psi_in, 
00166                                   gauge_field &U,
00167                                   coefficients &coeff, 
00168                                   mdp_real absolute_precision=staggered_inversion_precision,
00169                                   mdp_real relative_precision=0,
00170                                   int max_steps=2000) {
00171 
00172     site x(U.lattice());
00173     int i;
00174     staggered_field r(psi_in.lattice(),U.nc);
00175     mdp_real mass;
00176     inversion_stats stats;
00177     if(coeff.has_key("mass")) mass=coeff["mass"];
00178     else error("coefficient mass undeclared");
00179     coeff["mass"]=0;
00180     
00181     // set masses to zero (becasue use mul_q for d_slash)
00182     
00183     // It is important to initilize the output here
00184     // because staggered_BiCG_QQh uses it.
00185     forallsites(x)
00186     for(i=0; i<U.nc; i++) 
00187       psi_out(x,i)=0;
00188     psi_out.update();
00189     
00190     mul_Q(r,psi_in,U,coeff,EVEN);
00191     
00192     forallsitesofparity(x,EVEN) 
00193       for(i=0; i<U.nc; i++) 
00194         r(x,i)=-r(x,i)+2.0*mass*psi_in(x,i);
00195     
00196     stats=staggered_BiCG_QQh(psi_out,r,U,coeff,mass,EVEN,absolute_precision, relative_precision,max_steps);
00197     psi_out.update();
00198     mul_Q(r,psi_out,U,coeff,ODD);
00199     forallsitesofparity(x,ODD) 
00200       for(i=0; i<U.nc; i++) 
00201         psi_out(x,i)=1.0/(2.0*mass)*(psi_in(x,i)-r(x,i));
00202     
00203     // restore masses
00204     coeff["mass"]=mass;
00205     stats.mul_Q_steps++;
00206     return stats;
00207   }
00208 };

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