portelli/Grid
1
0
Fork 0
mirror of https://github.com/paboyle/Grid.git synced 2025-04-04 03:05:55 +01:00
Code Releases Activity

VPGCR updates

Browse Source
This commit is contained in:
Peter Boyle 2015-06-30 15:02:27 +01:00
parent 490042f8e1
commit 7cfe432ee2
1 changed files with 29 additions and 2 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

31
lib/algorithms/iterative/PrecGeneralisedConjugateResidual.h
View File

@ -5,6 +5,8 @@
//VPGCR Abe and Zhang, 2005. //VPGCR Abe and Zhang, 2005.
//INTERNATIONAL JOURNAL OF NUMERICAL ANALYSIS AND MODELING //INTERNATIONAL JOURNAL OF NUMERICAL ANALYSIS AND MODELING
//Computing and Information Volume 2, Number 2, Pages 147-161 //Computing and Information Volume 2, Number 2, Pages 147-161
//NB. Likely not original reference since they are focussing on a preconditioner variant.
// but VPGCR was nicely written up in their paper
/////////////////////////////////////////////////////////////////////////////////////////////////////// ///////////////////////////////////////////////////////////////////////////////////////////////////////
namespace Grid { namespace Grid {
@ -19,7 +21,7 @@ namespace Grid {
int steps; int steps;
LinearFunction<Field> &Preconditioner; LinearFunction<Field> &Preconditioner;
PrecGeneralisedConjugateResidual(RealD tol,Integer maxit,LinearFunction<Field> &Prec,int _mmax,int _nstep) : PrecGeneralisedConjugateResidual(RealD tol,Integer maxit,LinearFunction<Field> &Prec,int _mmax,int _nstep) :
Tolerance(tol), Tolerance(tol),
MaxIterations(maxit), MaxIterations(maxit),
Preconditioner(Prec), Preconditioner(Prec),
@ -71,6 +73,8 @@ namespace Grid {
Field r(grid); Field r(grid);
Field z(grid); Field z(grid);
Field tmp(grid);
Field ttmp(grid);
Field Az(grid); Field Az(grid);
//////////////////////////////// ////////////////////////////////
@ -91,6 +95,24 @@ namespace Grid {
// p = Prec(r) // p = Prec(r)
///////////////////// /////////////////////
Preconditioner(r,z); Preconditioner(r,z);
std::cout<< " Preconditioner in " << norm2(r)<<std::endl;
std::cout<< " Preconditioner out " << norm2(z)<<std::endl;
Linop.HermOp(z,tmp);
std::cout<< " Preconditioner Aout " << norm2(tmp)<<std::endl;
ttmp=tmp;
tmp=tmp-r;
std::cout<< " Preconditioner resid " << std::sqrt(norm2(tmp)/norm2(r))<<std::endl;
/*
std::cout<<r<<std::endl;
std::cout<<z<<std::endl;
std::cout<<ttmp<<std::endl;
std::cout<<tmp<<std::endl;
*/
Linop.HermOpAndNorm(z,Az,zAz,zAAz); Linop.HermOpAndNorm(z,Az,zAz,zAAz);
//p[0],q[0],qq[0] //p[0],q[0],qq[0]
@ -115,14 +137,19 @@ namespace Grid {
cp = axpy_norm(r,-a,q[peri_k],r); cp = axpy_norm(r,-a,q[peri_k],r);
std::cout<< " VPCG_step resid" <<sqrt(cp/rsq)<<std::endl;
if((k==nstep-1)||(cp<rsq)){ if((k==nstep-1)||(cp<rsq)){
return cp; return cp;
} }
Preconditioner(r,z);// solve Az = r Preconditioner(r,z);// solve Az = r
Linop.HermOpAndNorm(z,Az,zAz,zAAz); Linop.HermOpAndNorm(z,Az,zAz,zAAz);
Linop.HermOp(z,tmp);
tmp=tmp-r;
std::cout<< " 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;