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

Added Staggered Type Preconditioned operator

Browse Source
This commit is contained in:
Chulwoo Jung 2017-08-18 01:44:31 -04:00
parent e73e4b4002
commit 0145685f96
2 changed files with 56 additions and 4 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

52
lib/algorithms/LinearOperator.h
View File

@ -300,6 +300,58 @@ namespace Grid {
} }
}; };
template<class Matrix,class Field>
class SchurStagOperator : public LinearOperatorBase<Field> {
protected:
Matrix &_Mat;
public:
SchurStagOperator (Matrix &Mat): _Mat(Mat){};
virtual RealD Mpc (const Field &in, Field &out) {
Field tmp(in._grid);
Field tmp2(in._grid);
_Mat.Mooee(in,out);
_Mat.Mooee(out,tmp);
_Mat.Meooe(in,out);
_Mat.Meooe(out,tmp2);
return axpy_norm(out,-1.0,tmp2,tmp);
}
virtual RealD MpcDag (const Field &in, Field &out){
return Mpc(in,out);
}
#if 0
virtual void MpcDagMpc(const Field &in, Field &out,RealD &ni,RealD &no) {
Field tmp(in._grid);
ni=Mpc(in,tmp);
no=MpcDag(tmp,out);
}
#endif
void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){
n1 = Mpc(in,out);n2=0.;
}
void HermOp(const Field &in, Field &out){
RealD n1,n2;
HermOpAndNorm(in,out,n1,n2);
}
void Op (const Field &in, Field &out){
Mpc(in,out);
}
void AdjOp (const Field &in, Field &out){
MpcDag(in,out);
}
// Support for coarsening to a multigrid
void OpDiag (const Field &in, Field &out) {
assert(0); // must coarsen the unpreconditioned system
}
void OpDir (const Field &in, Field &out,int dir,int disp) {
assert(0);
}
};
#if 0 #if 0
// This is specific to (Z)mobius fermions // This is specific to (Z)mobius fermions
template<class Matrix, class Field> template<class Matrix, class Field>

8
lib/algorithms/iterative/ImplicitlyRestartedLanczosCJ.h
View File

@ -661,12 +661,12 @@ PARALLEL_FOR_LOOP
RealD vv0 = norm2(v); RealD vv0 = norm2(v);
eval2[i] = vnum/vden; eval2[i] = vnum/vden;
v -= eval2[i]*B[i]; v -= eval2[i]*B[i];
RealD vv = norm2(v); RealD vv = norm2(v)/vv0;
std::cout.precision(13); std::cout.precision(13);
std::cout<<GridLogMessage << "[" << std::setw(3)<< std::setiosflags(std::ios_base::right) <<i<<"] "; std::cout<<GridLogMessage << "[" << std::setw(3)<< std::setiosflags(std::ios_base::right) <<i<<"] ";
std::cout<<"eval = "<<std::setw(25)<< std::setiosflags(std::ios_base::left)<< eval2[i]; std::cout<<"eval = "<<std::setw(25)<< std::setiosflags(std::ios_base::left)<< eval2[i];
std::cout<<"|H B[i] - eval[i]B[i]|^2 "<< std::setw(25)<< std::setiosflags(std::ios_base::right)<< vv; std::cout<<"|H B[i] - eval[i]B[i]|^2/|H B[i]|^2 "<< std::setw(25)<< std::setiosflags(std::ios_base::right)<< vv;
std::cout<<" "<< vnum/(sqrt(vden)*sqrt(vv0)) << std::endl; std::cout<<" "<< vnum/(sqrt(vden)*sqrt(vv0)) << std::endl;
// change the criteria as evals are supposed to be sorted, all evals smaller(larger) than Nstop should have converged // change the criteria as evals are supposed to be sorted, all evals smaller(larger) than Nstop should have converged
@ -696,12 +696,12 @@ void FinalCheck( int _Nk,
RealD vv0 = norm2(v); RealD vv0 = norm2(v);
RealD eval2 = vnum/vden; RealD eval2 = vnum/vden;
v -= eval2*evec[j]; v -= eval2*evec[j];
RealD vv = norm2(v); RealD vv = norm2(v)/vv0;
std::cout.precision(13); std::cout.precision(13);
std::cout<<GridLogMessage << "[" << std::setw(3)<< std::setiosflags(std::ios_base::right) <<j<<"] "; std::cout<<GridLogMessage << "[" << std::setw(3)<< std::setiosflags(std::ios_base::right) <<j<<"] ";
std::cout<<"eval = "<<std::setw(25)<< std::setiosflags(std::ios_base::left)<< eval2; std::cout<<"eval = "<<std::setw(25)<< std::setiosflags(std::ios_base::left)<< eval2;
std::cout<<"|H B[i] - eval[i]B[i]|^2 "<< std::setw(25)<< std::setiosflags(std::ios_base::right)<< vv; std::cout<<"|H B[i] - eval[i]B[i]|^2/|H B[i]|^2 "<< std::setw(25)<< std::setiosflags(std::ios_base::right)<< vv;
std::cout<<" "<< vnum/(sqrt(vden)*sqrt(vv0)) << std::endl; std::cout<<" "<< vnum/(sqrt(vden)*sqrt(vv0)) << std::endl;
eval[j] = eval2; eval[j] = eval2;