mirror of
https://github.com/paboyle/Grid.git
synced 2025-04-10 14:10:46 +01:00
Save two more files
This commit is contained in:
Daniel Richtmann
Daniel Richtmann
parent
789e892865
commit
1ab8d5cc13
197
lib/algorithms/iterative/MinimalResidual.h
Normal file
197
lib/algorithms/iterative/MinimalResidual.h
Normal file
@ -0,0 +1,197 @@
|
||||
/*************************************************************************************
|
||||
|
||||
Grid physics library, www.github.com/paboyle/Grid
|
||||
|
||||
Source file: ./lib/algorithms/iterative/MinimalResidual.h
|
||||
|
||||
Copyright (C) 2015
|
||||
|
||||
Author: Azusa Yamaguchi <ayamaguc@staffmail.ed.ac.uk>
|
||||
Author: Peter Boyle <paboyle@ph.ed.ac.uk>
|
||||
Author: paboyle <paboyle@ph.ed.ac.uk>
|
||||
|
||||
This program is free software; you can redistribute it and/or modify
|
||||
it under the terms of the GNU General Public License as published by
|
||||
the Free Software Foundation; either version 2 of the License, or
|
||||
(at your option) any later version.
|
||||
|
||||
This program is distributed in the hope that it will be useful,
|
||||
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||
GNU General Public License for more details.
|
||||
|
||||
You should have received a copy of the GNU General Public License along
|
||||
with this program; if not, write to the Free Software Foundation, Inc.,
|
||||
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
|
||||
|
||||
See the full license in the file "LICENSE" in the top level distribution
|
||||
directory
|
||||
*************************************************************************************/
|
||||
/* END LEGAL */
|
||||
#ifndef GRID_MINIMAL_RESIDUAL_H
|
||||
#define GRID_MINIMAL_RESIDUAL_H
|
||||
|
||||
namespace Grid {
|
||||
|
||||
/////////////////////////////////////////////////////////////
|
||||
// Base classes for iterative processes based on operators
|
||||
// single input vec, single output vec.
|
||||
/////////////////////////////////////////////////////////////
|
||||
|
||||
template <class Field>
|
||||
class MinimalResidual : public OperatorFunction<Field> {
|
||||
public:
|
||||
bool ErrorOnNoConverge; // throw an assert when the MR fails to converge.
|
||||
// Defaults true.
|
||||
RealD Tolerance;
|
||||
Integer MaxIterations;
|
||||
Integer IterationsToComplete; //Number of iterations the MR took to finish. Filled in upon completion
|
||||
|
||||
MinimalResidual(RealD tol, Integer maxit, bool err_on_no_conv = true)
|
||||
: Tolerance(tol),
|
||||
MaxIterations(maxit),
|
||||
ErrorOnNoConverge(err_on_no_conv){};
|
||||
|
||||
void operator()(LinearOperatorBase<Field> &Linop, const Field &src,
|
||||
Field &psi) {
|
||||
psi.checkerboard = src.checkerboard; // Check
|
||||
conformable(psi, src);
|
||||
|
||||
Field p {src};
|
||||
Field matrixTimesPsi {src};
|
||||
Field r {src};
|
||||
|
||||
RealD alpha {};
|
||||
|
||||
// Initial residual computation & set up
|
||||
RealD guess = norm2(psi);
|
||||
assert(std::isnan(guess) == 0);
|
||||
|
||||
Linop.HermOp(psi, matrixTimesPsi);
|
||||
|
||||
r = src - matrixTimesPsi;
|
||||
|
||||
Linop.HermOp(r, p);
|
||||
|
||||
alpha = innerProduct(p,r) / innerProduct(p,p);
|
||||
psi = psi + alpha * r;
|
||||
r = r - alpha * p;
|
||||
|
||||
Linop.HermOp(r, p);
|
||||
|
||||
|
||||
////////////////////////////////////////////////////////////////////////////////
|
||||
////////////////////////////////////////////////////////////////////////////////
|
||||
|
||||
// RealD cp, c, a, d, b, ssq, qq, b_pred;
|
||||
|
||||
Field p(src);
|
||||
Field matrixTimesPsi(src);
|
||||
// Field r(src);
|
||||
|
||||
// Initial residual computation & set up
|
||||
RealD guess = norm2(psi);
|
||||
assert(std::isnan(guess) == 0);
|
||||
|
||||
|
||||
Linop.HermOpAndNorm(psi, matrixTimesPsi, d, b);
|
||||
|
||||
|
||||
r = src - matrixTimesPsi;
|
||||
p = matrixTimesPsi;
|
||||
|
||||
a = norm2(p);
|
||||
cp = a;
|
||||
ssq = norm2(src);
|
||||
|
||||
std::cout << GridLogIterative << std::setprecision(4)
|
||||
<< "MinimalResidual: guess " << guess << std::endl;
|
||||
std::cout << GridLogIterative << std::setprecision(4)
|
||||
<< "MinimalResidual: src " << ssq << std::endl;
|
||||
std::cout << GridLogIterative << std::setprecision(4)
|
||||
<< "MinimalResidual: mp " << d << std::endl;
|
||||
std::cout << GridLogIterative << std::setprecision(4)
|
||||
<< "MinimalResidual: matrixTimesPsi " << b << std::endl;
|
||||
std::cout << GridLogIterative << std::setprecision(4)
|
||||
<< "MinimalResidual: cp,r " << cp << std::endl;
|
||||
std::cout << GridLogIterative << std::setprecision(4)
|
||||
<< "MinimalResidual: p " << a << std::endl;
|
||||
|
||||
RealD rsq = Tolerance * Tolerance * ssq;
|
||||
|
||||
// Check if guess is really REALLY good :)
|
||||
if (cp <= rsq) {
|
||||
return;
|
||||
}
|
||||
|
||||
std::cout << GridLogIterative << std::setprecision(4)
|
||||
<< "MinimalResidual: k=0 residual " << cp << " target " << rsq
|
||||
<< std::endl;
|
||||
|
||||
GridStopWatch LinalgTimer;
|
||||
GridStopWatch MatrixTimer;
|
||||
GridStopWatch SolverTimer;
|
||||
|
||||
SolverTimer.Start();
|
||||
int k;
|
||||
for (k = 1; k <= MaxIterations; k++) {
|
||||
c = cp;
|
||||
|
||||
MatrixTimer.Start();
|
||||
Linop.HermOpAndNorm(p, matrixTimesPsi, d, qq);
|
||||
MatrixTimer.Stop();
|
||||
|
||||
LinalgTimer.Start();
|
||||
// RealD qqck = norm2(matrixTimesPsi);
|
||||
// ComplexD dck = innerProduct(p,matrixTimesPsi);
|
||||
|
||||
a = c / d;
|
||||
b_pred = a * (a * qq - d) / c;
|
||||
|
||||
cp = axpy_norm(r, -a, matrixTimesPsi, r);
|
||||
b = cp / c;
|
||||
|
||||
// Fuse these loops ; should be really easy
|
||||
psi = a * p + psi;
|
||||
p = p * b + r;
|
||||
|
||||
LinalgTimer.Stop();
|
||||
std::cout << GridLogIterative << "MinimalResidual: Iteration " << k
|
||||
<< " residual " << cp << " target " << rsq << std::endl;
|
||||
|
||||
// Stopping condition
|
||||
if (cp <= rsq) {
|
||||
SolverTimer.Stop();
|
||||
Linop.HermOpAndNorm(psi, matrixTimesPsi, d, qq);
|
||||
p = matrixTimesPsi - src;
|
||||
|
||||
RealD matrixTimesPsiNorm = sqrt(norm2(matrixTimesPsi));
|
||||
RealD psinorm = sqrt(norm2(psi));
|
||||
RealD srcnorm = sqrt(norm2(src));
|
||||
RealD resnorm = sqrt(norm2(p));
|
||||
RealD true_residual = resnorm / srcnorm;
|
||||
|
||||
std::cout << GridLogMessage
|
||||
<< "MinimalResidual: Converged on iteration " << k << std::endl;
|
||||
std::cout << GridLogMessage << "Computed residual " << sqrt(cp / ssq)
|
||||
<< " true residual " << true_residual << " target "
|
||||
<< Tolerance << std::endl;
|
||||
std::cout << GridLogMessage << "Time elapsed: Iterations "
|
||||
<< SolverTimer.Elapsed() << " Matrix "
|
||||
<< MatrixTimer.Elapsed() << " Linalg "
|
||||
<< LinalgTimer.Elapsed();
|
||||
std::cout << std::endl;
|
||||
|
||||
if (ErrorOnNoConverge) assert(true_residual / Tolerance < 10000.0);
|
||||
IterationsToComplete = k;
|
||||
return;
|
||||
}
|
||||
}
|
||||
std::cout << GridLogMessage << "MinimalResidual did NOT converge"
|
||||
<< std::endl;
|
||||
if (ErrorOnNoConverge) assert(0);
|
||||
IterationsToComplete = k;
|
||||
}
|
||||
};
|
||||
}
|
||||
#endif
|
||||
806
tests/solver/Test_wilson_ddalphaamg.cc
Normal file
806
tests/solver/Test_wilson_ddalphaamg.cc
Normal file
@ -0,0 +1,806 @@
|
||||
/*************************************************************************************
|
||||
|
||||
Grid physics library, www.github.com/paboyle/Grid
|
||||
|
||||
Source file: ./tests/Test_dwf_hdcr.cc
|
||||
|
||||
Copyright (C) 2015
|
||||
|
||||
Author: Daniel Richtmann <daniel.richtmann@ur.de>
|
||||
|
||||
This program is free software; you can redistribute it and/or modify
|
||||
it under the terms of the GNU General Public License as published by
|
||||
the Free Software Foundation; either version 2 of the License, or
|
||||
(at your option) any later version.
|
||||
|
||||
This program is distributed in the hope that it will be useful,
|
||||
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||
GNU General Public License for more details.
|
||||
|
||||
You should have received a copy of the GNU General Public License along
|
||||
with this program; if not, write to the Free Software Foundation, Inc.,
|
||||
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
|
||||
|
||||
See the full license in the file "LICENSE" in the top level distribution directory
|
||||
*************************************************************************************/
|
||||
/* END LEGAL */
|
||||
#include <Grid/Grid.h>
|
||||
#include <Grid/algorithms/iterative/PrecGeneralisedConjugateResidual.h>
|
||||
//#include <algorithms/iterative/PrecConjugateResidual.h>
|
||||
|
||||
using namespace std;
|
||||
using namespace Grid;
|
||||
using namespace Grid::QCD;
|
||||
|
||||
class myclass: Serializable {
|
||||
public:
|
||||
|
||||
GRID_SERIALIZABLE_CLASS_MEMBERS(myclass,
|
||||
int, domaindecompose,
|
||||
int, domainsize,
|
||||
int, order,
|
||||
int, Ls,
|
||||
double, mq,
|
||||
double, lo,
|
||||
double, hi,
|
||||
int, steps);
|
||||
|
||||
myclass(){};
|
||||
|
||||
};
|
||||
myclass params;
|
||||
|
||||
RealD InverseApproximation(RealD x){
|
||||
return 1.0/x;
|
||||
}
|
||||
|
||||
template<class Fobj,class CComplex,int nbasis, class Matrix>
|
||||
class MultiGridPreconditioner : public LinearFunction< Lattice<Fobj> > {
|
||||
public:
|
||||
|
||||
typedef Aggregation<Fobj,CComplex,nbasis> Aggregates;
|
||||
typedef CoarsenedMatrix<Fobj,CComplex,nbasis> CoarseOperator;
|
||||
|
||||
typedef typename Aggregation<Fobj,CComplex,nbasis>::siteVector siteVector;
|
||||
typedef typename Aggregation<Fobj,CComplex,nbasis>::CoarseScalar CoarseScalar;
|
||||
typedef typename Aggregation<Fobj,CComplex,nbasis>::CoarseVector CoarseVector;
|
||||
typedef typename Aggregation<Fobj,CComplex,nbasis>::CoarseMatrix CoarseMatrix;
|
||||
typedef typename Aggregation<Fobj,CComplex,nbasis>::FineField FineField;
|
||||
typedef LinearOperatorBase<FineField> FineOperator;
|
||||
|
||||
Aggregates & _Aggregates;
|
||||
CoarseOperator & _CoarseOperator;
|
||||
Matrix & _FineMatrix;
|
||||
FineOperator & _FineOperator;
|
||||
Matrix & _SmootherMatrix;
|
||||
FineOperator & _SmootherOperator;
|
||||
|
||||
// Constructor
|
||||
MultiGridPreconditioner(Aggregates &Agg, CoarseOperator &Coarse,
|
||||
FineOperator &Fine,Matrix &FineMatrix,
|
||||
FineOperator &Smooth,Matrix &SmootherMatrix)
|
||||
: _Aggregates(Agg),
|
||||
_CoarseOperator(Coarse),
|
||||
_FineOperator(Fine),
|
||||
_FineMatrix(FineMatrix),
|
||||
_SmootherOperator(Smooth),
|
||||
_SmootherMatrix(SmootherMatrix)
|
||||
{
|
||||
}
|
||||
|
||||
void PowerMethod(const FineField &in) {
|
||||
|
||||
FineField p1(in._grid);
|
||||
FineField p2(in._grid);
|
||||
|
||||
MdagMLinearOperator<Matrix,FineField> fMdagMOp(_FineMatrix);
|
||||
|
||||
p1=in;
|
||||
RealD absp2;
|
||||
for(int i=0;i<20;i++){
|
||||
RealD absp1=std::sqrt(norm2(p1));
|
||||
fMdagMOp.HermOp(p1,p2);// this is the G5 herm bit
|
||||
// _FineOperator.Op(p1,p2);// this is the G5 herm bit
|
||||
RealD absp2=std::sqrt(norm2(p2));
|
||||
if(i%10==9)
|
||||
std::cout<<GridLogMessage << "Power method on mdagm "<<i<<" " << absp2/absp1<<std::endl;
|
||||
p1=p2*(1.0/std::sqrt(absp2));
|
||||
}
|
||||
}
|
||||
|
||||
void operator()(const FineField &in, FineField & out) {
|
||||
if ( params.domaindecompose ) {
|
||||
operatorSAP(in,out);
|
||||
} else {
|
||||
operatorCheby(in,out);
|
||||
}
|
||||
}
|
||||
|
||||
////////////////////////////////////////////////////////////////////////
|
||||
// ADEF2: [PTM+Q] in = [1 - Q A] M in + Q in = Min + Q [ in -A Min]
|
||||
// ADEF1: [MP+Q ] in =M [1 - A Q] in + Q in
|
||||
////////////////////////////////////////////////////////////////////////
|
||||
#if 1
|
||||
void operatorADEF2(const FineField &in, FineField & out) {
|
||||
|
||||
CoarseVector Csrc(_CoarseOperator.Grid());
|
||||
CoarseVector Ctmp(_CoarseOperator.Grid());
|
||||
CoarseVector Csol(_CoarseOperator.Grid());
|
||||
|
||||
ConjugateGradient<CoarseVector> CG(1.0e-10,100000);
|
||||
ConjugateGradient<FineField> fCG(3.0e-2,1000);
|
||||
|
||||
HermitianLinearOperator<CoarseOperator,CoarseVector> HermOp(_CoarseOperator);
|
||||
MdagMLinearOperator<CoarseOperator,CoarseVector> MdagMOp(_CoarseOperator);
|
||||
MdagMLinearOperator<Matrix,FineField> fMdagMOp(_FineMatrix);
|
||||
|
||||
FineField tmp(in._grid);
|
||||
FineField res(in._grid);
|
||||
FineField Min(in._grid);
|
||||
|
||||
// Monitor completeness of low mode space
|
||||
_Aggregates.ProjectToSubspace (Csrc,in);
|
||||
_Aggregates.PromoteFromSubspace(Csrc,out);
|
||||
std::cout<<GridLogMessage<<"Coarse Grid Preconditioner\nCompleteness in: "<<std::sqrt(norm2(out)/norm2(in))<<std::endl;
|
||||
|
||||
// [PTM+Q] in = [1 - Q A] M in + Q in = Min + Q [ in -A Min]
|
||||
_FineOperator.Op(in,tmp);// this is the G5 herm bit
|
||||
fCG(fMdagMOp,tmp,Min); // solves MdagM = g5 M g5M
|
||||
|
||||
// Monitor completeness of low mode space
|
||||
_Aggregates.ProjectToSubspace (Csrc,Min);
|
||||
_Aggregates.PromoteFromSubspace(Csrc,out);
|
||||
std::cout<<GridLogMessage<<"Completeness Min: "<<std::sqrt(norm2(out)/norm2(Min))<<std::endl;
|
||||
|
||||
_FineOperator.Op(Min,tmp);
|
||||
tmp = in - tmp; // in - A Min
|
||||
|
||||
Csol=zero;
|
||||
_Aggregates.ProjectToSubspace (Csrc,tmp);
|
||||
HermOp.AdjOp(Csrc,Ctmp);// Normal equations
|
||||
CG(MdagMOp,Ctmp,Csol);
|
||||
|
||||
HermOp.Op(Csol,Ctmp);
|
||||
Ctmp=Ctmp-Csrc;
|
||||
std::cout<<GridLogMessage<<"coarse space true residual "<<std::sqrt(norm2(Ctmp)/norm2(Csrc))<<std::endl;
|
||||
_Aggregates.PromoteFromSubspace(Csol,out);
|
||||
|
||||
_FineOperator.Op(out,res);
|
||||
res=res-tmp;
|
||||
std::cout<<GridLogMessage<<"promoted sol residual "<<std::sqrt(norm2(res)/norm2(tmp))<<std::endl;
|
||||
_Aggregates.ProjectToSubspace (Csrc,res);
|
||||
std::cout<<GridLogMessage<<"coarse space proj of residual "<<norm2(Csrc)<<std::endl;
|
||||
|
||||
|
||||
out = out+Min; // additive coarse space correction
|
||||
// out = Min; // no additive coarse space correction
|
||||
|
||||
_FineOperator.Op(out,tmp);
|
||||
tmp=tmp-in; // tmp is new residual
|
||||
|
||||
std::cout<<GridLogMessage<< " Preconditioner in " << norm2(in)<<std::endl;
|
||||
std::cout<<GridLogMessage<< " Preconditioner out " << norm2(out)<<std::endl;
|
||||
std::cout<<GridLogMessage<<"preconditioner thinks residual is "<<std::sqrt(norm2(tmp)/norm2(in))<<std::endl;
|
||||
|
||||
}
|
||||
#endif
|
||||
// ADEF1: [MP+Q ] in =M [1 - A Q] in + Q in
|
||||
#if 1
|
||||
void operatorADEF1(const FineField &in, FineField & out) {
|
||||
|
||||
CoarseVector Csrc(_CoarseOperator.Grid());
|
||||
CoarseVector Ctmp(_CoarseOperator.Grid());
|
||||
CoarseVector Csol(_CoarseOperator.Grid()); Csol=zero;
|
||||
|
||||
ConjugateGradient<CoarseVector> CG(1.0e-10,100000);
|
||||
ConjugateGradient<FineField> fCG(3.0e-2,1000);
|
||||
|
||||
HermitianLinearOperator<CoarseOperator,CoarseVector> HermOp(_CoarseOperator);
|
||||
MdagMLinearOperator<CoarseOperator,CoarseVector> MdagMOp(_CoarseOperator);
|
||||
ShiftedMdagMLinearOperator<Matrix,FineField> fMdagMOp(_FineMatrix,0.1);
|
||||
|
||||
FineField tmp(in._grid);
|
||||
FineField res(in._grid);
|
||||
FineField Qin(in._grid);
|
||||
|
||||
// Monitor completeness of low mode space
|
||||
// _Aggregates.ProjectToSubspace (Csrc,in);
|
||||
// _Aggregates.PromoteFromSubspace(Csrc,out);
|
||||
// std::cout<<GridLogMessage<<"Coarse Grid Preconditioner\nCompleteness in: "<<std::sqrt(norm2(out)/norm2(in))<<std::endl;
|
||||
|
||||
_Aggregates.ProjectToSubspace (Csrc,in);
|
||||
HermOp.AdjOp(Csrc,Ctmp);// Normal equations
|
||||
CG(MdagMOp,Ctmp,Csol);
|
||||
_Aggregates.PromoteFromSubspace(Csol,Qin);
|
||||
|
||||
// Qin=0;
|
||||
_FineOperator.Op(Qin,tmp);// A Q in
|
||||
tmp = in - tmp; // in - A Q in
|
||||
|
||||
_FineOperator.Op(tmp,res);// this is the G5 herm bit
|
||||
fCG(fMdagMOp,res,out); // solves MdagM = g5 M g5M
|
||||
|
||||
out = out + Qin;
|
||||
|
||||
_FineOperator.Op(out,tmp);
|
||||
tmp=tmp-in; // tmp is new residual
|
||||
|
||||
std::cout<<GridLogMessage<<"preconditioner thinks residual is "<<std::sqrt(norm2(tmp)/norm2(in))<<std::endl;
|
||||
|
||||
}
|
||||
#endif
|
||||
|
||||
void SAP (const FineField & src,FineField & psi){
|
||||
|
||||
Lattice<iScalar<vInteger> > coor(src._grid);
|
||||
Lattice<iScalar<vInteger> > subset(src._grid);
|
||||
|
||||
FineField r(src._grid);
|
||||
FineField zz(src._grid); zz=zero;
|
||||
FineField vec1(src._grid);
|
||||
FineField vec2(src._grid);
|
||||
|
||||
const Integer block=params.domainsize;
|
||||
|
||||
subset=zero;
|
||||
for(int mu=0;mu<Nd;mu++){
|
||||
LatticeCoordinate(coor,mu+1);
|
||||
coor = div(coor,block);
|
||||
subset = subset+coor;
|
||||
}
|
||||
subset = mod(subset,(Integer)2);
|
||||
|
||||
ShiftedMdagMLinearOperator<Matrix,FineField> fMdagMOp(_SmootherMatrix,0.0);
|
||||
Chebyshev<FineField> Cheby (params.lo,params.hi,params.order,InverseApproximation);
|
||||
|
||||
RealD resid;
|
||||
for(int i=0;i<params.steps;i++){
|
||||
|
||||
// Even domain residual
|
||||
_FineOperator.Op(psi,vec1);// this is the G5 herm bit
|
||||
r= src - vec1 ;
|
||||
resid = norm2(r) /norm2(src);
|
||||
std::cout << "SAP "<<i<<" resid "<<resid<<std::endl;
|
||||
|
||||
// Even domain solve
|
||||
r= where(subset==(Integer)0,r,zz);
|
||||
_SmootherOperator.AdjOp(r,vec1);
|
||||
Cheby(fMdagMOp,vec1,vec2); // solves MdagM = g5 M g5M
|
||||
psi = psi + vec2;
|
||||
|
||||
// Odd domain residual
|
||||
_FineOperator.Op(psi,vec1);// this is the G5 herm bit
|
||||
r= src - vec1 ;
|
||||
r= where(subset==(Integer)1,r,zz);
|
||||
|
||||
resid = norm2(r) /norm2(src);
|
||||
std::cout << "SAP "<<i<<" resid "<<resid<<std::endl;
|
||||
|
||||
// Odd domain solve
|
||||
_SmootherOperator.AdjOp(r,vec1);
|
||||
Cheby(fMdagMOp,vec1,vec2); // solves MdagM = g5 M g5M
|
||||
psi = psi + vec2;
|
||||
|
||||
_FineOperator.Op(psi,vec1);// this is the G5 herm bit
|
||||
r= src - vec1 ;
|
||||
resid = norm2(r) /norm2(src);
|
||||
std::cout << "SAP "<<i<<" resid "<<resid<<std::endl;
|
||||
|
||||
}
|
||||
|
||||
};
|
||||
|
||||
void SmootherTest (const FineField & in){
|
||||
|
||||
FineField vec1(in._grid);
|
||||
FineField vec2(in._grid);
|
||||
RealD lo[3] = { 0.5, 1.0, 2.0};
|
||||
|
||||
// MdagMLinearOperator<Matrix,FineField> fMdagMOp(_FineMatrix);
|
||||
ShiftedMdagMLinearOperator<Matrix,FineField> fMdagMOp(_SmootherMatrix,0.0);
|
||||
|
||||
RealD Ni,r;
|
||||
|
||||
Ni = norm2(in);
|
||||
|
||||
for(int ilo=0;ilo<3;ilo++){
|
||||
for(int ord=5;ord<50;ord*=2){
|
||||
|
||||
_SmootherOperator.AdjOp(in,vec1);
|
||||
|
||||
Chebyshev<FineField> Cheby (lo[ilo],70.0,ord,InverseApproximation);
|
||||
Cheby(fMdagMOp,vec1,vec2); // solves MdagM = g5 M g5M
|
||||
|
||||
_FineOperator.Op(vec2,vec1);// this is the G5 herm bit
|
||||
vec1 = in - vec1; // tmp = in - A Min
|
||||
r=norm2(vec1);
|
||||
std::cout<<GridLogMessage << "Smoother resid "<<std::sqrt(r/Ni)<<std::endl;
|
||||
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
void operatorCheby(const FineField &in, FineField & out) {
|
||||
|
||||
CoarseVector Csrc(_CoarseOperator.Grid());
|
||||
CoarseVector Ctmp(_CoarseOperator.Grid());
|
||||
CoarseVector Csol(_CoarseOperator.Grid()); Csol=zero;
|
||||
|
||||
ConjugateGradient<CoarseVector> CG(3.0e-3,100000);
|
||||
// ConjugateGradient<FineField> fCG(3.0e-2,1000);
|
||||
|
||||
HermitianLinearOperator<CoarseOperator,CoarseVector> HermOp(_CoarseOperator);
|
||||
MdagMLinearOperator<CoarseOperator,CoarseVector> MdagMOp(_CoarseOperator);
|
||||
// MdagMLinearOperator<Matrix,FineField> fMdagMOp(_FineMatrix);
|
||||
ShiftedMdagMLinearOperator<Matrix,FineField> fMdagMOp(_SmootherMatrix,0.0);
|
||||
|
||||
FineField vec1(in._grid);
|
||||
FineField vec2(in._grid);
|
||||
|
||||
// Chebyshev<FineField> Cheby (0.5,70.0,30,InverseApproximation);
|
||||
// Chebyshev<FineField> ChebyAccu(0.5,70.0,30,InverseApproximation);
|
||||
Chebyshev<FineField> Cheby (params.lo,params.hi,params.order,InverseApproximation);
|
||||
Chebyshev<FineField> ChebyAccu(params.lo,params.hi,params.order,InverseApproximation);
|
||||
// Cheby.JacksonSmooth();
|
||||
// ChebyAccu.JacksonSmooth();
|
||||
|
||||
// _Aggregates.ProjectToSubspace (Csrc,in);
|
||||
// _Aggregates.PromoteFromSubspace(Csrc,out);
|
||||
// std::cout<<GridLogMessage<<"Completeness: "<<std::sqrt(norm2(out)/norm2(in))<<std::endl;
|
||||
|
||||
// ofstream fout("smoother");
|
||||
// Cheby.csv(fout);
|
||||
|
||||
// V11 multigrid.
|
||||
// Use a fixed chebyshev and hope hermiticity helps.
|
||||
|
||||
// To make a working smoother for indefinite operator
|
||||
// must multiply by "Mdag" (ouch loses all low mode content)
|
||||
// and apply to poly approx of (mdagm)^-1.
|
||||
// so that we end up with an odd polynomial.
|
||||
|
||||
RealD Ni = norm2(in);
|
||||
|
||||
_SmootherOperator.AdjOp(in,vec1);// this is the G5 herm bit
|
||||
ChebyAccu(fMdagMOp,vec1,out); // solves MdagM = g5 M g5M
|
||||
|
||||
std::cout<<GridLogMessage << "Smoother norm "<<norm2(out)<<std::endl;
|
||||
|
||||
// Update with residual for out
|
||||
_FineOperator.Op(out,vec1);// this is the G5 herm bit
|
||||
vec1 = in - vec1; // tmp = in - A Min
|
||||
|
||||
RealD r = norm2(vec1);
|
||||
|
||||
std::cout<<GridLogMessage << "Smoother resid "<<std::sqrt(r/Ni)<< " " << r << " " << Ni <<std::endl;
|
||||
|
||||
_Aggregates.ProjectToSubspace (Csrc,vec1);
|
||||
HermOp.AdjOp(Csrc,Ctmp);// Normal equations
|
||||
CG(MdagMOp,Ctmp,Csol);
|
||||
_Aggregates.PromoteFromSubspace(Csol,vec1); // Ass^{-1} [in - A Min]_s
|
||||
// Q = Q[in - A Min]
|
||||
out = out+vec1;
|
||||
|
||||
// Three preconditioner smoothing -- hermitian if C3 = C1
|
||||
// Recompute error
|
||||
_FineOperator.Op(out,vec1);// this is the G5 herm bit
|
||||
vec1 = in - vec1; // tmp = in - A Min
|
||||
r=norm2(vec1);
|
||||
|
||||
std::cout<<GridLogMessage << "Coarse resid "<<std::sqrt(r/Ni)<<std::endl;
|
||||
|
||||
// Reapply smoother
|
||||
_SmootherOperator.Op(vec1,vec2); // this is the G5 herm bit
|
||||
ChebyAccu(fMdagMOp,vec2,vec1); // solves MdagM = g5 M g5M
|
||||
|
||||
out =out+vec1;
|
||||
vec1 = in - vec1; // tmp = in - A Min
|
||||
r=norm2(vec1);
|
||||
std::cout<<GridLogMessage << "Smoother resid "<<std::sqrt(r/Ni)<<std::endl;
|
||||
|
||||
}
|
||||
|
||||
void operatorSAP(const FineField &in, FineField & out) {
|
||||
|
||||
CoarseVector Csrc(_CoarseOperator.Grid());
|
||||
CoarseVector Ctmp(_CoarseOperator.Grid());
|
||||
CoarseVector Csol(_CoarseOperator.Grid()); Csol=zero;
|
||||
|
||||
ConjugateGradient<CoarseVector> CG(1.0e-3,100000);
|
||||
|
||||
HermitianLinearOperator<CoarseOperator,CoarseVector> HermOp(_CoarseOperator);
|
||||
MdagMLinearOperator<CoarseOperator,CoarseVector> MdagMOp(_CoarseOperator);
|
||||
|
||||
FineField vec1(in._grid);
|
||||
FineField vec2(in._grid);
|
||||
|
||||
_Aggregates.ProjectToSubspace (Csrc,in);
|
||||
_Aggregates.PromoteFromSubspace(Csrc,out);
|
||||
std::cout<<GridLogMessage<<"Completeness: "<<std::sqrt(norm2(out)/norm2(in))<<std::endl;
|
||||
|
||||
|
||||
// To make a working smoother for indefinite operator
|
||||
// must multiply by "Mdag" (ouch loses all low mode content)
|
||||
// and apply to poly approx of (mdagm)^-1.
|
||||
// so that we end up with an odd polynomial.
|
||||
SAP(in,out);
|
||||
|
||||
// Update with residual for out
|
||||
_FineOperator.Op(out,vec1);// this is the G5 herm bit
|
||||
vec1 = in - vec1; // tmp = in - A Min
|
||||
|
||||
RealD r = norm2(vec1);
|
||||
RealD Ni = norm2(in);
|
||||
std::cout<<GridLogMessage << "SAP resid "<<std::sqrt(r/Ni)<< " " << r << " " << Ni <<std::endl;
|
||||
|
||||
_Aggregates.ProjectToSubspace (Csrc,vec1);
|
||||
HermOp.AdjOp(Csrc,Ctmp);// Normal equations
|
||||
CG(MdagMOp,Ctmp,Csol);
|
||||
_Aggregates.PromoteFromSubspace(Csol,vec1); // Ass^{-1} [in - A Min]_s
|
||||
// Q = Q[in - A Min]
|
||||
out = out+vec1;
|
||||
|
||||
// Three preconditioner smoothing -- hermitian if C3 = C1
|
||||
// Recompute error
|
||||
_FineOperator.Op(out,vec1);// this is the G5 herm bit
|
||||
vec1 = in - vec1; // tmp = in - A Min
|
||||
r=norm2(vec1);
|
||||
|
||||
std::cout<<GridLogMessage << "Coarse resid "<<std::sqrt(r/Ni)<<std::endl;
|
||||
|
||||
// Reapply smoother
|
||||
SAP(vec1,vec2);
|
||||
out =out+vec2;
|
||||
|
||||
|
||||
// Update with residual for out
|
||||
_FineOperator.Op(out,vec1);// this is the G5 herm bit
|
||||
vec1 = in - vec1; // tmp = in - A Min
|
||||
|
||||
r = norm2(vec1);
|
||||
Ni = norm2(in);
|
||||
std::cout<<GridLogMessage << "SAP resid(post) "<<std::sqrt(r/Ni)<< " " << r << " " << Ni <<std::endl;
|
||||
|
||||
}
|
||||
|
||||
};
|
||||
|
||||
#if 0
|
||||
int main (int argc, char ** argv)
|
||||
{
|
||||
Grid_init(&argc,&argv);
|
||||
|
||||
const int Ls=params.Ls;
|
||||
|
||||
GridCartesian * UGrid = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi());
|
||||
GridRedBlackCartesian * UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
|
||||
|
||||
///////////////////////////////////////////////////
|
||||
// Construct a coarsened grid; utility for this?
|
||||
///////////////////////////////////////////////////
|
||||
std::vector<int> block ({2,2,2,2});
|
||||
const int nbasis= 32;
|
||||
|
||||
std::vector<int> clatt = GridDefaultLatt();
|
||||
for(int d=0;d<clatt.size();d++){
|
||||
clatt[d] = clatt[d]/block[d];
|
||||
}
|
||||
GridCartesian *Coarse4d = SpaceTimeGrid::makeFourDimGrid(clatt, GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi());;
|
||||
|
||||
std::vector<int> seeds4({1,2,3,4});
|
||||
std::vector<int> seeds5({5,6,7,8});
|
||||
std::vector<int> cseeds({5,6,7,8});
|
||||
GridParallelRNG RNG5(FGrid); RNG5.SeedFixedIntegers(seeds5);
|
||||
GridParallelRNG RNG4(UGrid); RNG4.SeedFixedIntegers(seeds4);
|
||||
GridParallelRNG CRNG(Coarse5d);CRNG.SeedFixedIntegers(cseeds);
|
||||
|
||||
Gamma g5(Gamma::Algebra::Gamma5);
|
||||
|
||||
LatticeFermion src(FGrid); gaussian(RNG5,src);// src=src+g5*src;
|
||||
LatticeFermion result(FGrid); result=zero;
|
||||
LatticeFermion ref(FGrid); ref=zero;
|
||||
LatticeFermion tmp(FGrid);
|
||||
LatticeFermion err(FGrid);
|
||||
LatticeGaugeField Umu(UGrid);
|
||||
LatticeGaugeField UmuDD(UGrid);
|
||||
LatticeColourMatrix U(UGrid);
|
||||
LatticeColourMatrix zz(UGrid);
|
||||
|
||||
FieldMetaData header;
|
||||
std::string file("./ckpoint_lat.4000");
|
||||
NerscIO::readConfiguration(Umu,header,file);
|
||||
|
||||
|
||||
if ( params.domaindecompose ) {
|
||||
Lattice<iScalar<vInteger> > coor(UGrid);
|
||||
zz=zero;
|
||||
for(int mu=0;mu<Nd;mu++){
|
||||
LatticeCoordinate(coor,mu);
|
||||
U = PeekIndex<LorentzIndex>(Umu,mu);
|
||||
U = where(mod(coor,params.domainsize)==(Integer)0,zz,U);
|
||||
PokeIndex<LorentzIndex>(UmuDD,U,mu);
|
||||
}
|
||||
} else {
|
||||
UmuDD = Umu;
|
||||
}
|
||||
// SU3::ColdConfiguration(RNG4,Umu);
|
||||
// SU3::TepidConfiguration(RNG4,Umu);
|
||||
// SU3::HotConfiguration(RNG4,Umu);
|
||||
// Umu=zero;
|
||||
|
||||
RealD mass=params.mq;
|
||||
RealD M5=1.8;
|
||||
|
||||
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
|
||||
std::cout<<GridLogMessage << "Building g5R5 hermitian DWF operator" <<std::endl;
|
||||
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
|
||||
DomainWallFermionR Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5);
|
||||
DomainWallFermionR DdwfDD(UmuDD,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5);
|
||||
|
||||
typedef Aggregation<vSpinColourVector,vTComplex,nbasis> Subspace;
|
||||
typedef CoarsenedMatrix<vSpinColourVector,vTComplex,nbasis> CoarseOperator;
|
||||
typedef CoarseOperator::CoarseVector CoarseVector;
|
||||
|
||||
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
|
||||
std::cout<<GridLogMessage << "Calling Aggregation class to build subspace" <<std::endl;
|
||||
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
|
||||
MdagMLinearOperator<DomainWallFermionR,LatticeFermion> HermDefOp(Ddwf);
|
||||
Subspace Aggregates(Coarse5d,FGrid);
|
||||
// Aggregates.CreateSubspace(RNG5,HermDefOp,nbasis);
|
||||
assert ( (nbasis & 0x1)==0);
|
||||
int nb=nbasis/2;
|
||||
std::cout<<GridLogMessage << " nbasis/2 = "<<nb<<std::endl;
|
||||
// Aggregates.CreateSubspace(RNG5,HermDefOp,nb);
|
||||
Aggregates.CreateSubspaceLanczos(RNG5,HermDefOp,nb);
|
||||
for(int n=0;n<nb;n++){
|
||||
G5R5(Aggregates.subspace[n+nb],Aggregates.subspace[n]);
|
||||
std::cout<<GridLogMessage<<n<<" subspace "<<norm2(Aggregates.subspace[n+nb])<<" "<<norm2(Aggregates.subspace[n]) <<std::endl;
|
||||
}
|
||||
for(int n=0;n<nbasis;n++){
|
||||
std::cout<<GridLogMessage << "vec["<<n<<"] = "<<norm2(Aggregates.subspace[n]) <<std::endl;
|
||||
}
|
||||
|
||||
// for(int i=0;i<nbasis;i++){
|
||||
// result = Aggregates.subspace[i];
|
||||
// Aggregates.subspace[i]=result+g5*result;
|
||||
// }
|
||||
result=zero;
|
||||
|
||||
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
|
||||
std::cout<<GridLogMessage << "Building coarse representation of Indef operator" <<std::endl;
|
||||
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
|
||||
Gamma5R5HermitianLinearOperator<DomainWallFermionR,LatticeFermion> HermIndefOp(Ddwf);
|
||||
Gamma5R5HermitianLinearOperator<DomainWallFermionR,LatticeFermion> HermIndefOpDD(DdwfDD);
|
||||
CoarsenedMatrix<vSpinColourVector,vTComplex,nbasis> LDOp(*Coarse5d);
|
||||
LDOp.CoarsenOperator(FGrid,HermIndefOp,Aggregates);
|
||||
|
||||
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
|
||||
std::cout<<GridLogMessage << "Testing some coarse space solvers " <<std::endl;
|
||||
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
|
||||
CoarseVector c_src (Coarse5d);
|
||||
CoarseVector c_res (Coarse5d);
|
||||
gaussian(CRNG,c_src);
|
||||
c_res=zero;
|
||||
|
||||
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
|
||||
std::cout<<GridLogMessage << "Solving posdef-CG on coarse space "<< std::endl;
|
||||
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
|
||||
MdagMLinearOperator<CoarseOperator,CoarseVector> PosdefLdop(LDOp);
|
||||
ConjugateGradient<CoarseVector> CG(1.0e-6,100000);
|
||||
// CG(PosdefLdop,c_src,c_res);
|
||||
|
||||
// std::cout<<GridLogMessage << "**************************************************"<< std::endl;
|
||||
// std::cout<<GridLogMessage << "Solving indef-MCR on coarse space "<< std::endl;
|
||||
// std::cout<<GridLogMessage << "**************************************************"<< std::endl;
|
||||
// HermitianLinearOperator<CoarseOperator,CoarseVector> HermIndefLdop(LDOp);
|
||||
// ConjugateResidual<CoarseVector> MCR(1.0e-6,100000);
|
||||
//MCR(HermIndefLdop,c_src,c_res);
|
||||
|
||||
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
|
||||
std::cout<<GridLogMessage << "Building deflation preconditioner "<< std::endl;
|
||||
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
|
||||
|
||||
MultiGridPreconditioner <vSpinColourVector,vTComplex,nbasis,DomainWallFermionR> Precon (Aggregates, LDOp,
|
||||
HermIndefOp,Ddwf,
|
||||
HermIndefOp,Ddwf);
|
||||
|
||||
MultiGridPreconditioner <vSpinColourVector,vTComplex,nbasis,DomainWallFermionR> PreconDD(Aggregates, LDOp,
|
||||
HermIndefOp,Ddwf,
|
||||
HermIndefOpDD,DdwfDD);
|
||||
TrivialPrecon<LatticeFermion> simple;
|
||||
|
||||
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
|
||||
std::cout<<GridLogMessage << "Testing smoother efficacy"<< std::endl;
|
||||
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
|
||||
// Precon.SmootherTest(src);
|
||||
|
||||
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
|
||||
std::cout<<GridLogMessage << "Testing DD smoother efficacy"<< std::endl;
|
||||
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
|
||||
// PreconDD.SmootherTest(src);
|
||||
|
||||
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
|
||||
std::cout<<GridLogMessage << "Testing SAP smoother efficacy"<< std::endl;
|
||||
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
|
||||
// PreconDD.SAP(src,result);
|
||||
|
||||
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
|
||||
std::cout<<GridLogMessage << "Unprec CG "<< std::endl;
|
||||
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
|
||||
|
||||
// TrivialPrecon<LatticeFermion> simple;
|
||||
// ConjugateGradient<LatticeFermion> fCG(1.0e-8,100000);
|
||||
// fCG(HermDefOp,src,result);
|
||||
// exit(0);
|
||||
|
||||
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
|
||||
std::cout<<GridLogMessage << "Testing GCR on indef matrix "<< std::endl;
|
||||
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
|
||||
// PrecGeneralisedConjugateResidual<LatticeFermion> UPGCR(1.0e-8,100000,simple,8,128);
|
||||
// UPGCR(HermIndefOp,src,result);
|
||||
|
||||
|
||||
/// Get themax eval
|
||||
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
|
||||
std::cout<<GridLogMessage <<" Applying power method to find spectral range "<<std::endl;
|
||||
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
|
||||
Precon.PowerMethod(src);
|
||||
|
||||
|
||||
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
|
||||
std::cout<<GridLogMessage << "Building a two level DDPGCR "<< std::endl;
|
||||
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
|
||||
// PrecGeneralisedConjugateResidual<LatticeFermion> PGCRDD(1.0e-8,100000,PreconDD,8,128);
|
||||
// result=zero;
|
||||
// std::cout<<GridLogMessage<<"checking norm src "<<norm2(src)<<std::endl;
|
||||
// PGCRDD(HermIndefOp,src,result);
|
||||
|
||||
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
|
||||
std::cout<<GridLogMessage << "Building a two level PGCR "<< std::endl;
|
||||
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
|
||||
PrecGeneralisedConjugateResidual<LatticeFermion> PGCR(1.0e-8,100000,Precon,8,8);
|
||||
std::cout<<GridLogMessage<<"checking norm src "<<norm2(src)<<std::endl;
|
||||
result=zero;
|
||||
PGCR(HermIndefOp,src,result);
|
||||
|
||||
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
|
||||
std::cout<<GridLogMessage << "Red Black Prec CG "<< std::endl;
|
||||
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
|
||||
SchurDiagMooeeOperator<DomainWallFermionR,LatticeFermion> HermOpEO(Ddwf);
|
||||
ConjugateGradient<LatticeFermion> pCG(1.0e-8,10000);
|
||||
|
||||
LatticeFermion src_o(FrbGrid);
|
||||
LatticeFermion result_o(FrbGrid);
|
||||
pickCheckerboard(Odd,src_o,src);
|
||||
result_o=zero;
|
||||
|
||||
pCG(HermOpEO,src_o,result_o);
|
||||
|
||||
|
||||
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
|
||||
std::cout<<GridLogMessage << "Done "<< std::endl;
|
||||
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
|
||||
Grid_finalize();
|
||||
}
|
||||
|
||||
#else
|
||||
int main (int argc, char ** argv)
|
||||
{
|
||||
Grid_init(&argc,&argv);
|
||||
|
||||
params.domaindecompose = 1;
|
||||
params.domainsize= 1;
|
||||
params.order = 1;
|
||||
params.Ls = 1;
|
||||
params.mq = 1;
|
||||
params.lo = 1;
|
||||
params.hi = 1;
|
||||
params.steps = 1;
|
||||
|
||||
const int Ls=params.Ls;
|
||||
|
||||
GridCartesian * UGrid = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi());
|
||||
GridRedBlackCartesian * UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
|
||||
|
||||
///////////////////////////////////////////////////
|
||||
// Construct a coarsened grid; utility for this?
|
||||
///////////////////////////////////////////////////
|
||||
std::vector<int> block ({4,4,4,4});
|
||||
const int nbasis= 32;
|
||||
|
||||
std::vector<int> clatt = GridDefaultLatt();
|
||||
for(int d=0;d<clatt.size();d++){
|
||||
clatt[d] = clatt[d]/block[d];
|
||||
}
|
||||
GridCartesian *Coarse4d = SpaceTimeGrid::makeFourDimGrid(clatt, GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi());;
|
||||
|
||||
std::vector<int> seedsFine({1,2,3,4});
|
||||
std::vector<int> seedsCoarse({5,6,7,8});
|
||||
|
||||
GridParallelRNG pRNGFine(UGrid); pRNGFine.SeedFixedIntegers(seedsFine);
|
||||
GridParallelRNG pRNGCoarse(Coarse4d); pRNGCoarse.SeedFixedIntegers(seedsCoarse);
|
||||
|
||||
Gamma g5(Gamma::Algebra::Gamma5);
|
||||
|
||||
LatticeFermion src(UGrid); gaussian(pRNGFine,src);// src=src+g5*src;
|
||||
LatticeFermion result(UGrid); result=zero;
|
||||
LatticeFermion ref(UGrid); ref=zero;
|
||||
LatticeFermion tmp(UGrid);
|
||||
LatticeFermion err(UGrid);
|
||||
LatticeGaugeField Umu(UGrid); SU3::HotConfiguration(pRNGFine,Umu);
|
||||
LatticeGaugeField UmuDD(UGrid);
|
||||
LatticeColourMatrix U(UGrid);
|
||||
LatticeColourMatrix zz(UGrid);
|
||||
|
||||
if ( params.domaindecompose ) {
|
||||
Lattice<iScalar<vInteger> > coor(UGrid);
|
||||
zz=zero;
|
||||
for(int mu=0;mu<Nd;mu++){
|
||||
LatticeCoordinate(coor,mu);
|
||||
U = PeekIndex<LorentzIndex>(Umu,mu);
|
||||
U = where(mod(coor,params.domainsize)==(Integer)0,zz,U);
|
||||
PokeIndex<LorentzIndex>(UmuDD,U,mu);
|
||||
}
|
||||
} else {
|
||||
UmuDD = Umu;
|
||||
}
|
||||
|
||||
RealD mass=params.mq;
|
||||
RealD M5=1.8;
|
||||
|
||||
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
|
||||
std::cout<<GridLogMessage << "Hello "<< std::endl;
|
||||
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
|
||||
|
||||
std::cout << params << std::endl;
|
||||
|
||||
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
|
||||
std::cout<<GridLogMessage << "Building the wilson operator" <<std::endl;
|
||||
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
|
||||
|
||||
WilsonFermionR Dw(Umu,*UGrid,*UrbGrid,mass);
|
||||
WilsonFermionR DwDD(UmuDD,*UGrid,*UrbGrid,mass);
|
||||
|
||||
typedef Aggregation<vSpinColourVector,vTComplex,nbasis> Subspace;
|
||||
typedef CoarsenedMatrix<vSpinColourVector,vTComplex,nbasis> CoarseOperator;
|
||||
typedef CoarseOperator::CoarseVector CoarseVector;
|
||||
|
||||
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
|
||||
std::cout<<GridLogMessage << "Calling Aggregation class to build subspace" <<std::endl;
|
||||
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
|
||||
|
||||
Subspace Aggregates(Coarse4d,UGrid);
|
||||
assert ( (nbasis & 0x1)==0);
|
||||
int nb=nbasis/2;
|
||||
std::cout<<GridLogMessage << " nbasis/2 = "<<nb<<std::endl;
|
||||
|
||||
Aggregates.CreateSubspaceRandom(pRNGFine);
|
||||
|
||||
for(int n=0;n<nb;n++){
|
||||
G5R5(Aggregates.subspace[n+nb],Aggregates.subspace[n]);
|
||||
std::cout<<GridLogMessage<<n<<" subspace "<<norm2(Aggregates.subspace[n+nb])<<" "<<norm2(Aggregates.subspace[n]) <<std::endl;
|
||||
}
|
||||
for(int n=0;n<nbasis;n++){
|
||||
std::cout<<GridLogMessage << "vec["<<n<<"] = "<<norm2(Aggregates.subspace[n]) <<std::endl;
|
||||
}
|
||||
|
||||
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
|
||||
std::cout<<GridLogMessage << "Building coarse representation of Dirac operator" <<std::endl;
|
||||
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
|
||||
|
||||
CoarsenedMatrix<vSpinColourVector,vTComplex,nbasis> LDOp(*Coarse4d);
|
||||
// LDOp.CoarsenOperator(UGrid,Dw,Aggregates); // problem with this line
|
||||
|
||||
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
|
||||
std::cout<<GridLogMessage << "Testing some coarse space solvers " <<std::endl;
|
||||
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
|
||||
|
||||
CoarseVector c_src (Coarse4d);
|
||||
CoarseVector c_res (Coarse4d);
|
||||
gaussian(pRNGCoarse,c_src);
|
||||
c_res=zero;
|
||||
|
||||
Grid_finalize();
|
||||
}
|
||||
#endif
|
||||
Loading…
x
Reference in New Issue
Block a user