diff --git a/Grid/algorithms/LinearOperator.h b/Grid/algorithms/LinearOperator.h index 4a5eeeec..b70e2d6e 100644 --- a/Grid/algorithms/LinearOperator.h +++ b/Grid/algorithms/LinearOperator.h @@ -322,7 +322,7 @@ public: }; template - class HermitianSchurOperatorBase : public LinearOperatorBase + class NonHermitianSchurOperatorBase : public LinearOperatorBase { public: virtual RealD Mpc (const Field& in, Field& out) = 0; @@ -334,16 +334,10 @@ public: no = MpcDag(tmp,out); } virtual void HermOpAndNorm(const Field& in, Field& out, RealD& n1, RealD& n2) { - out.Checkerboard() = in.Checkerboard(); - - Mpc(in, out); - - ComplexD dot = innerProduct(in,out); n1 = real(dot); - n2 = norm2(out); + assert(0); } virtual void HermOp(const Field& in, Field& out) { - RealD n1, n2; - HermOpAndNorm(in, out, n1, n2); + assert(0); } void Op(const Field& in, Field& out) { Mpc(in, out); @@ -361,11 +355,11 @@ public: }; template - class HermitianSchurDiagMooeeOperator : public HermitianSchurOperatorBase + class NonHermitianSchurDiagMooeeOperator : public NonHermitianSchurOperatorBase { public: Matrix& _Mat; - HermitianSchurDiagMooeeOperator(Matrix& Mat): _Mat(Mat){}; + NonHermitianSchurDiagMooeeOperator(Matrix& Mat): _Mat(Mat){}; virtual RealD Mpc(const Field& in, Field& out) { Field tmp(in.Grid()); tmp.Checkerboard() = !in.Checkerboard(); @@ -392,13 +386,13 @@ public: }; template - class HermitianSchurDiagOneOperator : public HermitianSchurOperatorBase + class NonHermitianSchurDiagOneOperator : public NonHermitianSchurOperatorBase { protected: Matrix &_Mat; public: - HermitianSchurDiagOneOperator (Matrix& Mat): _Mat(Mat){}; + NonHermitianSchurDiagOneOperator (Matrix& Mat): _Mat(Mat){}; virtual RealD Mpc(const Field& in, Field& out) { Field tmp(in.Grid()); @@ -422,13 +416,13 @@ public: }; template - class HermitianSchurDiagTwoOperator : public HermitianSchurOperatorBase + class NonHermitianSchurDiagTwoOperator : public NonHermitianSchurOperatorBase { protected: Matrix& _Mat; public: - HermitianSchurDiagTwoOperator(Matrix& Mat): _Mat(Mat){}; + NonHermitianSchurDiagTwoOperator(Matrix& Mat): _Mat(Mat){}; virtual RealD Mpc(const Field& in, Field& out) { Field tmp(in.Grid()); diff --git a/Grid/algorithms/iterative/BiCGSTAB.h b/Grid/algorithms/iterative/BiCGSTAB.h index 3e609241..3a7be1ef 100644 --- a/Grid/algorithms/iterative/BiCGSTAB.h +++ b/Grid/algorithms/iterative/BiCGSTAB.h @@ -62,7 +62,7 @@ class BiCGSTAB : public OperatorFunction conformable(psi, src); RealD cp(0), rho(1), rho_prev(0), alpha(1), beta(0), omega(1); - RealD a(0), bo(0), d(0), b(0), ssq(0), qq(0); + RealD a(0), bo(0), b(0), ssq(0); Field p(src); Field r(src); @@ -79,7 +79,8 @@ class BiCGSTAB : public OperatorFunction RealD guess = norm2(psi); assert(std::isnan(guess) == 0); - Linop.HermOpAndNorm(psi, v, d, b); + Linop.Op(psi, v); + b = norm2(v); r = src - v; rhat = r; @@ -131,7 +132,7 @@ class BiCGSTAB : public OperatorFunction LinalgTimer.Stop(); MatrixTimer.Start(); - Linop.HermOp(p,v); + Linop.Op(p,v); MatrixTimer.Stop(); LinalgTimer.Start(); @@ -155,7 +156,7 @@ class BiCGSTAB : public OperatorFunction LinalgTimer.Stop(); MatrixTimer.Start(); - Linop.HermOp(s,t); + Linop.Op(s,t); MatrixTimer.Stop(); LinalgTimer.Start(); @@ -181,7 +182,7 @@ class BiCGSTAB : public OperatorFunction if(cp <= rsq) { SolverTimer.Stop(); - Linop.HermOpAndNorm(psi, v, d, qq); + Linop.Op(psi, v); p = v - src; RealD srcnorm = sqrt(norm2(src)); diff --git a/Grid/algorithms/iterative/BiCGSTABMixedPrec.h b/Grid/algorithms/iterative/BiCGSTABMixedPrec.h index ebdb1071..d9b2a9d7 100644 --- a/Grid/algorithms/iterative/BiCGSTABMixedPrec.h +++ b/Grid/algorithms/iterative/BiCGSTABMixedPrec.h @@ -105,7 +105,7 @@ class MixedPrecisionBiCGSTAB : public LinearFunction for(outer_iter = 0; outer_iter < MaxOuterIterations; outer_iter++) { // Compute double precision rsd and also new RHS vector. - Linop_d.HermOp(sol_d, tmp_d); + Linop_d.Op(sol_d, tmp_d); RealD norm = axpy_norm(src_d, -1., tmp_d, src_d_in); //src_d is residual vector std::cout << GridLogMessage << "MixedPrecisionBiCGSTAB: Outer iteration " << outer_iter << " residual " << norm << " target " << stop << std::endl; diff --git a/Grid/algorithms/iterative/SchurRedBlack.h b/Grid/algorithms/iterative/SchurRedBlack.h index 92f39cdd..d0b133a3 100644 --- a/Grid/algorithms/iterative/SchurRedBlack.h +++ b/Grid/algorithms/iterative/SchurRedBlack.h @@ -405,14 +405,14 @@ namespace Grid { } }; - template class HermitianSchurRedBlackDiagMooeeSolve : public SchurRedBlackBase + template class NonHermitianSchurRedBlackDiagMooeeSolve : public SchurRedBlackBase { public: typedef CheckerBoardedSparseMatrixBase Matrix; - HermitianSchurRedBlackDiagMooeeSolve(OperatorFunction& HermitianRBSolver, const bool initSubGuess = false, + NonHermitianSchurRedBlackDiagMooeeSolve(OperatorFunction& RBSolver, const bool initSubGuess = false, const bool _solnAsInitGuess = false) - : SchurRedBlackBase(HermitianRBSolver, initSubGuess, _solnAsInitGuess) {}; + : SchurRedBlackBase(RBSolver, initSubGuess, _solnAsInitGuess) {}; ////////////////////////////////////////////////////// // Override RedBlack specialisation @@ -458,14 +458,14 @@ namespace Grid { virtual void RedBlackSolve(Matrix& _Matrix, const Field& src_o, Field& sol_o) { - HermitianSchurDiagMooeeOperator _HermOpEO(_Matrix); - this->_HermitianRBSolver(_HermOpEO, src_o, sol_o); assert(sol_o.Checkerboard() == Odd); + NonHermitianSchurDiagMooeeOperator _OpEO(_Matrix); + this->_HermitianRBSolver(_OpEO, src_o, sol_o); assert(sol_o.Checkerboard() == Odd); } virtual void RedBlackSolve(Matrix& _Matrix, const std::vector& src_o, std::vector& sol_o) { - HermitianSchurDiagMooeeOperator _HermOpEO(_Matrix); - this->_HermitianRBSolver(_HermOpEO, src_o, sol_o); + NonHermitianSchurDiagMooeeOperator _OpEO(_Matrix); + this->_HermitianRBSolver(_OpEO, src_o, sol_o); } }; @@ -547,7 +547,7 @@ namespace Grid { } }; - template class HermitianSchurRedBlackDiagTwoSolve : public SchurRedBlackBase + template class NonHermitianSchurRedBlackDiagTwoSolve : public SchurRedBlackBase { public: typedef CheckerBoardedSparseMatrixBase Matrix; @@ -555,9 +555,9 @@ namespace Grid { ///////////////////////////////////////////////////// // Wrap the usual normal equations Schur trick ///////////////////////////////////////////////////// - HermitianSchurRedBlackDiagTwoSolve(OperatorFunction& HermitianRBSolver, const bool initSubGuess = false, + NonHermitianSchurRedBlackDiagTwoSolve(OperatorFunction& RBSolver, const bool initSubGuess = false, const bool _solnAsInitGuess = false) - : SchurRedBlackBase(HermitianRBSolver, initSubGuess, _solnAsInitGuess) {}; + : SchurRedBlackBase(RBSolver, initSubGuess, _solnAsInitGuess) {}; virtual void RedBlackSource(Matrix& _Matrix, const Field& src, Field& src_e, Field& src_o) { @@ -606,14 +606,14 @@ namespace Grid { virtual void RedBlackSolve(Matrix& _Matrix, const Field& src_o, Field& sol_o) { - HermitianSchurDiagTwoOperator _HermOpEO(_Matrix); - this->_HermitianRBSolver(_HermOpEO, src_o, sol_o); + NonHermitianSchurDiagTwoOperator _OpEO(_Matrix); + this->_HermitianRBSolver(_OpEO, src_o, sol_o); }; virtual void RedBlackSolve(Matrix& _Matrix, const std::vector& src_o, std::vector& sol_o) { - SchurDiagTwoOperator _HermOpEO(_Matrix); - this->_HermitianRBSolver(_HermOpEO, src_o, sol_o); + NonHermitianSchurDiagTwoOperator _OpEO(_Matrix); + this->_HermitianRBSolver(_OpEO, src_o, sol_o); } }; } diff --git a/tests/solver/Test_wilsonclover_bicgstab_prec.cc b/tests/solver/Test_wilsonclover_bicgstab_prec.cc index 708cc6f6..c1905400 100644 --- a/tests/solver/Test_wilsonclover_bicgstab_prec.cc +++ b/tests/solver/Test_wilsonclover_bicgstab_prec.cc @@ -77,7 +77,7 @@ int main (int argc, char ** argv) pickCheckerboard(Odd, src_o, src); result_o = Zero(); - HermitianSchurDiagMooeeOperator HermOp(Dw); + NonHermitianSchurDiagMooeeOperator HermOp(Dw); BiCGSTAB CG(1.0e-8,10000); CG(HermOp, src_o, result_o); diff --git a/tests/solver/Test_wilsonclover_bicgstab_schur.cc b/tests/solver/Test_wilsonclover_bicgstab_schur.cc index c949748f..af021267 100644 --- a/tests/solver/Test_wilsonclover_bicgstab_schur.cc +++ b/tests/solver/Test_wilsonclover_bicgstab_schur.cc @@ -73,7 +73,7 @@ int main (int argc, char ** argv) WilsonCloverFermionR Dw(Umu, Grid, RBGrid, mass, csw_r, csw_t); BiCGSTAB CG(1.0e-8,10000); - HermitianSchurRedBlackDiagMooeeSolve SchurSolver(CG); + NonHermitianSchurRedBlackDiagMooeeSolve SchurSolver(CG); SchurSolver(Dw, src, result); diff --git a/tests/solver/Test_wilsonclover_bicgstab_unprec.cc b/tests/solver/Test_wilsonclover_bicgstab_unprec.cc index 021b3247..21b7521d 100644 --- a/tests/solver/Test_wilsonclover_bicgstab_unprec.cc +++ b/tests/solver/Test_wilsonclover_bicgstab_unprec.cc @@ -72,7 +72,7 @@ int main (int argc, char ** argv) RealD csw_t = 1.0; WilsonCloverFermionR Dw(Umu, Grid, RBGrid, mass, csw_r, csw_t); - HermitianLinearOperator HermOp(Dw); + NonHermitianLinearOperator HermOp(Dw); BiCGSTAB CG(1.0e-8,10000); CG(HermOp,src,result); diff --git a/tests/solver/Test_wilsonclover_mixedbicgstab_prec.cc b/tests/solver/Test_wilsonclover_mixedbicgstab_prec.cc index 848191c5..0af83f8b 100644 --- a/tests/solver/Test_wilsonclover_mixedbicgstab_prec.cc +++ b/tests/solver/Test_wilsonclover_mixedbicgstab_prec.cc @@ -80,16 +80,16 @@ int main (int argc, char ** argv) result_o_2.Checkerboard() = Odd; result_o_2 = Zero(); - HermitianSchurDiagMooeeOperator HermOpEO_d(Dw_d); - HermitianSchurDiagMooeeOperator HermOpEO_f(Dw_f); + NonHermitianSchurDiagMooeeOperator NonHermOpEO_d(Dw_d); + NonHermitianSchurDiagMooeeOperator NonHermOpEO_f(Dw_f); std::cout << GridLogMessage << "::::::::::::: Starting mixed CG" << std::endl; - MixedPrecisionBiCGSTAB mCG(1.0e-8, 10000, 50, FrbGrid_f, HermOpEO_f, HermOpEO_d); + MixedPrecisionBiCGSTAB mCG(1.0e-8, 10000, 50, FrbGrid_f, NonHermOpEO_f, NonHermOpEO_d); mCG(src_o, result_o); std::cout << GridLogMessage << "::::::::::::: Starting regular CG" << std::endl; BiCGSTAB CG(1.0e-8, 10000); - CG(HermOpEO_d, src_o, result_o_2); + CG(NonHermOpEO_d, src_o, result_o_2); LatticeFermionD diff_o(FrbGrid_d); RealD diff = axpy_norm(diff_o, -1.0, result_o, result_o_2);