portelli/Grid
1
0
Fork 0
mirror of https://github.com/paboyle/Grid.git synced 2024-11-12 16:55:37 +00:00
Code Releases Activity

ADEF1 and ADEF2 2 level CG

Browse Source
This commit is contained in:
Peter Boyle 2023-10-05 16:55:19 -04:00
parent 59b9d0e030
commit dd557af84b
1 changed files with 129 additions and 97 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

226
Grid/algorithms/iterative/AdefGeneric.h
View File

@ -42,46 +42,30 @@ Author: Peter Boyle <paboyle@ph.ed.ac.uk>
*/ */
NAMESPACE_BEGIN(Grid); NAMESPACE_BEGIN(Grid);
template<class Field, class CoarseField, class Aggregation> template<class Field>
class TwoLevelADEF2 : public LinearFunction<Field> class TwoLevelCG : public LinearFunction<Field>
{ {
public: public:
RealD Tolerance; RealD Tolerance;
Integer MaxIterations; Integer MaxIterations;
const int mmax = 1;
GridBase *grid; GridBase *grid;
GridBase *coarsegrid;
// Fine operator, Smoother, CoarseSolver // Fine operator, Smoother, CoarseSolver
LinearOperatorBase<Field> &_FineLinop; LinearOperatorBase<Field> &_FineLinop;
LinearFunction<Field> &_Smoother; LinearFunction<Field> &_Smoother;
LinearFunction<CoarseField> &_CoarseSolver;
LinearFunction<CoarseField> &_CoarseSolverPrecise;
// Need something that knows how to get from Coarse to fine and back again
// void ProjectToSubspace(CoarseVector &CoarseVec,const FineField &FineVec){
// void PromoteFromSubspace(const CoarseVector &CoarseVec,FineField &FineVec){
Aggregation &_Aggregates;
// more most opertor functions // more most opertor functions
TwoLevelADEF2(RealD tol, TwoLevelCG(RealD tol,
Integer maxit, Integer maxit,
LinearOperatorBase<Field> &FineLinop, LinearOperatorBase<Field> &FineLinop,
LinearFunction<Field> &Smoother, LinearFunction<Field> &Smoother,
LinearFunction<CoarseField> &CoarseSolver, GridBase *fine) :
LinearFunction<CoarseField> &CoarseSolverPrecise,
Aggregation &Aggregates
) :
Tolerance(tol), Tolerance(tol),
MaxIterations(maxit), MaxIterations(maxit),
_FineLinop(FineLinop), _FineLinop(FineLinop),
_Smoother(Smoother), _Smoother(Smoother)
_CoarseSolver(CoarseSolver),
_CoarseSolverPrecise(CoarseSolverPrecise),
_Aggregates(Aggregates)
{ {
coarsegrid = Aggregates.CoarseGrid; grid = fine;
grid = Aggregates.FineGrid;
}; };
virtual void operator() (const Field &src, Field &psi) virtual void operator() (const Field &src, Field &psi)
@ -195,47 +179,8 @@ class TwoLevelADEF2 : public LinearFunction<Field>
public: public:
virtual void PcgM1(Field & in, Field & out) virtual void PcgM1(Field & in, Field & out) =0;
{ virtual void Vstart(Field & x,const Field & src)=0;
// [PTM+Q] in = [1 - Q A] M in + Q in = Min + Q [ in -A Min]
Field tmp(grid);
Field Min(grid);
CoarseField PleftProj(coarsegrid);
CoarseField PleftMss_proj(coarsegrid);
GridStopWatch SmootherTimer;
GridStopWatch MatrixTimer;
SmootherTimer.Start();
_Smoother(in,Min);
SmootherTimer.Stop();
MatrixTimer.Start();
_FineLinop.HermOp(Min,out);
MatrixTimer.Stop();
axpy(tmp,-1.0,out,in); // tmp = in - A Min
GridStopWatch ProjTimer;
GridStopWatch CoarseTimer;
GridStopWatch PromTimer;
ProjTimer.Start();
_Aggregates.ProjectToSubspace(PleftProj,tmp);
ProjTimer.Stop();
CoarseTimer.Start();
_CoarseSolver(PleftProj,PleftMss_proj); // Ass^{-1} [in - A Min]_s
CoarseTimer.Stop();
PromTimer.Start();
_Aggregates.PromoteFromSubspace(PleftMss_proj,tmp);// tmp = Q[in - A Min]
PromTimer.Stop();
std::cout << GridLogMessage << "PcgM1 breakdown "<<std::endl;
std::cout << GridLogMessage << "\tSmoother " << SmootherTimer.Elapsed() <<std::endl;
std::cout << GridLogMessage << "\tMatrix " << MatrixTimer.Elapsed() <<std::endl;
std::cout << GridLogMessage << "\tProj " << ProjTimer.Elapsed() <<std::endl;
std::cout << GridLogMessage << "\tCoarse " << CoarseTimer.Elapsed() <<std::endl;
std::cout << GridLogMessage << "\tProm " << PromTimer.Elapsed() <<std::endl;
axpy(out,1.0,Min,tmp); // Min+tmp
}
virtual void PcgM2(const Field & in, Field & out) { virtual void PcgM2(const Field & in, Field & out) {
out=in; out=in;
@ -249,6 +194,89 @@ class TwoLevelADEF2 : public LinearFunction<Field>
return dd; return dd;
} }
/////////////////////////////////////////////////////////////////////
// Only Def1 has non-trivial Vout.
/////////////////////////////////////////////////////////////////////
virtual void Vout (Field & in, Field & out,Field & src){
out = in;
}
};
template<class Field, class CoarseField, class Aggregation>
class TwoLevelADEF2 : public TwoLevelCG<Field>
{
public:
///////////////////////////////////////////////////////////////////////////////////
// Need something that knows how to get from Coarse to fine and back again
// void ProjectToSubspace(CoarseVector &CoarseVec,const FineField &FineVec){
// void PromoteFromSubspace(const CoarseVector &CoarseVec,FineField &FineVec){
///////////////////////////////////////////////////////////////////////////////////
GridBase *coarsegrid;
Aggregation &_Aggregates;
LinearFunction<CoarseField> &_CoarseSolver;
LinearFunction<CoarseField> &_CoarseSolverPrecise;
///////////////////////////////////////////////////////////////////////////////////
// more most opertor functions
TwoLevelADEF2(RealD tol,
Integer maxit,
LinearOperatorBase<Field> &FineLinop,
LinearFunction<Field> &Smoother,
LinearFunction<CoarseField> &CoarseSolver,
LinearFunction<CoarseField> &CoarseSolverPrecise,
Aggregation &Aggregates
) :
TwoLevelCG<Field>(tol,maxit,FineLinop,Smoother,Aggregates.FineGrid),
_CoarseSolver(CoarseSolver),
_CoarseSolverPrecise(CoarseSolverPrecise),
_Aggregates(Aggregates)
{
coarsegrid = Aggregates.CoarseGrid;
};
virtual void PcgM1(Field & in, Field & out)
{
// [PTM+Q] in = [1 - Q A] M in + Q in = Min + Q [ in -A Min]
Field tmp(this->grid);
Field Min(this->grid);
CoarseField PleftProj(this->coarsegrid);
CoarseField PleftMss_proj(this->coarsegrid);
GridStopWatch SmootherTimer;
GridStopWatch MatrixTimer;
SmootherTimer.Start();
this->_Smoother(in,Min);
SmootherTimer.Stop();
MatrixTimer.Start();
this->_FineLinop.HermOp(Min,out);
MatrixTimer.Stop();
axpy(tmp,-1.0,out,in); // tmp = in - A Min
GridStopWatch ProjTimer;
GridStopWatch CoarseTimer;
GridStopWatch PromTimer;
ProjTimer.Start();
this->_Aggregates.ProjectToSubspace(PleftProj,tmp);
ProjTimer.Stop();
CoarseTimer.Start();
this->_CoarseSolver(PleftProj,PleftMss_proj); // Ass^{-1} [in - A Min]_s
CoarseTimer.Stop();
PromTimer.Start();
this->_Aggregates.PromoteFromSubspace(PleftMss_proj,tmp);// tmp = Q[in - A Min]
PromTimer.Stop();
std::cout << GridLogPerformance << "PcgM1 breakdown "<<std::endl;
std::cout << GridLogPerformance << "\tSmoother " << SmootherTimer.Elapsed() <<std::endl;
std::cout << GridLogPerformance << "\tMatrix " << MatrixTimer.Elapsed() <<std::endl;
std::cout << GridLogPerformance << "\tProj " << ProjTimer.Elapsed() <<std::endl;
std::cout << GridLogPerformance << "\tCoarse " << CoarseTimer.Elapsed() <<std::endl;
std::cout << GridLogPerformance << "\tProm " << PromTimer.Elapsed() <<std::endl;
axpy(out,1.0,Min,tmp); // Min+tmp
}
virtual void Vstart(Field & x,const Field & src) virtual void Vstart(Field & x,const Field & src)
{ {
/////////////////////////////////// ///////////////////////////////////
@ -262,65 +290,69 @@ class TwoLevelADEF2 : public LinearFunction<Field>
// = src_s - (A guess)_s - src_s + (A guess)_s // = src_s - (A guess)_s - src_s + (A guess)_s
// = 0 // = 0
/////////////////////////////////// ///////////////////////////////////
Field r(grid); Field r(this->grid);
Field mmp(grid); Field mmp(this->grid);
CoarseField PleftProj(coarsegrid); CoarseField PleftProj(this->coarsegrid);
CoarseField PleftMss_proj(coarsegrid); CoarseField PleftMss_proj(this->coarsegrid);
_Aggregates.ProjectToSubspace(PleftProj,src); this->_Aggregates.ProjectToSubspace(PleftProj,src);
_CoarseSolverPrecise(PleftProj,PleftMss_proj); // Ass^{-1} r_s this->_CoarseSolverPrecise(PleftProj,PleftMss_proj); // Ass^{-1} r_s
_Aggregates.PromoteFromSubspace(PleftMss_proj,x); this->_Aggregates.PromoteFromSubspace(PleftMss_proj,x);
} }
/////////////////////////////////////////////////////////////////////
// Only Def1 has non-trivial Vout.
/////////////////////////////////////////////////////////////////////
virtual void Vout (Field & in, Field & out,Field & src){
out = in;
}
}; };
template<class Field, class CoarseField, class Aggregation> template<class Field>
class TwoLevelADEF1ev : public TwoLevelADEF2<Field,CoarseField,Aggregation> class TwoLevelADEF1defl : public TwoLevelCG<Field>
{ {
TwoLevelADEF1ev(RealD tol, public:
const std::vector<Field> &evec;
const std::vector<RealD> &eval;
TwoLevelADEF1defl(RealD tol,
Integer maxit, Integer maxit,
LinearOperatorBase<Field> &FineLinop, LinearOperatorBase<Field> &FineLinop,
LinearFunction<Field> &Smoother, LinearFunction<Field> &Smoother,
LinearFunction<CoarseField> &CoarseSolver, std::vector<Field> &_evec,
LinearFunction<CoarseField> &CoarseSolverPrecise, std::vector<RealD> &_eval) :
Aggregation &Aggregates ) : TwoLevelCG<Field>(tol,maxit,FineLinop,Smoother,_evec[0].Grid()),
TwoLevelADEF2<Field,CoarseField,Aggregation>(tol,maxit,FineLinop,Smoother,CoarseSolver,CoarseSolverPrecise,Aggregates) evec(_evec),
eval(_eval)
{}; {};
// Can just inherit existing Vstart
// Can just inherit existing Vout // Can just inherit existing Vout
// Can just inherit existing M2 // Can just inherit existing M2
// Can just inherit existing M3 // Can just inherit existing M3
// Simple vstart - do nothing
virtual void Vstart(Field & x,const Field & src){};
// Override PcgM1 // Override PcgM1
virtual void PcgM1(Field & in, Field & out) virtual void PcgM1(Field & in, Field & out)
{ {
CoarseField PleftProj(this->coarsegrid); int N=evec.size();
CoarseField PleftMss_proj(this->coarsegrid);
Field tmp(this->grid);
Field Pin(this->grid); Field Pin(this->grid);
Field Qin(this->grid);
//MP + Q = M(1-AQ) + Q = M //MP + Q = M(1-AQ) + Q = M
// // If we are eigenvector deflating in coarse space // // If we are eigenvector deflating in coarse space
// // Q = Sum_i |phi_i> 1/lambda_i <phi_i| // // Q = Sum_i |phi_i> 1/lambda_i <phi_i|
// // A Q = Sum_i |phi_i> <phi_i| // // A Q = Sum_i |phi_i> <phi_i|
// // M(1-AQ) = M(1-proj) + Q // // M(1-AQ) = M(1-proj) + Q
this->_Aggregates.ProjectToSubspace(PleftProj,in); Qin.Checkerboard()=in.Checkerboard();
this->_Aggregates.PromoteFromSubspace(PleftProj,tmp);// tmp = Qin Qin = Zero();
Pin = in;
for (int i=0;i<N;i++) {
const Field& tmp = evec[i];
auto ip = TensorRemove(innerProduct(tmp,in));
axpy(Qin, ip / eval[i],tmp,Qin);
axpy(Pin, -ip ,tmp,Pin);
}
Pin = in - tmp;
this->_Smoother(Pin,out); this->_Smoother(Pin,out);
this->_CoarseSolver(PleftProj,PleftMss_proj); out = out + Qin;
this->_Aggregates.PromoteFromSubspace(PleftMss_proj,tmp);// tmp = Qin
out = out + tmp;
} }
}; };