diff --git a/lib/algorithms/iterative/ImplicitlyRestartedLanczos.h b/lib/algorithms/iterative/ImplicitlyRestartedLanczos.h index 7b85c095..7d5a1889 100644 --- a/lib/algorithms/iterative/ImplicitlyRestartedLanczos.h +++ b/lib/algorithms/iterative/ImplicitlyRestartedLanczos.h @@ -181,6 +181,7 @@ enum IRLdiagonalisation { template class ImplicitlyRestartedLanczosHermOpTester : public ImplicitlyRestartedLanczosTester { public: + LinearFunction &_HermOp; ImplicitlyRestartedLanczosHermOpTester(LinearFunction &HermOp) : _HermOp(HermOp) { }; int ReconstructEval(int j,RealD resid,Field &B, RealD &eval,RealD evalMaxApprox) @@ -243,6 +244,13 @@ class ImplicitlyRestartedLanczos { ///////////////////////// public: + + static void Deflate(const std::vector &_v, + const std::vector& eval, + const Field& src_orig,Field& result) { + basisDeflate(_v,eval,src_orig,result); + } + ////////////////////////////////////////////////////////////////// // PAB: ////////////////////////////////////////////////////////////////// diff --git a/lib/algorithms/iterative/LocalCoherenceLanczos.h b/lib/algorithms/iterative/LocalCoherenceLanczos.h index d5d1bbc2..4c05f4c7 100644 --- a/lib/algorithms/iterative/LocalCoherenceLanczos.h +++ b/lib/algorithms/iterative/LocalCoherenceLanczos.h @@ -28,7 +28,9 @@ Author: paboyle /* END LEGAL */ #ifndef GRID_LOCAL_COHERENCE_IRL_H #define GRID_LOCAL_COHERENCE_IRL_H + namespace Grid { + struct LanczosParams : Serializable { public: GRID_SERIALIZABLE_CLASS_MEMBERS(LanczosParams, @@ -70,21 +72,24 @@ public: typedef Lattice FineField; LinearOperatorBase &_Linop; - Aggregation &_Aggregate; + std::vector &subspace; - ProjectedHermOp(LinearOperatorBase& linop, Aggregation &aggregate) : - _Linop(linop), - _Aggregate(aggregate) { }; + ProjectedHermOp(LinearOperatorBase& linop, std::vector & _subspace) : + _Linop(linop), subspace(_subspace) + { + assert(subspace.size() >0); + }; void operator()(const CoarseField& in, CoarseField& out) { + GridBase *FineGrid = subspace[0]._grid; + int checkerboard = subspace[0].checkerboard; + + FineField fin (FineGrid); fin.checkerboard= checkerboard; + FineField fout(FineGrid); fout.checkerboard = checkerboard; - GridBase *FineGrid = _Aggregate.FineGrid; - FineField fin(FineGrid); - FineField fout(FineGrid); - - _Aggregate.PromoteFromSubspace(in,fin); std::cout< & _poly; LinearOperatorBase &_Linop; - Aggregation &_Aggregate; + std::vector &subspace; - ProjectedFunctionHermOp(OperatorFunction & poly,LinearOperatorBase& linop, - Aggregation &aggregate) : + ProjectedFunctionHermOp(OperatorFunction & poly, + LinearOperatorBase& linop, + std::vector & _subspace) : _poly(poly), _Linop(linop), - _Aggregate(aggregate) { }; + subspace(_subspace) + { }; void operator()(const CoarseField& in, CoarseField& out) { - - GridBase *FineGrid = _Aggregate.FineGrid; - - FineField fin(FineGrid) ;fin.checkerboard =_Aggregate.checkerboard; - FineField fout(FineGrid);fout.checkerboard =_Aggregate.checkerboard; - _Aggregate.PromoteFromSubspace(in,fin); std::cout< & _Poly; OperatorFunction & _smoother; LinearOperatorBase &_Linop; - Aggregation &_Aggregate; - RealD _coarse_relax_tol; + RealD _coarse_relax_tol; + std::vector &_subspace; + ImplicitlyRestartedLanczosSmoothedTester(LinearFunction &Poly, OperatorFunction &smoother, LinearOperatorBase &Linop, - Aggregation &Aggregate, + std::vector &subspace, RealD coarse_relax_tol=5.0e3) - : _smoother(smoother), _Linop(Linop),_Aggregate(Aggregate), _Poly(Poly), _coarse_relax_tol(coarse_relax_tol) { }; + : _smoother(smoother), _Linop(Linop), _Poly(Poly), _subspace(subspace), + _coarse_relax_tol(coarse_relax_tol) + { }; int TestConvergence(int j,RealD eresid,CoarseField &B, RealD &eval,RealD evalMaxApprox) { CoarseField v(B); RealD eval_poly = eval; + // Apply operator _Poly(B,v); @@ -168,14 +180,13 @@ class ImplicitlyRestartedLanczosSmoothedTester : public ImplicitlyRestartedLanc } int ReconstructEval(int j,RealD eresid,CoarseField &B, RealD &eval,RealD evalMaxApprox) { - GridBase *FineGrid = _Aggregate.FineGrid; - - int checkerboard = _Aggregate.checkerboard; - + GridBase *FineGrid = _subspace[0]._grid; + int checkerboard = _subspace[0].checkerboard; FineField fB(FineGrid);fB.checkerboard =checkerboard; FineField fv(FineGrid);fv.checkerboard =checkerboard; - _Aggregate.PromoteFromSubspace(B,fv); + blockPromote(B,fv,_subspace); + _smoother(_Linop,fv,fB); RealD eval_poly = eval; @@ -217,27 +228,80 @@ protected: int _checkerboard; LinearOperatorBase & _FineOp; - // FIXME replace Aggregation with vector of fine; the code reuse is too small for - // the hassle and complexity of cross coupling. - Aggregation _Aggregate; - std::vector evals_fine; - std::vector evals_coarse; - std::vector evec_coarse; + std::vector &evals_fine; + std::vector &evals_coarse; + std::vector &subspace; + std::vector &evec_coarse; + +private: + std::vector _evals_fine; + std::vector _evals_coarse; + std::vector _subspace; + std::vector _evec_coarse; + public: + static void Deflate(std::vector subspace, + std::vector evec_coarse, + std::vector eval_coarse, + const FineField& src_orig,FineField& result) + { + int N = (int)evec_coarse.size(); + CoarseField src_coarse(evec_coarse[0]._grid); + CoarseField res_coarse(evec_coarse[0]._grid); res_coarse = zero; + blockProject(src_orig,src_coarse,subspace); + for (int i=0;i &FineOp, - int checkerboard) : + GridBase *CoarseGrid, + LinearOperatorBase &FineOp, + int checkerboard) : _CoarseGrid(CoarseGrid), _FineGrid(FineGrid), - _Aggregate(CoarseGrid,FineGrid,checkerboard), _FineOp(FineOp), - _checkerboard(checkerboard) + _checkerboard(checkerboard), + evals_fine (_evals_fine), + evals_coarse(_evals_coarse), + subspace (_subspace), + evec_coarse(_evec_coarse) { evals_fine.resize(0); evals_coarse.resize(0); }; - void Orthogonalise(void ) { _Aggregate.Orthogonalise(); } + ////////////////////////////////////////////////////////////////////////// + // Alternate constructore, external storage for use by Hadrons module + ////////////////////////////////////////////////////////////////////////// + LocalCoherenceLanczos(GridBase *FineGrid, + GridBase *CoarseGrid, + LinearOperatorBase &FineOp, + int checkerboard, + std::vector &ext_subspace, + std::vector &ext_coarse, + std::vector &ext_eval_fine, + std::vector &ext_eval_coarse + ) : + _CoarseGrid(CoarseGrid), + _FineGrid(FineGrid), + _FineOp(FineOp), + _checkerboard(checkerboard), + evals_fine (ext_eval_fine), + evals_coarse(ext_eval_coarse), + subspace (ext_subspace), + evec_coarse (ext_coarse) + { + evals_fine.resize(0); + evals_coarse.resize(0); + }; + + void Orthogonalise(void ) { + CoarseScalar InnerProd(_CoarseGrid); + blockOrthogonalise(InnerProd,subspace);std::cout << GridLogMessage <<" Gramm-Schmidt pass 1"< static RealD normalise(T& v) { @@ -246,43 +310,44 @@ public: v = v * (1.0/nn); return nn; } - + /* void fakeFine(void) { int Nk = nbasis; - _Aggregate.subspace.resize(Nk,_FineGrid); - _Aggregate.subspace[0]=1.0; - _Aggregate.subspace[0].checkerboard=_checkerboard; - normalise(_Aggregate.subspace[0]); + subspace.resize(Nk,_FineGrid); + subspace[0]=1.0; + subspace[0].checkerboard=_checkerboard; + normalise(subspace[0]); PlainHermOp Op(_FineOp); for(int k=1;k Op(_FineOp); ImplicitlyRestartedLanczosHermOpTester SimpleTester(Op); for(int k=0;k ChebySmooth(cheby_smooth); - ProjectedFunctionHermOp ChebyOp (ChebySmooth,_FineOp,_Aggregate); - ImplicitlyRestartedLanczosSmoothedTester ChebySmoothTester(ChebyOp,ChebySmooth,_FineOp,_Aggregate,relax); + ProjectedFunctionHermOp ChebyOp (ChebySmooth,_FineOp,_subspace); + ImplicitlyRestartedLanczosSmoothedTester ChebySmoothTester(ChebyOp,ChebySmooth,_FineOp,subspace,relax); for(int k=0;k Op(_FineOp); evals_fine.resize(Nm); - _Aggregate.subspace.resize(Nm,_FineGrid); + subspace.resize(Nm,_FineGrid); ImplicitlyRestartedLanczos IRL(ChebyOp,Op,Nstop,Nk,Nm,resid,MaxIt,betastp,MinRes); FineField src(_FineGrid); src=1.0; src.checkerboard = _checkerboard; int Nconv; - IRL.calc(evals_fine,_Aggregate.subspace,src,Nconv,false); + IRL.calc(evals_fine,subspace,src,Nconv,false); // Shrink down to number saved assert(Nstop>=nbasis); assert(Nconv>=nbasis); evals_fine.resize(nbasis); - _Aggregate.subspace.resize(nbasis,_FineGrid); + subspace.resize(nbasis,_FineGrid); } void calcCoarse(ChebyParams cheby_op,ChebyParams cheby_smooth,RealD relax, int Nstop, int Nk, int Nm,RealD resid, RealD MaxIt, RealD betastp, int MinRes) { Chebyshev Cheby(cheby_op); - ProjectedHermOp Op(_FineOp,_Aggregate); - ProjectedFunctionHermOp ChebyOp (Cheby,_FineOp,_Aggregate); + ProjectedHermOp Op(_FineOp,_subspace); + ProjectedFunctionHermOp ChebyOp (Cheby,_FineOp,_subspace); ////////////////////////////////////////////////////////////////////////////////////////////////// // create a smoother and see if we can get a cheap convergence test and smooth inside the IRL ////////////////////////////////////////////////////////////////////////////////////////////////// Chebyshev ChebySmooth(cheby_smooth); - ImplicitlyRestartedLanczosSmoothedTester ChebySmoothTester(ChebyOp,ChebySmooth,_FineOp,_Aggregate,relax); + ImplicitlyRestartedLanczosSmoothedTester ChebySmoothTester(ChebyOp,ChebySmooth,_FineOp,_subspace,relax); evals_coarse.resize(Nm); evec_coarse.resize(Nm,_CoarseGrid); diff --git a/lib/communicator/Communicator_mpi3.cc b/lib/communicator/Communicator_mpi3.cc index ef47d617..6732dcdf 100644 --- a/lib/communicator/Communicator_mpi3.cc +++ b/lib/communicator/Communicator_mpi3.cc @@ -49,6 +49,7 @@ void CartesianCommunicator::Init(int *argc, char ***argv) Grid_quiesce_nodes(); + // Never clean up as done once. MPI_Comm_dup (MPI_COMM_WORLD,&communicator_world); GlobalSharedMemory::Init(communicator_world); @@ -88,6 +89,8 @@ CartesianCommunicator::CartesianCommunicator(const std::vector &processors) GlobalSharedMemory::OptimalCommunicator (processors,optimal_comm); // Remap using the shared memory optimising routine InitFromMPICommunicator(processors,optimal_comm); SetCommunicator(optimal_comm); + // Free the temp communicator + MPI_Comm_free(&optimal_comm); } ////////////////////////////////// @@ -183,8 +186,8 @@ CartesianCommunicator::CartesianCommunicator(const std::vector &processors, } else { srank = 0; - comm_split = parent.communicator; - // std::cout << " Inherited communicator " < &processors, ////////////////////////////////////////////////////////////////////////////////////////////////////// SetCommunicator(comm_split); + // Free the temp communicator + MPI_Comm_free(&comm_split); + if(0){ std::cout << " ndim " <<_ndimension<<" " << parent._ndimension << std::endl; for(int d=0;d &processors, void CartesianCommunicator::InitFromMPICommunicator(const std::vector &processors, MPI_Comm communicator_base) { + //////////////////////////////////////////////////// + // Creates communicator, and the communicator_halo + //////////////////////////////////////////////////// _ndimension = processors.size(); _processor_coor.resize(_ndimension); diff --git a/lib/qcd/action/fermion/CayleyFermion5D.cc b/lib/qcd/action/fermion/CayleyFermion5D.cc index eace6484..e053b98c 100644 --- a/lib/qcd/action/fermion/CayleyFermion5D.cc +++ b/lib/qcd/action/fermion/CayleyFermion5D.cc @@ -73,7 +73,7 @@ void CayleyFermion5D::DminusDag(const FermionField &psi, FermionField &chi this->DW(psi,tmp_f,DaggerYes); for(int s=0;s subspace(nbasis,FGrid); @@ -119,7 +120,7 @@ int main (int argc, char ** argv) MdagMLinearOperator HermDefOp(Ddwf); typedef Aggregation Subspace; - Subspace Aggregates(Coarse5d,FGrid); + Subspace Aggregates(Coarse5d,FGrid,cb); Aggregates.CreateSubspaceRandom(RNG5); subspace=Aggregates.subspace; diff --git a/tests/debug/Test_cayley_ldop_cr.cc b/tests/debug/Test_cayley_ldop_cr.cc index cbefdd46..c6005fd0 100644 --- a/tests/debug/Test_cayley_ldop_cr.cc +++ b/tests/debug/Test_cayley_ldop_cr.cc @@ -78,6 +78,7 @@ int main (int argc, char ** argv) RealD mass=0.1; RealD M5=1.5; + int cb=0; std::cout< HermDefOp(Ddwf); - Subspace Aggregates(Coarse5d,FGrid); + Subspace Aggregates(Coarse5d,FGrid,cb); Aggregates.CreateSubspace(RNG5,HermDefOp); diff --git a/tests/lanczos/Test_dwf_compressed_lanczos_reorg.cc b/tests/lanczos/Test_dwf_compressed_lanczos_reorg.cc index 4c702a33..3dff4b90 100644 --- a/tests/lanczos/Test_dwf_compressed_lanczos_reorg.cc +++ b/tests/lanczos/Test_dwf_compressed_lanczos_reorg.cc @@ -56,12 +56,12 @@ public: void checkpointFine(std::string evecs_file,std::string evals_file) { - assert(this->_Aggregate.subspace.size()==nbasis); + assert(this->subspace.size()==nbasis); emptyUserRecord record; Grid::QCD::ScidacWriter WR; WR.open(evecs_file); for(int k=0;k_Aggregate.subspace[k],record); + WR.writeScidacFieldRecord(this->subspace[k],record); } WR.close(); @@ -72,7 +72,7 @@ public: void checkpointFineRestore(std::string evecs_file,std::string evals_file) { this->evals_fine.resize(nbasis); - this->_Aggregate.subspace.resize(nbasis,this->_FineGrid); + this->subspace.resize(nbasis,this->_FineGrid); std::cout << GridLogIRL<< "checkpointFineRestore: Reading evals from "<_Aggregate.subspace[k].checkerboard=this->_checkerboard; - RD.readScidacFieldRecord(this->_Aggregate.subspace[k],record); + this->subspace[k].checkerboard=this->_checkerboard; + RD.readScidacFieldRecord(this->subspace[k],record); } RD.close(); @@ -221,7 +221,9 @@ int main (int argc, char ** argv) { std::cout << GridLogIRL<<"Checkpointing Fine evecs"<