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

Hide internal data

Browse Source
This commit is contained in:
paboyle 2018-01-26 23:07:34 +00:00
parent 8dccffdfd5
commit 1ebd56c3fb
10 changed files with 86 additions and 86 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

18
lib/algorithms/CoarsenedMatrix.h
View File

@ -129,7 +129,7 @@ public:
blockProject(iProj,subspace[i],subspace); blockProject(iProj,subspace[i],subspace);
eProj=zero; eProj=zero;
parallel_for(int ss=0;ss<CoarseGrid->oSites();ss++){ parallel_for(int ss=0;ss<CoarseGrid->oSites();ss++){
eProj._odata[ss](i)=CComplex(1.0); eProj[ss](i)=CComplex(1.0);
} }
eProj=eProj - iProj; eProj=eProj - iProj;
std::cout<<GridLogMessage<<"Orthog check error "<<i<<" " << norm2(eProj)<<std::endl; std::cout<<GridLogMessage<<"Orthog check error "<<i<<" " << norm2(eProj)<<std::endl;
@ -140,7 +140,7 @@ public:
blockProject(CoarseVec,FineVec,subspace); blockProject(CoarseVec,FineVec,subspace);
} }
void PromoteFromSubspace(const CoarseVector &CoarseVec,FineField &FineVec){ void PromoteFromSubspace(const CoarseVector &CoarseVec,FineField &FineVec){
FineVec.checkerboard = subspace[0].checkerboard; FineVec.Checkerboard() = subspace[0].Checkerboard();
blockPromote(CoarseVec,FineVec,subspace); blockPromote(CoarseVec,FineVec,subspace);
} }
void CreateSubspaceRandom(GridParallelRNG &RNG){ void CreateSubspaceRandom(GridParallelRNG &RNG){
@ -282,15 +282,15 @@ public:
SE=Stencil.GetEntry(ptype,point,ss); SE=Stencil.GetEntry(ptype,point,ss);
if(SE->_is_local&&SE->_permute) { if(SE->_is_local&&SE->_permute) {
permute(nbr,in._odata[SE->_offset],ptype); permute(nbr,in[SE->_offset],ptype);
} else if(SE->_is_local) { } else if(SE->_is_local) {
nbr = in._odata[SE->_offset]; nbr = in[SE->_offset];
} else { } else {
nbr = Stencil.CommBuf()[SE->_offset]; nbr = Stencil.CommBuf()[SE->_offset];
} }
res = res + A[point]._odata[ss]*nbr; res = res + A[point][ss]*nbr;
} }
vstream(out._odata[ss],res); vstream(out[ss],res);
} }
return norm2(out); return norm2(out);
}; };
@ -387,9 +387,9 @@ public:
parallel_for(int ss=0;ss<Grid()->oSites();ss++){ parallel_for(int ss=0;ss<Grid()->oSites();ss++){
for(int j=0;j<nbasis;j++){ for(int j=0;j<nbasis;j++){
if( disp!= 0 ) { if( disp!= 0 ) {
A[p]._odata[ss](j,i) = oProj._odata[ss](j); A[p][ss](j,i) = oProj[ss](j);
} }
A[self_stencil]._odata[ss](j,i) = A[self_stencil]._odata[ss](j,i) + iProj._odata[ss](j); A[self_stencil][ss](j,i) = A[self_stencil][ss](j,i) + iProj[ss](j);
} }
} }
} }
@ -443,7 +443,7 @@ public:
A[8] = val*ident; A[8] = val*ident;
// for(int s=0;s<Grid()->oSites();s++) { // for(int s=0;s<Grid()->oSites();s++) {
// A[8]._odata[s]=val._odata[s]; // A[8][s]=val[s];
// } // }
} }
void ForceHermitian(void) { void ForceHermitian(void) {

8
lib/algorithms/FFT.h
View File

@ -217,8 +217,8 @@ public:
FFTW_plan p; FFTW_plan p;
{ {
FFTW_scalar *in = (FFTW_scalar *)&pgbuf._odata[0]; FFTW_scalar *in = (FFTW_scalar *)&pgbuf[0];
FFTW_scalar *out= (FFTW_scalar *)&pgbuf._odata[0]; FFTW_scalar *out= (FFTW_scalar *)&pgbuf[0];
p = FFTW<scalar>::fftw_plan_many_dft(rank,n,howmany, p = FFTW<scalar>::fftw_plan_many_dft(rank,n,howmany,
in,inembed, in,inembed,
istride,idist, istride,idist,
@ -261,8 +261,8 @@ public:
thread_loop_in_region( (int idx=0;idx<NN;idx++), { thread_loop_in_region( (int idx=0;idx<NN;idx++), {
pencil_g.LocalIndexToLocalCoor(idx, cbuf); pencil_g.LocalIndexToLocalCoor(idx, cbuf);
if ( cbuf[dim] == 0 ) { // restricts loop to plane at lcoor[dim]==0 if ( cbuf[dim] == 0 ) { // restricts loop to plane at lcoor[dim]==0
FFTW_scalar *in = (FFTW_scalar *)&pgbuf._odata[idx]; FFTW_scalar *in = (FFTW_scalar *)&pgbuf[idx];
FFTW_scalar *out= (FFTW_scalar *)&pgbuf._odata[idx]; FFTW_scalar *out= (FFTW_scalar *)&pgbuf[idx];
FFTW<scalar>::fftw_execute_dft(p,in,out); FFTW<scalar>::fftw_execute_dft(p,in,out);
} }
}); });

2
lib/algorithms/iterative/AdefGeneric.h
View File

@ -90,7 +90,7 @@ class TwoLevelFlexiblePcg : public LinearFunction<Field>
void operator() (const Field &src, Field &psi){ void operator() (const Field &src, Field &psi){
void operator() (const Field &src, Field &psi){ void operator() (const Field &src, Field &psi){
psi.checkerboard = src.checkerboard; psi.Checkerboard() = src.Checkerboard();
grid = src._grid; grid = src._grid;
RealD f; RealD f;

6
lib/algorithms/iterative/BlockConjugateGradient.h
View File

@ -143,7 +143,7 @@ public:
std::cout<<GridLogMessage<<" Block Conjugate Gradient : Orthog "<<Orthog<<" Nblock "<<Nblock<<std::endl; std::cout<<GridLogMessage<<" Block Conjugate Gradient : Orthog "<<Orthog<<" Nblock "<<Nblock<<std::endl;
X.checkerboard = B.checkerboard; X.Checkerboard() = B.Checkerboard();
conformable(X, B); conformable(X, B);
Field tmp(B); Field tmp(B);
@ -327,7 +327,7 @@ public:
std::cout<<GridLogMessage<<" Block Conjugate Gradient : Orthog "<<Orthog<<" Nblock "<<Nblock<<std::endl; std::cout<<GridLogMessage<<" Block Conjugate Gradient : Orthog "<<Orthog<<" Nblock "<<Nblock<<std::endl;
Psi.checkerboard = Src.checkerboard; Psi.Checkerboard() = Src.Checkerboard();
conformable(Psi, Src); conformable(Psi, Src);
Field P(Src); Field P(Src);
@ -474,7 +474,7 @@ public:
std::cout<<GridLogMessage<<"MultiRHS Conjugate Gradient : Orthog "<<Orthog<<" Nblock "<<Nblock<<std::endl; std::cout<<GridLogMessage<<"MultiRHS Conjugate Gradient : Orthog "<<Orthog<<" Nblock "<<Nblock<<std::endl;
Psi.checkerboard = Src.checkerboard; Psi.Checkerboard() = Src.Checkerboard();
conformable(Psi, Src); conformable(Psi, Src);
Field P(Src); Field P(Src);

2
lib/algorithms/iterative/ConjugateGradient.h
View File

@ -54,7 +54,7 @@ public:
void operator()(LinearOperatorBase<Field> &Linop, const Field &src, Field &psi) { void operator()(LinearOperatorBase<Field> &Linop, const Field &src, Field &psi) {
psi.checkerboard = src.checkerboard; psi.Checkerboard() = src.Checkerboard();
conformable(psi, src); conformable(psi, src);
RealD cp, c, a, d, b, ssq, qq, b_pred; RealD cp, c, a, d, b, ssq, qq, b_pred;

12
lib/algorithms/iterative/ConjugateGradientMixedPrec.h
View File

@ -67,18 +67,18 @@ public:
GridStopWatch TotalTimer; GridStopWatch TotalTimer;
TotalTimer.Start(); TotalTimer.Start();
int cb = src_d_in.checkerboard; int cb = src_d_in.Checkerboard();
sol_d.checkerboard = cb; sol_d.Checkerboard() = cb;
RealD src_norm = norm2(src_d_in); RealD src_norm = norm2(src_d_in);
RealD stop = src_norm * Tolerance*Tolerance; RealD stop = src_norm * Tolerance*Tolerance;
GridBase* DoublePrecGrid = src_d_in._grid; GridBase* DoublePrecGrid = src_d_in._grid;
FieldD tmp_d(DoublePrecGrid); FieldD tmp_d(DoublePrecGrid);
tmp_d.checkerboard = cb; tmp_d.Checkerboard() = cb;
FieldD tmp2_d(DoublePrecGrid); FieldD tmp2_d(DoublePrecGrid);
tmp2_d.checkerboard = cb; tmp2_d.Checkerboard() = cb;
FieldD src_d(DoublePrecGrid); FieldD src_d(DoublePrecGrid);
src_d = src_d_in; //source for next inner iteration, computed from residual during operation src_d = src_d_in; //source for next inner iteration, computed from residual during operation
@ -86,10 +86,10 @@ public:
RealD inner_tol = InnerTolerance; RealD inner_tol = InnerTolerance;
FieldF src_f(SinglePrecGrid); FieldF src_f(SinglePrecGrid);
src_f.checkerboard = cb; src_f.Checkerboard() = cb;
FieldF sol_f(SinglePrecGrid); FieldF sol_f(SinglePrecGrid);
sol_f.checkerboard = cb; sol_f.Checkerboard() = cb;
ConjugateGradient<FieldF> CG_f(inner_tol, MaxInnerIterations); ConjugateGradient<FieldF> CG_f(inner_tol, MaxInnerIterations);
CG_f.ErrorOnNoConverge = false; CG_f.ErrorOnNoConverge = false;

4
lib/algorithms/iterative/ConjugateGradientReliableUpdate.h
View File

@ -76,7 +76,7 @@ public:
LinearOperatorBase<FieldF> *Linop_f_use = &Linop_f; LinearOperatorBase<FieldF> *Linop_f_use = &Linop_f;
bool using_fallback = false; bool using_fallback = false;
psi.checkerboard = src.checkerboard; psi.Checkerboard() = src.Checkerboard();
conformable(psi, src); conformable(psi, src);
RealD cp, c, a, d, b, ssq, qq, b_pred; RealD cp, c, a, d, b, ssq, qq, b_pred;
@ -116,7 +116,7 @@ public:
//Single prec initialization //Single prec initialization
FieldF r_f(SinglePrecGrid); FieldF r_f(SinglePrecGrid);
r_f.checkerboard = r.checkerboard; r_f.Checkerboard() = r.Checkerboard();
precisionChange(r_f, r); precisionChange(r_f, r);
FieldF psi_f(r_f); FieldF psi_f(r_f);

14
lib/algorithms/iterative/ImplicitlyRestartedLanczos.h
View File

@ -64,11 +64,11 @@ void basisRotate(std::vector<Field> &basis,Eigen::MatrixXd& Qt,int j0, int j1, i
for(int j=j0; j<j1; ++j){ for(int j=j0; j<j1; ++j){
for(int k=k0; k<k1; ++k){ for(int k=k0; k<k1; ++k){
B[j] +=Qt(j,k) * basis[k]._odata[ss]; B[j] +=Qt(j,k) * basis[k][ss];
} }
} }
for(int j=j0; j<j1; ++j){ for(int j=j0; j<j1; ++j){
basis[j]._odata[ss] = B[j]; basis[j][ss] = B[j];
} }
} }
} }
@ -81,13 +81,13 @@ void basisRotateJ(Field &result,std::vector<Field> &basis,Eigen::MatrixXd& Qt,in
typedef typename Field::vector_object vobj; typedef typename Field::vector_object vobj;
GridBase* grid = basis[0]._grid; GridBase* grid = basis[0]._grid;
result.checkerboard = basis[0].checkerboard; result.Checkerboard() = basis[0].Checkerboard();
parallel_for(int ss=0;ss < grid->oSites();ss++){ parallel_for(int ss=0;ss < grid->oSites();ss++){
vobj B = zero; vobj B = zero;
for(int k=k0; k<k1; ++k){ for(int k=k0; k<k1; ++k){
B +=Qt(j,k) * basis[k]._odata[ss]; B +=Qt(j,k) * basis[k][ss];
} }
result._odata[ss] = B; result[ss] = B;
} }
} }
@ -487,7 +487,7 @@ until convergence
std::cout << GridLogIRL << "Test convergence: rotate subset of vectors to test convergence " << std::endl; std::cout << GridLogIRL << "Test convergence: rotate subset of vectors to test convergence " << std::endl;
Field B(grid); B.checkerboard = evec[0].checkerboard; Field B(grid); B.Checkerboard() = evec[0].Checkerboard();
// power of two search pattern; not every evalue in eval2 is assessed. // power of two search pattern; not every evalue in eval2 is assessed.
for(int jj = 1; jj<=Nstop; jj*=2){ for(int jj = 1; jj<=Nstop; jj*=2){
@ -525,7 +525,7 @@ until convergence
converged: converged:
{ {
Field B(grid); B.checkerboard = evec[0].checkerboard; Field B(grid); B.Checkerboard() = evec[0].Checkerboard();
basisRotate(evec,Qt,0,Nk,0,Nk,Nm); basisRotate(evec,Qt,0,Nk,0,Nk,Nm);
std::cout << GridLogIRL << " Rotated basis"<<std::endl; std::cout << GridLogIRL << " Rotated basis"<<std::endl;
Nconv=0; Nconv=0;

16
lib/algorithms/iterative/LocalCoherenceLanczos.h
View File

@ -113,8 +113,8 @@ public:
GridBase *FineGrid = _Aggregate.FineGrid; GridBase *FineGrid = _Aggregate.FineGrid;
FineField fin(FineGrid) ;fin.checkerboard =_Aggregate.checkerboard; FineField fin(FineGrid) ;fin.Checkerboard() =_Aggregate.Checkerboard();
FineField fout(FineGrid);fout.checkerboard =_Aggregate.checkerboard; FineField fout(FineGrid);fout.Checkerboard() =_Aggregate.Checkerboard();
_Aggregate.PromoteFromSubspace(in,fin); std::cout<<GridLogIRL<<"ProjectedFunctionHermop : Promote to fine"<<std::endl; _Aggregate.PromoteFromSubspace(in,fin); std::cout<<GridLogIRL<<"ProjectedFunctionHermop : Promote to fine"<<std::endl;
_poly(_Linop,fin,fout); std::cout<<GridLogIRL<<"ProjectedFunctionHermop : Poly "<<std::endl; _poly(_Linop,fin,fout); std::cout<<GridLogIRL<<"ProjectedFunctionHermop : Poly "<<std::endl;
@ -172,10 +172,10 @@ public:
{ {
GridBase *FineGrid = _Aggregate.FineGrid; GridBase *FineGrid = _Aggregate.FineGrid;
int checkerboard = _Aggregate.checkerboard; int checkerboard = _Aggregate.Checkerboard();
FineField fB(FineGrid);fB.checkerboard =checkerboard; FineField fB(FineGrid);fB.Checkerboard() =checkerboard;
FineField fv(FineGrid);fv.checkerboard =checkerboard; FineField fv(FineGrid);fv.Checkerboard() =checkerboard;
_Aggregate.PromoteFromSubspace(B,fv); _Aggregate.PromoteFromSubspace(B,fv);
_smoother(_Linop,fv,fB); _smoother(_Linop,fv,fB);
@ -254,11 +254,11 @@ public:
int Nk = nbasis; int Nk = nbasis;
_Aggregate.subspace.resize(Nk,_FineGrid); _Aggregate.subspace.resize(Nk,_FineGrid);
_Aggregate.subspace[0]=1.0; _Aggregate.subspace[0]=1.0;
_Aggregate.subspace[0].checkerboard=_checkerboard; _Aggregate.subspace[0].Checkerboard()=_checkerboard;
normalise(_Aggregate.subspace[0]); normalise(_Aggregate.subspace[0]);
PlainHermOp<FineField> Op(_FineOp); PlainHermOp<FineField> Op(_FineOp);
for(int k=1;k<Nk;k++){ for(int k=1;k<Nk;k++){
_Aggregate.subspace[k].checkerboard=_checkerboard; _Aggregate.subspace[k].Checkerboard()=_checkerboard;
Op(_Aggregate.subspace[k-1],_Aggregate.subspace[k]); Op(_Aggregate.subspace[k-1],_Aggregate.subspace[k]);
normalise(_Aggregate.subspace[k]); normalise(_Aggregate.subspace[k]);
} }
@ -308,7 +308,7 @@ public:
ImplicitlyRestartedLanczos<FineField> IRL(ChebyOp,Op,Nstop,Nk,Nm,resid,MaxIt,betastp,MinRes); ImplicitlyRestartedLanczos<FineField> IRL(ChebyOp,Op,Nstop,Nk,Nm,resid,MaxIt,betastp,MinRes);
FineField src(_FineGrid); src=1.0; src.checkerboard = _checkerboard; FineField src(_FineGrid); src=1.0; src.Checkerboard() = _checkerboard;
int Nconv; int Nconv;
IRL.calc(evals_fine,_Aggregate.subspace,src,Nconv,false); IRL.calc(evals_fine,_Aggregate.subspace,src,Nconv,false);

90
lib/algorithms/iterative/SchurRedBlack.h
View File

@ -135,30 +135,30 @@ public:
///////////////////////////////////////////////////// /////////////////////////////////////////////////////
// src_o = (source_o - Moe MeeInv source_e) // src_o = (source_o - Moe MeeInv source_e)
///////////////////////////////////////////////////// /////////////////////////////////////////////////////
_Matrix.MooeeInv(src_e,tmp); assert( tmp.checkerboard ==Even); _Matrix.MooeeInv(src_e,tmp); assert( tmp.Checkerboard() ==Even);
_Matrix.Meooe (tmp,Mtmp); assert( Mtmp.checkerboard ==Odd); _Matrix.Meooe (tmp,Mtmp); assert( Mtmp.Checkerboard() ==Odd);
tmp=src_o-Mtmp; assert( tmp.checkerboard ==Odd); tmp=src_o-Mtmp; assert( tmp.Checkerboard() ==Odd);
//src_o = tmp; assert(src_o.checkerboard ==Odd); //src_o = tmp; assert(src_o.Checkerboard() ==Odd);
_Matrix.Mooee(tmp,src_o); // Extra factor of "m" in source from dumb choice of matrix norm. _Matrix.Mooee(tmp,src_o); // Extra factor of "m" in source from dumb choice of matrix norm.
////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////
// Call the red-black solver // Call the red-black solver
////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////
std::cout<<GridLogMessage << "SchurRedBlackStaggeredSolver calling the Mpc solver" <<std::endl; std::cout<<GridLogMessage << "SchurRedBlackStaggeredSolver calling the Mpc solver" <<std::endl;
_HermitianRBSolver(_HermOpEO,src_o,sol_o); assert(sol_o.checkerboard==Odd); _HermitianRBSolver(_HermOpEO,src_o,sol_o); assert(sol_o.Checkerboard()==Odd);
std::cout<<GridLogMessage << "SchurRedBlackStaggeredSolver called the Mpc solver" <<std::endl; std::cout<<GridLogMessage << "SchurRedBlackStaggeredSolver called the Mpc solver" <<std::endl;
/////////////////////////////////////////////////// ///////////////////////////////////////////////////
// sol_e = M_ee^-1 * ( src_e - Meo sol_o )... // sol_e = M_ee^-1 * ( src_e - Meo sol_o )...
/////////////////////////////////////////////////// ///////////////////////////////////////////////////
_Matrix.Meooe(sol_o,tmp); assert( tmp.checkerboard ==Even); _Matrix.Meooe(sol_o,tmp); assert( tmp.Checkerboard() ==Even);
src_e = src_e-tmp; assert( src_e.checkerboard ==Even); src_e = src_e-tmp; assert( src_e.Checkerboard() ==Even);
_Matrix.MooeeInv(src_e,sol_e); assert( sol_e.checkerboard ==Even); _Matrix.MooeeInv(src_e,sol_e); assert( sol_e.Checkerboard() ==Even);
std::cout<<GridLogMessage << "SchurRedBlackStaggeredSolver reconstructed other CB" <<std::endl; std::cout<<GridLogMessage << "SchurRedBlackStaggeredSolver reconstructed other CB" <<std::endl;
setCheckerboard(out,sol_e); assert( sol_e.checkerboard ==Even); setCheckerboard(out,sol_e); assert( sol_e.Checkerboard() ==Even);
setCheckerboard(out,sol_o); assert( sol_o.checkerboard ==Odd ); setCheckerboard(out,sol_o); assert( sol_o.Checkerboard() ==Odd );
std::cout<<GridLogMessage << "SchurRedBlackStaggeredSolver inserted solution" <<std::endl; std::cout<<GridLogMessage << "SchurRedBlackStaggeredSolver inserted solution" <<std::endl;
// Verify the unprec residual // Verify the unprec residual
@ -214,28 +214,28 @@ public:
///////////////////////////////////////////////////// /////////////////////////////////////////////////////
// src_o = Mdag * (source_o - Moe MeeInv source_e) // src_o = Mdag * (source_o - Moe MeeInv source_e)
///////////////////////////////////////////////////// /////////////////////////////////////////////////////
_Matrix.MooeeInv(src_e,tmp); assert( tmp.checkerboard ==Even); _Matrix.MooeeInv(src_e,tmp); assert( tmp.Checkerboard() ==Even);
_Matrix.Meooe (tmp,Mtmp); assert( Mtmp.checkerboard ==Odd); _Matrix.Meooe (tmp,Mtmp); assert( Mtmp.Checkerboard() ==Odd);
tmp=src_o-Mtmp; assert( tmp.checkerboard ==Odd); tmp=src_o-Mtmp; assert( tmp.Checkerboard() ==Odd);
// get the right MpcDag // get the right MpcDag
_HermOpEO.MpcDag(tmp,src_o); assert(src_o.checkerboard ==Odd); _HermOpEO.MpcDag(tmp,src_o); assert(src_o.Checkerboard() ==Odd);
////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////
// Call the red-black solver // Call the red-black solver
////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////
std::cout<<GridLogMessage << "SchurRedBlack solver calling the MpcDagMp solver" <<std::endl; std::cout<<GridLogMessage << "SchurRedBlack solver calling the MpcDagMp solver" <<std::endl;
_HermitianRBSolver(_HermOpEO,src_o,sol_o); assert(sol_o.checkerboard==Odd); _HermitianRBSolver(_HermOpEO,src_o,sol_o); assert(sol_o.Checkerboard()==Odd);
/////////////////////////////////////////////////// ///////////////////////////////////////////////////
// sol_e = M_ee^-1 * ( src_e - Meo sol_o )... // sol_e = M_ee^-1 * ( src_e - Meo sol_o )...
/////////////////////////////////////////////////// ///////////////////////////////////////////////////
_Matrix.Meooe(sol_o,tmp); assert( tmp.checkerboard ==Even); _Matrix.Meooe(sol_o,tmp); assert( tmp.Checkerboard() ==Even);
src_e = src_e-tmp; assert( src_e.checkerboard ==Even); src_e = src_e-tmp; assert( src_e.Checkerboard() ==Even);
_Matrix.MooeeInv(src_e,sol_e); assert( sol_e.checkerboard ==Even); _Matrix.MooeeInv(src_e,sol_e); assert( sol_e.Checkerboard() ==Even);
setCheckerboard(out,sol_e); assert( sol_e.checkerboard ==Even); setCheckerboard(out,sol_e); assert( sol_e.Checkerboard() ==Even);
setCheckerboard(out,sol_o); assert( sol_o.checkerboard ==Odd ); setCheckerboard(out,sol_o); assert( sol_o.Checkerboard() ==Odd );
// Verify the unprec residual // Verify the unprec residual
_Matrix.M(out,resid); _Matrix.M(out,resid);
@ -293,30 +293,30 @@ public:
///////////////////////////////////////////////////// /////////////////////////////////////////////////////
// src_o = Mdag * (source_o - Moe MeeInv source_e) // src_o = Mdag * (source_o - Moe MeeInv source_e)
///////////////////////////////////////////////////// /////////////////////////////////////////////////////
_Matrix.MooeeInv(src_e,tmp); assert( tmp.checkerboard ==Even); _Matrix.MooeeInv(src_e,tmp); assert( tmp.Checkerboard() ==Even);
_Matrix.Meooe (tmp,Mtmp); assert( Mtmp.checkerboard ==Odd); _Matrix.Meooe (tmp,Mtmp); assert( Mtmp.Checkerboard() ==Odd);
tmp=src_o-Mtmp; assert( tmp.checkerboard ==Odd); tmp=src_o-Mtmp; assert( tmp.Checkerboard() ==Odd);
// get the right MpcDag // get the right MpcDag
_HermOpEO.MpcDag(tmp,src_o); assert(src_o.checkerboard ==Odd); _HermOpEO.MpcDag(tmp,src_o); assert(src_o.Checkerboard() ==Odd);
////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////
// Call the red-black solver // Call the red-black solver
////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////
std::cout<<GridLogMessage << "SchurRedBlack solver calling the MpcDagMp solver" <<std::endl; std::cout<<GridLogMessage << "SchurRedBlack solver calling the MpcDagMp solver" <<std::endl;
// _HermitianRBSolver(_HermOpEO,src_o,sol_o); assert(sol_o.checkerboard==Odd); // _HermitianRBSolver(_HermOpEO,src_o,sol_o); assert(sol_o.Checkerboard()==Odd);
_HermitianRBSolver(_HermOpEO,src_o,tmp); assert(tmp.checkerboard==Odd); _HermitianRBSolver(_HermOpEO,src_o,tmp); assert(tmp.Checkerboard()==Odd);
_Matrix.MooeeInv(tmp,sol_o); assert( sol_o.checkerboard ==Odd); _Matrix.MooeeInv(tmp,sol_o); assert( sol_o.Checkerboard() ==Odd);
/////////////////////////////////////////////////// ///////////////////////////////////////////////////
// sol_e = M_ee^-1 * ( src_e - Meo sol_o )... // sol_e = M_ee^-1 * ( src_e - Meo sol_o )...
/////////////////////////////////////////////////// ///////////////////////////////////////////////////
_Matrix.Meooe(sol_o,tmp); assert( tmp.checkerboard ==Even); _Matrix.Meooe(sol_o,tmp); assert( tmp.Checkerboard() ==Even);
src_e = src_e-tmp; assert( src_e.checkerboard ==Even); src_e = src_e-tmp; assert( src_e.Checkerboard() ==Even);
_Matrix.MooeeInv(src_e,sol_e); assert( sol_e.checkerboard ==Even); _Matrix.MooeeInv(src_e,sol_e); assert( sol_e.Checkerboard() ==Even);
setCheckerboard(out,sol_e); assert( sol_e.checkerboard ==Even); setCheckerboard(out,sol_e); assert( sol_e.Checkerboard() ==Even);
setCheckerboard(out,sol_o); assert( sol_o.checkerboard ==Odd ); setCheckerboard(out,sol_o); assert( sol_o.Checkerboard() ==Odd );
// Verify the unprec residual // Verify the unprec residual
_Matrix.M(out,resid); _Matrix.M(out,resid);
@ -372,31 +372,31 @@ public:
///////////////////////////////////////////////////// /////////////////////////////////////////////////////
// src_o = Mdag * (source_o - Moe MeeInv source_e) // src_o = Mdag * (source_o - Moe MeeInv source_e)
///////////////////////////////////////////////////// /////////////////////////////////////////////////////
_Matrix.MooeeInv(src_e,tmp); assert( tmp.checkerboard ==Even); _Matrix.MooeeInv(src_e,tmp); assert( tmp.Checkerboard() ==Even);
_Matrix.Meooe (tmp,Mtmp); assert( Mtmp.checkerboard ==Odd); _Matrix.Meooe (tmp,Mtmp); assert( Mtmp.Checkerboard() ==Odd);
tmp=src_o-Mtmp; assert( tmp.checkerboard ==Odd); tmp=src_o-Mtmp; assert( tmp.Checkerboard() ==Odd);
// get the right MpcDag // get the right MpcDag
_HermOpEO.MpcDag(tmp,src_o); assert(src_o.checkerboard ==Odd); _HermOpEO.MpcDag(tmp,src_o); assert(src_o.Checkerboard() ==Odd);
////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////
// Call the red-black solver // Call the red-black solver
////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////
std::cout<<GridLogMessage << "SchurRedBlack solver calling the MpcDagMp solver" <<std::endl; std::cout<<GridLogMessage << "SchurRedBlack solver calling the MpcDagMp solver" <<std::endl;
// _HermitianRBSolver(_HermOpEO,src_o,sol_o); assert(sol_o.checkerboard==Odd); // _HermitianRBSolver(_HermOpEO,src_o,sol_o); assert(sol_o.Checkerboard()==Odd);
// _HermitianRBSolver(_HermOpEO,src_o,tmp); assert(tmp.checkerboard==Odd); // _HermitianRBSolver(_HermOpEO,src_o,tmp); assert(tmp.Checkerboard()==Odd);
_HermitianRBSolver(src_o,tmp); assert(tmp.checkerboard==Odd); _HermitianRBSolver(src_o,tmp); assert(tmp.Checkerboard()==Odd);
_Matrix.MooeeInv(tmp,sol_o); assert( sol_o.checkerboard ==Odd); _Matrix.MooeeInv(tmp,sol_o); assert( sol_o.Checkerboard() ==Odd);
/////////////////////////////////////////////////// ///////////////////////////////////////////////////
// sol_e = M_ee^-1 * ( src_e - Meo sol_o )... // sol_e = M_ee^-1 * ( src_e - Meo sol_o )...
/////////////////////////////////////////////////// ///////////////////////////////////////////////////
_Matrix.Meooe(sol_o,tmp); assert( tmp.checkerboard ==Even); _Matrix.Meooe(sol_o,tmp); assert( tmp.Checkerboard() ==Even);
src_e = src_e-tmp; assert( src_e.checkerboard ==Even); src_e = src_e-tmp; assert( src_e.Checkerboard() ==Even);
_Matrix.MooeeInv(src_e,sol_e); assert( sol_e.checkerboard ==Even); _Matrix.MooeeInv(src_e,sol_e); assert( sol_e.Checkerboard() ==Even);
setCheckerboard(out,sol_e); assert( sol_e.checkerboard ==Even); setCheckerboard(out,sol_e); assert( sol_e.Checkerboard() ==Even);
setCheckerboard(out,sol_o); assert( sol_o.checkerboard ==Odd ); setCheckerboard(out,sol_o); assert( sol_o.Checkerboard() ==Odd );
// Verify the unprec residual // Verify the unprec residual
_Matrix.M(out,resid); _Matrix.M(out,resid);