portelli/Grid
1
0
Fork 0
mirror of https://github.com/paboyle/Grid.git synced 2025-04-09 21:50:45 +01:00
Code Releases Activity

Logging

Browse Source
This commit is contained in:
Peter Boyle 2019-12-09 02:54:52 -05:00
parent 2912071f83
commit 0dfdf80407
1 changed files with 26 additions and 12 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

38
Grid/algorithms/iterative/PrecGeneralisedConjugateResidual.h
View File

@ -86,7 +86,7 @@ public:
cp=GCRnStep(Linop,src,psi,rsq); cp=GCRnStep(Linop,src,psi,rsq);
std::cout<<GridLogMessage<<"VPGCR("<<mmax<<","<<nstep<<") "<< steps <<" steps cp = "<<cp<<std::endl; std::cout<<GridLogMessage<<"VPGCR("<<mmax<<","<<nstep<<") "<< steps <<" steps cp = "<<cp<<" target "<<rsq <<std::endl;
if(cp<rsq) { if(cp<rsq) {
@ -134,6 +134,9 @@ public:
std::vector<Field> p(mmax,grid); std::vector<Field> p(mmax,grid);
std::vector<RealD> qq(mmax); std::vector<RealD> qq(mmax);
std::cout<<GridLogIterative<< " ************** "<< std::endl;
std::cout<<GridLogIterative<< " GCRnStep("<<nstep<<")"<<std::endl;
std::cout<<GridLogIterative<< " ************** "<< std::endl;
////////////////////////////////// //////////////////////////////////
// initial guess x0 is taken as nonzero. // initial guess x0 is taken as nonzero.
@ -143,25 +146,32 @@ public:
Linop.HermOpAndNorm(psi,Az,zAz,zAAz); Linop.HermOpAndNorm(psi,Az,zAz,zAAz);
MatTimer.Stop(); MatTimer.Stop();
LinalgTimer.Start(); LinalgTimer.Start();
r=src-Az; r=src-Az;
LinalgTimer.Stop(); LinalgTimer.Stop();
std::cout<<GridLogIterative<< " GCRnStep true residual r = src - A psi "<<norm2(r) <<std::endl;
///////////////////// /////////////////////
// p = Prec(r) // p = Prec(r)
///////////////////// /////////////////////
std::cout<<GridLogIterative<< " GCRnStep apply preconditioner z= M^-1 r "<< std::endl;
std::cout<<GridLogIterative<< " --------------------------------------- "<< std::endl;
PrecTimer.Start(); PrecTimer.Start();
Preconditioner(r,z); Preconditioner(r,z);
PrecTimer.Stop(); PrecTimer.Stop();
std::cout<<GridLogIterative<< " --------------------------------------- "<< std::endl;
std::cout<<GridLogIterative<< " GCRnStep called Preconditioner z "<< norm2(z) <<std::endl;
MatTimer.Start(); // MatTimer.Start();
Linop.HermOp(z,tmp); // Linop.HermOp(z,tmp);
MatTimer.Stop(); // MatTimer.Stop();
LinalgTimer.Start(); // LinalgTimer.Start();
ttmp=tmp; // ttmp=tmp;
tmp=tmp-r; // tmp=tmp-r;
LinalgTimer.Stop(); // LinalgTimer.Stop();
/* /*
std::cout<<GridLogMessage<<r<<std::endl; std::cout<<GridLogMessage<<r<<std::endl;
@ -175,10 +185,12 @@ public:
MatTimer.Stop(); MatTimer.Stop();
LinalgTimer.Start(); LinalgTimer.Start();
//p[0],q[0],qq[0] //p[0],q[0],qq[0]
p[0]= z; p[0]= z;
q[0]= Az; q[0]= Az;
qq[0]= zAAz; qq[0]= zAAz;
std::cout<<GridLogIterative<< " GCRnStep p0=z, q0 = A p0 " <<std::endl;
cp =norm2(r); cp =norm2(r);
LinalgTimer.Stop(); LinalgTimer.Stop();
@ -200,24 +212,26 @@ public:
cp = axpy_norm(r,-a,q[peri_k],r); cp = axpy_norm(r,-a,q[peri_k],r);
LinalgTimer.Stop(); LinalgTimer.Stop();
std::cout<<GridLogMessage<< " VPGCR_step["<<steps<<"] resid " << cp << " target " <<rsq<<std::endl;
if((k==nstep-1)||(cp<rsq)){ if((k==nstep-1)||(cp<rsq)){
return cp; return cp;
} }
std::cout<<GridLogMessage<< " VPGCR_step["<<steps<<"] resid " <<sqrt(cp/rsq)<<std::endl;
std::cout<<GridLogIterative<< " GCRnStep apply preconditioner z= M^-1 r "<< std::endl;
std::cout<<GridLogIterative<< " --------------------------------------- "<< std::endl;
PrecTimer.Start(); PrecTimer.Start();
Preconditioner(r,z);// solve Az = r Preconditioner(r,z);// solve Az = r
PrecTimer.Stop(); PrecTimer.Stop();
std::cout<<GridLogIterative<< " --------------------------------------- "<< std::endl;
std::cout<<GridLogIterative<< " GCRnStep called Preconditioner z "<< norm2(z) <<std::endl;
MatTimer.Start(); MatTimer.Start();
Linop.HermOpAndNorm(z,Az,zAz,zAAz); Linop.HermOpAndNorm(z,Az,zAz,zAAz);
Linop.HermOp(z,tmp);
MatTimer.Stop(); MatTimer.Stop();
LinalgTimer.Start(); LinalgTimer.Start();
tmp=tmp-r;
std::cout<<GridLogMessage<< " Preconditioner resid " <<sqrt(norm2(tmp)/norm2(r))<<std::endl;
q[peri_kp]=Az; q[peri_kp]=Az;
p[peri_kp]=z; p[peri_kp]=z;