diff --git a/TODO b/TODO index c37cbf8b..83bfda5e 100644 --- a/TODO +++ b/TODO @@ -3,19 +3,19 @@ TODO: Large item work list: -1)- BG/Q port and check +1)- BG/Q port and check ; Andrew says ok. 2)- Christoph's local basis expansion Lanczos -3)- Precision conversion and sort out localConvert <-- partial - - - Consistent linear solver flop count/rate -- PARTIAL, time but no flop/s yet +-- +3a)- RNG I/O in ILDG/SciDAC (minor) +3b)- Precision conversion and sort out localConvert <-- partial/easy +3c)- Consistent linear solver flop count/rate -- PARTIAL, time but no flop/s yet 4)- Physical propagator interface 5)- Conserved currents 6)- Multigrid Wilson and DWF, compare to other Multigrid implementations 7)- HDCR resume Recent DONE - --- MultiRHS with spread out extra dim -- Go through filesystem with SciDAC I/O. <--- DONE +-- MultiRHS with spread out extra dim -- Go through filesystem with SciDAC I/O ; <-- DONE ; bmark cori -- Lanczos Remove DenseVector, DenseMatrix; Use Eigen instead. <-- DONE -- GaugeFix into central location <-- DONE -- Scidac and Ildg metadata handling <-- DONE diff --git a/benchmarks/Benchmark_staggered.cc b/benchmarks/Benchmark_staggered.cc index dc2dcf91..f5325b28 100644 --- a/benchmarks/Benchmark_staggered.cc +++ b/benchmarks/Benchmark_staggered.cc @@ -40,7 +40,7 @@ int main (int argc, char ** argv) std::vector simd_layout = GridDefaultSimd(Nd,vComplex::Nsimd()); std::vector mpi_layout = GridDefaultMpi(); GridCartesian Grid(latt_size,simd_layout,mpi_layout); - GridRedBlackCartesian RBGrid(latt_size,simd_layout,mpi_layout); + GridRedBlackCartesian RBGrid(&Grid); int threads = GridThread::GetThreads(); std::cout< simd_layout = GridDefaultSimd(Nd,vComplex::Nsimd()); std::vector mpi_layout = GridDefaultMpi(); GridCartesian Grid(latt_size,simd_layout,mpi_layout); - GridRedBlackCartesian RBGrid(latt_size,simd_layout,mpi_layout); + GridRedBlackCartesian RBGrid(&Grid); int threads = GridThread::GetThreads(); std::cout< #include #include #include +#include #include - -// Lanczos support -//#include #include #include #include -// Eigen/lanczos // EigCg -// MCR // Pcg -// Multishift CG // Hdcg // GCR // etc.. -// integrator/Leapfrog -// integrator/Omelyan -// integrator/ForceGradient - -// montecarlo/hmc -// montecarlo/rhmc -// montecarlo/metropolis -// etc... - - #endif diff --git a/lib/algorithms/iterative/ConjugateGradient.h b/lib/algorithms/iterative/ConjugateGradient.h index ed453161..5c968e04 100644 --- a/lib/algorithms/iterative/ConjugateGradient.h +++ b/lib/algorithms/iterative/ConjugateGradient.h @@ -52,8 +52,8 @@ class ConjugateGradient : public OperatorFunction { MaxIterations(maxit), ErrorOnNoConverge(err_on_no_conv){}; - void operator()(LinearOperatorBase &Linop, const Field &src, - Field &psi) { + void operator()(LinearOperatorBase &Linop, const Field &src, Field &psi) { + psi.checkerboard = src.checkerboard; conformable(psi, src); diff --git a/lib/cartesian/Cartesian_base.h b/lib/cartesian/Cartesian_base.h index f4f9a269..324772c5 100644 --- a/lib/cartesian/Cartesian_base.h +++ b/lib/cartesian/Cartesian_base.h @@ -49,6 +49,8 @@ public: template friend class Lattice; GridBase(const std::vector & processor_grid) : CartesianCommunicator(processor_grid) {}; + GridBase(const std::vector & processor_grid, + const CartesianCommunicator &parent) : CartesianCommunicator(processor_grid,parent) {}; // Physics Grid information. std::vector _simd_layout;// Which dimensions get relayed out over simd lanes. @@ -210,9 +212,6 @@ public: assert(lidx & gcoor,int & gidx){ gidx=0; int mult=1; diff --git a/lib/cartesian/Cartesian_full.h b/lib/cartesian/Cartesian_full.h index 815e3b22..a6a85ab7 100644 --- a/lib/cartesian/Cartesian_full.h +++ b/lib/cartesian/Cartesian_full.h @@ -61,9 +61,29 @@ public: virtual int CheckerBoardShift(int source_cb,int dim,int shift, int osite){ return shift; } + ///////////////////////////////////////////////////////////////////////// + // Constructor takes a parent grid and possibly subdivides communicator. + ///////////////////////////////////////////////////////////////////////// GridCartesian(const std::vector &dimensions, - const std::vector &simd_layout, - const std::vector &processor_grid) : GridBase(processor_grid) + const std::vector &simd_layout, + const std::vector &processor_grid, + const GridCartesian &parent) : GridBase(processor_grid,parent) + { + Init(dimensions,simd_layout,processor_grid); + } + ///////////////////////////////////////////////////////////////////////// + // Construct from comm world + ///////////////////////////////////////////////////////////////////////// + GridCartesian(const std::vector &dimensions, + const std::vector &simd_layout, + const std::vector &processor_grid) : GridBase(processor_grid) + { + Init(dimensions,simd_layout,processor_grid); + } + + void Init(const std::vector &dimensions, + const std::vector &simd_layout, + const std::vector &processor_grid) { /////////////////////// // Grid information diff --git a/lib/cartesian/Cartesian_red_black.h b/lib/cartesian/Cartesian_red_black.h index b1a5b9ef..f89cacc5 100644 --- a/lib/cartesian/Cartesian_red_black.h +++ b/lib/cartesian/Cartesian_red_black.h @@ -112,24 +112,57 @@ public: } }; - GridRedBlackCartesian(const GridBase *base) : GridRedBlackCartesian(base->_fdimensions,base->_simd_layout,base->_processors) {}; + //////////////////////////////////////////////////////////// + // Create Redblack from original grid; require full grid pointer ? + //////////////////////////////////////////////////////////// + GridRedBlackCartesian(const GridBase *base) : GridBase(base->_processors,*base) + { + int dims = base->_ndimension; + std::vector checker_dim_mask(dims,1); + int checker_dim = 0; + Init(base->_fdimensions,base->_simd_layout,base->_processors,checker_dim_mask,checker_dim); + }; - GridRedBlackCartesian(const std::vector &dimensions, + //////////////////////////////////////////////////////////// + // Create redblack from original grid, with non-trivial checker dim mask + //////////////////////////////////////////////////////////// + GridRedBlackCartesian(const GridBase *base, + const std::vector &checker_dim_mask, + int checker_dim + ) : GridBase(base->_processors,*base) + { + Init(base->_fdimensions,base->_simd_layout,base->_processors,checker_dim_mask,checker_dim) ; + } +#if 0 + //////////////////////////////////////////////////////////// + // Create redblack grid ;; deprecate these. Should not + // need direct creation of redblack without a full grid to base on + //////////////////////////////////////////////////////////// + GridRedBlackCartesian(const GridBase *base, + const std::vector &dimensions, const std::vector &simd_layout, const std::vector &processor_grid, const std::vector &checker_dim_mask, int checker_dim - ) : GridBase(processor_grid) + ) : GridBase(processor_grid,*base) { Init(dimensions,simd_layout,processor_grid,checker_dim_mask,checker_dim); } - GridRedBlackCartesian(const std::vector &dimensions, + + //////////////////////////////////////////////////////////// + // Create redblack grid + //////////////////////////////////////////////////////////// + GridRedBlackCartesian(const GridBase *base, + const std::vector &dimensions, const std::vector &simd_layout, - const std::vector &processor_grid) : GridBase(processor_grid) + const std::vector &processor_grid) : GridBase(processor_grid,*base) { std::vector checker_dim_mask(dimensions.size(),1); - Init(dimensions,simd_layout,processor_grid,checker_dim_mask,0); + int checker_dim = 0; + Init(dimensions,simd_layout,processor_grid,checker_dim_mask,checker_dim); } +#endif + void Init(const std::vector &dimensions, const std::vector &simd_layout, const std::vector &processor_grid, diff --git a/lib/communicator/Communicator_base.cc b/lib/communicator/Communicator_base.cc index 20c310c0..bcf429ab 100644 --- a/lib/communicator/Communicator_base.cc +++ b/lib/communicator/Communicator_base.cc @@ -67,7 +67,7 @@ void CartesianCommunicator::ShmBufferFreeAll(void) { ///////////////////////////////// // Grid information queries ///////////////////////////////// -int CartesianCommunicator::Dimensions(void) { return _ndimension; }; +int CartesianCommunicator::Dimensions(void) { return _ndimension; }; int CartesianCommunicator::IsBoss(void) { return _processor==0; }; int CartesianCommunicator::BossRank(void) { return 0; }; int CartesianCommunicator::ThisRank(void) { return _processor; }; @@ -96,6 +96,105 @@ void CartesianCommunicator::GlobalSumVector(ComplexD *c,int N) GlobalSumVector((double *)c,2*N); } + +#if defined( GRID_COMMS_MPI) || defined (GRID_COMMS_MPIT) + +CartesianCommunicator::CartesianCommunicator(const std::vector &processors,const CartesianCommunicator &parent) +{ + _ndimension = processors.size(); + assert(_ndimension = parent._ndimension); + + ////////////////////////////////////////////////////////////////////////////////////////////////////// + // split the communicator + ////////////////////////////////////////////////////////////////////////////////////////////////////// + int Nparent; + MPI_Comm_size(parent.communicator,&Nparent); + + int childsize=1; + for(int d=0;d 1 ) { + + std::cout << GridLogMessage<<"Child communicator of "<< std::hex << parent.communicator << std::dec< &processors, MPI_Comm communicator_base) +{ + // if ( communicator_base != communicator_world ) { + // std::cout << "Cartesian communicator created with a non-world communicator"< periodic(_ndimension,1); + MPI_Cart_create(communicator_base, _ndimension,&_processors[0],&periodic[0],1,&communicator); + MPI_Comm_rank(communicator,&_processor); + MPI_Cart_coords(communicator,_processor,_ndimension,&_processor_coor[0]); + + int Size; + MPI_Comm_size(communicator,&Size); + +#ifdef GRID_COMMS_MPIT + communicator_halo.resize (2*_ndimension); + for(int i=0;i<_ndimension*2;i++){ + MPI_Comm_dup(communicator,&communicator_halo[i]); + } +#endif + + assert(Size==_Nprocessors); +} + +CartesianCommunicator::CartesianCommunicator(const std::vector &processors) +{ + InitFromMPICommunicator(processors,communicator_world); +} + +#endif + #if !defined( GRID_COMMS_MPI3) int CartesianCommunicator::NodeCount(void) { return ProcessorCount();}; @@ -147,8 +246,13 @@ void *CartesianCommunicator::ShmBufferTranslate(int rank,void * local_p) { } void CartesianCommunicator::ShmInitGeneric(void){ #if 1 - - int mmap_flag = MAP_SHARED | MAP_ANONYMOUS; + int mmap_flag =0; +#ifdef MAP_ANONYMOUS + mmap_flag = mmap_flag| MAP_SHARED | MAP_ANONYMOUS; +#endif +#ifdef MAP_ANON + mmap_flag = mmap_flag| MAP_SHARED | MAP_ANON; +#endif #ifdef MAP_HUGETLB if ( Hugepages ) mmap_flag |= MAP_HUGETLB; #endif diff --git a/lib/communicator/Communicator_base.h b/lib/communicator/Communicator_base.h index ac866ced..bfdb0da1 100644 --- a/lib/communicator/Communicator_base.h +++ b/lib/communicator/Communicator_base.h @@ -83,6 +83,7 @@ class CartesianCommunicator { std::vector communicator_halo; typedef MPI_Request CommsRequest_t; + #else typedef int CommsRequest_t; #endif @@ -147,11 +148,23 @@ class CartesianCommunicator { // Must call in Grid startup //////////////////////////////////////////////// static void Init(int *argc, char ***argv); - + //////////////////////////////////////////////// - // Constructor of any given grid + // Constructors to sub-divide a parent communicator + // and default to comm world //////////////////////////////////////////////// + CartesianCommunicator(const std::vector &processors,const CartesianCommunicator &parent); CartesianCommunicator(const std::vector &pdimensions_in); + + private: +#if defined (GRID_COMMS_MPI) || defined (GRID_COMMS_MPIT) + //////////////////////////////////////////////// + // Private initialise from an MPI communicator + // Can use after an MPI_Comm_split, but hidden from user so private + //////////////////////////////////////////////// + void InitFromMPICommunicator(const std::vector &processors, MPI_Comm communicator_base); +#endif + public: //////////////////////////////////////////////////////////////////////////////////////// // Wraps MPI_Cart routines, or implements equivalent on other impls diff --git a/lib/communicator/Communicator_mpi.cc b/lib/communicator/Communicator_mpi.cc index bd2a62fb..a55c0164 100644 --- a/lib/communicator/Communicator_mpi.cc +++ b/lib/communicator/Communicator_mpi.cc @@ -52,29 +52,6 @@ void CartesianCommunicator::Init(int *argc, char ***argv) { MPI_Comm_dup (MPI_COMM_WORLD,&communicator_world); ShmInitGeneric(); } - -CartesianCommunicator::CartesianCommunicator(const std::vector &processors) -{ - _ndimension = processors.size(); - std::vector periodic(_ndimension,1); - - _Nprocessors=1; - _processors = processors; - _processor_coor.resize(_ndimension); - - MPI_Cart_create(communicator_world, _ndimension,&_processors[0],&periodic[0],1,&communicator); - MPI_Comm_rank(communicator,&_processor); - MPI_Cart_coords(communicator,_processor,_ndimension,&_processor_coor[0]); - - for(int i=0;i<_ndimension;i++){ - _Nprocessors*=_processors[i]; - } - - int Size; - MPI_Comm_size(communicator,&Size); - - assert(Size==_Nprocessors); -} void CartesianCommunicator::GlobalSum(uint32_t &u){ int ierr=MPI_Allreduce(MPI_IN_PLACE,&u,1,MPI_UINT32_T,MPI_SUM,communicator); assert(ierr==0); diff --git a/lib/communicator/Communicator_mpi3.cc b/lib/communicator/Communicator_mpi3.cc index 44aa1024..dce9588a 100644 --- a/lib/communicator/Communicator_mpi3.cc +++ b/lib/communicator/Communicator_mpi3.cc @@ -450,6 +450,15 @@ void CartesianCommunicator::ProcessorCoorFromRank(int rank, std::vector &c assert(lr!=-1); Lexicographic::CoorFromIndex(coor,lr,_processors); } + +////////////////////////////////// +// Try to subdivide communicator +////////////////////////////////// +CartesianCommunicator::CartesianCommunicator(const std::vector &processors,const CartesianCommunicator &parent) + : CartesianCommunicator(processors) +{ + std::cout << "Attempts to split MPI3 communicators will fail until implemented" < &processors) { int ierr; diff --git a/lib/communicator/Communicator_mpit.cc b/lib/communicator/Communicator_mpit.cc index eb6ef87d..5137c27b 100644 --- a/lib/communicator/Communicator_mpit.cc +++ b/lib/communicator/Communicator_mpit.cc @@ -53,33 +53,6 @@ void CartesianCommunicator::Init(int *argc, char ***argv) { ShmInitGeneric(); } -CartesianCommunicator::CartesianCommunicator(const std::vector &processors) -{ - _ndimension = processors.size(); - std::vector periodic(_ndimension,1); - - _Nprocessors=1; - _processors = processors; - _processor_coor.resize(_ndimension); - - MPI_Cart_create(communicator_world, _ndimension,&_processors[0],&periodic[0],1,&communicator); - MPI_Comm_rank(communicator,&_processor); - MPI_Cart_coords(communicator,_processor,_ndimension,&_processor_coor[0]); - - for(int i=0;i<_ndimension;i++){ - _Nprocessors*=_processors[i]; - } - - communicator_halo.resize (2*_ndimension); - for(int i=0;i<_ndimension*2;i++){ - MPI_Comm_dup(communicator,&communicator_halo[i]); - } - - int Size; - MPI_Comm_size(communicator,&Size); - - assert(Size==_Nprocessors); -} void CartesianCommunicator::GlobalSum(uint32_t &u){ int ierr=MPI_Allreduce(MPI_IN_PLACE,&u,1,MPI_UINT32_T,MPI_SUM,communicator); assert(ierr==0); diff --git a/lib/communicator/Communicator_none.cc b/lib/communicator/Communicator_none.cc index 5319ab93..40feefec 100644 --- a/lib/communicator/Communicator_none.cc +++ b/lib/communicator/Communicator_none.cc @@ -38,6 +38,9 @@ void CartesianCommunicator::Init(int *argc, char *** arv) ShmInitGeneric(); } +CartesianCommunicator::CartesianCommunicator(const std::vector &processors,const CartesianCommunicator &parent) + : CartesianCommunicator(processors) {} + CartesianCommunicator::CartesianCommunicator(const std::vector &processors) { _processors = processors; diff --git a/lib/communicator/Communicator_shmem.cc b/lib/communicator/Communicator_shmem.cc index 3c76c808..ed49285d 100644 --- a/lib/communicator/Communicator_shmem.cc +++ b/lib/communicator/Communicator_shmem.cc @@ -75,6 +75,11 @@ void CartesianCommunicator::Init(int *argc, char ***argv) { ShmInitGeneric(); } +CartesianCommunicator::CartesianCommunicator(const std::vector &processors,const CartesianCommunicator &parent) + : CartesianCommunicator(processors) +{ + std::cout << "Attempts to split SHMEM communicators will fail " < &processors) { _ndimension = processors.size(); diff --git a/lib/lattice/Lattice_reduction.h b/lib/lattice/Lattice_reduction.h index db012c8c..8a3fbece 100644 --- a/lib/lattice/Lattice_reduction.h +++ b/lib/lattice/Lattice_reduction.h @@ -544,7 +544,6 @@ static void sliceInnerProductMatrix( Eigen::MatrixXcd &mat, const Lattice for(int i=0;i(); - // std::cout << " Lorentz N/S/V/M : " << _LorentzN<<" "<<_LorentzScalar<<"/"<<_LorentzVector<<"/"<<_LorentzMatrix<_gsites; + + // std::cout << "R sizeof(sobj)= " <_gsites< munge; - BinaryIO::readLatticeObject< sobj, sobj >(field, filename, munge, offset, format,nersc_csum,scidac_csuma,scidac_csumb); + BinaryIO::readLatticeObject< vobj, sobj >(field, filename, munge, offset, format,nersc_csum,scidac_csuma,scidac_csumb); ///////////////////////////////////////////// // Insist checksum is next record ///////////////////////////////////////////// - readLimeObject(scidacChecksum_,std::string("scidacChecksum"),record_name); + readLimeObject(scidacChecksum_,std::string("scidacChecksum"),std::string(SCIDAC_CHECKSUM)); ///////////////////////////////////////////// // Verify checksums @@ -242,11 +252,19 @@ class GridLimeReader : public BinaryIO { // should this be a do while; can we miss a first record?? while ( limeReaderNextRecord(LimeR) == LIME_SUCCESS ) { + // std::cout << GridLogMessage<< " readLimeObject seeking "<< record_name <<" found record :" < xmlc(nbytes+1,'\0'); limeReaderReadData((void *)&xmlc[0], &nbytes, LimeR); + + // std::cout << GridLogMessage<< " readLimeObject matches XML " << &xmlc[0] <(record_name.c_str()), nbytes); + assert(h!= NULL); err=limeWriteRecordHeader(h, LimeW); assert(err>=0); err=limeWriteRecordData(&xmlstring[0], &nbytes, LimeW); assert(err>=0); err=limeWriterCloseRecord(LimeW); assert(err>=0); limeDestroyHeader(h); + // std::cout << " File offset is now"<_gsites; createLimeRecordHeader(record_name, 0, 0, PayloadSize); + + // std::cout << "W sizeof(sobj)" <_gsites<(); BinarySimpleMunger munge; BinaryIO::writeLatticeObject(field, filename, munge, offset, format,nersc_csum,scidac_csuma,scidac_csumb); @@ -354,7 +383,7 @@ class GridLimeWriter : public BinaryIO { checksum.suma= streama.str(); checksum.sumb= streamb.str(); std::cout << GridLogMessage<<" writing scidac checksums "< + template void writeScidacFieldRecord(Lattice &field,userRecord _userRecord) { - typedef typename vobj::scalar_object sobj; - uint64_t nbytes; GridBase * grid = field._grid; //////////////////////////////////////// @@ -397,6 +424,66 @@ class ScidacWriter : public GridLimeWriter { } }; + +class ScidacReader : public GridLimeReader { + public: + + template + void readScidacFileRecord(GridBase *grid,SerialisableUserFile &_userFile) + { + scidacFile _scidacFile(grid); + readLimeObject(_scidacFile,_scidacFile.SerialisableClassName(),std::string(SCIDAC_PRIVATE_FILE_XML)); + readLimeObject(_userFile,_userFile.SerialisableClassName(),std::string(SCIDAC_FILE_XML)); + } + //////////////////////////////////////////////// + // Write generic lattice field in scidac format + //////////////////////////////////////////////// + template + void readScidacFieldRecord(Lattice &field,userRecord &_userRecord) + { + typedef typename vobj::scalar_object sobj; + GridBase * grid = field._grid; + + //////////////////////////////////////// + // fill the Grid header + //////////////////////////////////////// + FieldMetaData header; + scidacRecord _scidacRecord; + scidacFile _scidacFile; + + ////////////////////////////////////////////// + // Fill the Lime file record by record + ////////////////////////////////////////////// + readLimeObject(header ,std::string("FieldMetaData"),std::string(GRID_FORMAT)); // Open message + readLimeObject(_userRecord,_userRecord.SerialisableClassName(),std::string(SCIDAC_RECORD_XML)); + readLimeObject(_scidacRecord,_scidacRecord.SerialisableClassName(),std::string(SCIDAC_PRIVATE_RECORD_XML)); + readLimeLatticeBinaryObject(field,std::string(ILDG_BINARY_DATA)); + } + void skipPastBinaryRecord(void) { + std::string rec_name(ILDG_BINARY_DATA); + while ( limeReaderNextRecord(LimeR) == LIME_SUCCESS ) { + if ( !strncmp(limeReaderType(LimeR), rec_name.c_str(),strlen(rec_name.c_str()) ) ) { + skipPastObjectRecord(std::string(SCIDAC_CHECKSUM)); + return; + } + } + } + void skipPastObjectRecord(std::string rec_name) { + while ( limeReaderNextRecord(LimeR) == LIME_SUCCESS ) { + if ( !strncmp(limeReaderType(LimeR), rec_name.c_str(),strlen(rec_name.c_str()) ) ) { + return; + } + } + } + void skipScidacFieldRecord() { + skipPastObjectRecord(std::string(GRID_FORMAT)); + skipPastObjectRecord(std::string(SCIDAC_RECORD_XML)); + skipPastObjectRecord(std::string(SCIDAC_PRIVATE_RECORD_XML)); + skipPastBinaryRecord(); + } +}; + + class IldgWriter : public ScidacWriter { public: @@ -425,8 +512,6 @@ class IldgWriter : public ScidacWriter { typedef iLorentzColourMatrix vobj; typedef typename vobj::scalar_object sobj; - uint64_t nbytes; - //////////////////////////////////////// // fill the Grid header //////////////////////////////////////// diff --git a/lib/parallelIO/IldgIOtypes.h b/lib/parallelIO/IldgIOtypes.h index c3a5321c..5b397e14 100644 --- a/lib/parallelIO/IldgIOtypes.h +++ b/lib/parallelIO/IldgIOtypes.h @@ -64,6 +64,11 @@ namespace Grid { // file compatability, so should be correct to assume the undocumented but defacto file structure. ///////////////////////////////////////////////////////////////////////////////// +struct emptyUserRecord : Serializable { + GRID_SERIALIZABLE_CLASS_MEMBERS(emptyUserRecord,int,dummy); + emptyUserRecord() { dummy=0; }; +}; + //////////////////////// // Scidac private file xml // 1.1416 16 16 32 0 diff --git a/lib/parallelIO/MetaData.h b/lib/parallelIO/MetaData.h index 6d45d0a5..ccc8b18f 100644 --- a/lib/parallelIO/MetaData.h +++ b/lib/parallelIO/MetaData.h @@ -85,6 +85,9 @@ namespace Grid { nd=4; dimension.resize(4); boundary.resize(4); + scidac_checksuma=0; + scidac_checksumb=0; + checksum=0; } }; @@ -104,6 +107,7 @@ namespace Grid { header.nd = nd; header.dimension.resize(nd); header.boundary.resize(nd); + header.data_start = 0; for(int d=0;d_fdimensions[d]; } diff --git a/lib/qcd/action/fermion/CayleyFermion5D.cc b/lib/qcd/action/fermion/CayleyFermion5D.cc index 838b1c3d..eace6484 100644 --- a/lib/qcd/action/fermion/CayleyFermion5D.cc +++ b/lib/qcd/action/fermion/CayleyFermion5D.cc @@ -77,7 +77,6 @@ void CayleyFermion5D::DminusDag(const FermionField &psi, FermionField &chi } } - template void CayleyFermion5D::CayleyReport(void) { this->Report(); @@ -119,7 +118,6 @@ template void CayleyFermion5D::CayleyZeroCounters(void) MooeeInvTime=0; } - template void CayleyFermion5D::M5D (const FermionField &psi, FermionField &chi) { diff --git a/lib/qcd/utils/SpaceTimeGrid.cc b/lib/qcd/utils/SpaceTimeGrid.cc index 3ada4a3b..b2b5d9c8 100644 --- a/lib/qcd/utils/SpaceTimeGrid.cc +++ b/lib/qcd/utils/SpaceTimeGrid.cc @@ -60,7 +60,7 @@ GridCartesian *SpaceTimeGrid::makeFiveDimGrid(int Ls,const GridCartesian simd5.push_back(FourDimGrid->_simd_layout[d]); mpi5.push_back(FourDimGrid->_processors[d]); } - return new GridCartesian(latt5,simd5,mpi5); + return new GridCartesian(latt5,simd5,mpi5,*FourDimGrid); } @@ -68,18 +68,14 @@ GridRedBlackCartesian *SpaceTimeGrid::makeFiveDimRedBlackGrid(int Ls,const GridC { int N4=FourDimGrid->_ndimension; int cbd=1; - std::vector latt5(1,Ls); - std::vector simd5(1,1); - std::vector mpi5(1,1); std::vector cb5(1,0); - for(int d=0;d_fdimensions[d]); - simd5.push_back(FourDimGrid->_simd_layout[d]); - mpi5.push_back(FourDimGrid->_processors[d]); cb5.push_back( 1); - } - return new GridRedBlackCartesian(latt5,simd5,mpi5,cb5,cbd); + } + GridCartesian *tmp = makeFiveDimGrid(Ls,FourDimGrid); + GridRedBlackCartesian *ret = new GridRedBlackCartesian(tmp,cb5,cbd); + delete tmp; + return ret; } @@ -97,26 +93,24 @@ GridCartesian *SpaceTimeGrid::makeFiveDimDWFGrid(int Ls,const GridCartes simd5.push_back(1); mpi5.push_back(FourDimGrid->_processors[d]); } - return new GridCartesian(latt5,simd5,mpi5); + return new GridCartesian(latt5,simd5,mpi5,*FourDimGrid); } - +/////////////////////////////////////////////////// +// Interface is inefficient and forces the deletion +// Pass in the non-redblack grid +/////////////////////////////////////////////////// GridRedBlackCartesian *SpaceTimeGrid::makeFiveDimDWFRedBlackGrid(int Ls,const GridCartesian *FourDimGrid) { int N4=FourDimGrid->_ndimension; - int nsimd = FourDimGrid->Nsimd(); int cbd=1; - std::vector latt5(1,Ls); - std::vector simd5(1,nsimd); - std::vector mpi5(1,1); std::vector cb5(1,0); - for(int d=0;d_fdimensions[d]); - simd5.push_back(1); - mpi5.push_back(FourDimGrid->_processors[d]); cb5.push_back(1); - } - return new GridRedBlackCartesian(latt5,simd5,mpi5,cb5,cbd); + } + GridCartesian *tmp = makeFiveDimDWFGrid(Ls,FourDimGrid); + GridRedBlackCartesian *ret = new GridRedBlackCartesian(tmp,cb5,cbd); + delete tmp; + return ret; } diff --git a/lib/util/Lexicographic.h b/lib/util/Lexicographic.h index 2d4e5df5..b922dba5 100644 --- a/lib/util/Lexicographic.h +++ b/lib/util/Lexicographic.h @@ -7,7 +7,7 @@ namespace Grid{ class Lexicographic { public: - static inline void CoorFromIndex (std::vector& coor,int index,std::vector &dims){ + static inline void CoorFromIndex (std::vector& coor,int index,const std::vector &dims){ int nd= dims.size(); coor.resize(nd); for(int d=0;d& coor,int &index,std::vector &dims){ + static inline void IndexFromCoor (const std::vector& coor,int &index,const std::vector &dims){ int nd=dims.size(); int stride=1; index=0; diff --git a/tests/Test_stencil.cc b/tests/Test_stencil.cc index fa4b0b57..e9517446 100644 --- a/tests/Test_stencil.cc +++ b/tests/Test_stencil.cc @@ -48,7 +48,7 @@ int main(int argc, char ** argv) { double volume = latt_size[0]*latt_size[1]*latt_size[2]*latt_size[3]; GridCartesian Fine(latt_size,simd_layout,mpi_layout); - GridRedBlackCartesian rbFine(latt_size,simd_layout,mpi_layout); + GridRedBlackCartesian rbFine(&Fine); GridParallelRNG fRNG(&Fine); // fRNG.SeedFixedIntegers(std::vector({45,12,81,9}); diff --git a/tests/core/Test_cshift_red_black.cc b/tests/core/Test_cshift_red_black.cc index f9269709..c7b0c2f1 100644 --- a/tests/core/Test_cshift_red_black.cc +++ b/tests/core/Test_cshift_red_black.cc @@ -47,7 +47,7 @@ int main (int argc, char ** argv) mask[0]=0; GridCartesian Fine (latt_size,simd_layout,mpi_layout); - GridRedBlackCartesian RBFine(latt_size,simd_layout,mpi_layout,mask,1); + GridRedBlackCartesian RBFine(&Fine,mask,1); GridParallelRNG FineRNG(&Fine); FineRNG.SeedFixedIntegers(std::vector({45,12,81,9})); diff --git a/tests/core/Test_cshift_red_black_rotate.cc b/tests/core/Test_cshift_red_black_rotate.cc index 3ef1cd21..aa9e6104 100644 --- a/tests/core/Test_cshift_red_black_rotate.cc +++ b/tests/core/Test_cshift_red_black_rotate.cc @@ -47,7 +47,7 @@ int main (int argc, char ** argv) mask[0]=0; GridCartesian Fine (latt_size,simd_layout,mpi_layout); - GridRedBlackCartesian RBFine(latt_size,simd_layout,mpi_layout,mask,1); + GridRedBlackCartesian RBFine(&Fine,mask,1); GridParallelRNG FineRNG(&Fine); FineRNG.SeedFixedIntegers(std::vector({45,12,81,9})); diff --git a/tests/core/Test_fft.cc b/tests/core/Test_fft.cc index 877683f0..b2336cfa 100644 --- a/tests/core/Test_fft.cc +++ b/tests/core/Test_fft.cc @@ -47,7 +47,7 @@ int main (int argc, char ** argv) vol = vol * latt_size[d]; } GridCartesian GRID(latt_size,simd_layout,mpi_layout); - GridRedBlackCartesian RBGRID(latt_size,simd_layout,mpi_layout); + GridRedBlackCartesian RBGRID(&GRID); LatticeComplexD one(&GRID); LatticeComplexD zz(&GRID); diff --git a/tests/core/Test_gpwilson_even_odd.cc b/tests/core/Test_gpwilson_even_odd.cc index fc12fe75..2069eb40 100644 --- a/tests/core/Test_gpwilson_even_odd.cc +++ b/tests/core/Test_gpwilson_even_odd.cc @@ -40,7 +40,7 @@ int main (int argc, char ** argv) std::vector simd_layout = GridDefaultSimd(Nd,vComplex::Nsimd()); std::vector mpi_layout = GridDefaultMpi(); GridCartesian Grid(latt_size,simd_layout,mpi_layout); - GridRedBlackCartesian RBGrid(latt_size,simd_layout,mpi_layout); + GridRedBlackCartesian RBGrid(&Grid); int threads = GridThread::GetThreads(); std::cout< simd_layout = GridDefaultSimd(Nd,vComplex::Nsimd()); std::vector mpi_layout = GridDefaultMpi(); GridCartesian Grid(latt_size,simd_layout,mpi_layout); - GridRedBlackCartesian RBGrid(latt_size,simd_layout,mpi_layout); + GridRedBlackCartesian RBGrid(&Grid); int threads = GridThread::GetThreads(); std::cout< simd_layout = GridDefaultSimd(Nd,vComplex::Nsimd()); std::vector mpi_layout = GridDefaultMpi(); GridCartesian Grid(latt_size,simd_layout,mpi_layout); - GridRedBlackCartesian RBGrid(latt_size,simd_layout,mpi_layout); + GridRedBlackCartesian RBGrid(&Grid); int threads = GridThread::GetThreads(); std::cout< simd_layout = GridDefaultSimd(Nd,vComplex::Nsimd()); std::vector mpi_layout = GridDefaultMpi(); GridCartesian Grid(latt_size,simd_layout,mpi_layout); - GridRedBlackCartesian RBGrid(latt_size,simd_layout,mpi_layout); + GridRedBlackCartesian RBGrid(&Grid); int threads = GridThread::GetThreads(); std::cout< mpi_layout = GridDefaultMpi(); GridCartesian Grid(latt_size,simd_layout,mpi_layout); - GridRedBlackCartesian RBGrid(latt_size,simd_layout,mpi_layout); + GridRedBlackCartesian RBGrid(&Grid); int threads = GridThread::GetThreads(); std::cout< mpi_layout = GridDefaultMpi(); GridCartesian Grid(latt_size,simd_layout,mpi_layout); - GridRedBlackCartesian RBGrid(latt_size,simd_layout,mpi_layout); + GridRedBlackCartesian RBGrid(&Grid); int threads = GridThread::GetThreads(); std::cout< mpi_layout = GridDefaultMpi(); GridCartesian Grid(latt_size,simd_layout,mpi_layout); - GridRedBlackCartesian RBGrid(latt_size,simd_layout,mpi_layout); + GridRedBlackCartesian RBGrid(&Grid); int threads = GridThread::GetThreads(); std::cout< mpi_layout = GridDefaultMpi(); GridCartesian Grid(latt_size,simd_layout,mpi_layout); - GridRedBlackCartesian RBGrid(latt_size,simd_layout,mpi_layout); + GridRedBlackCartesian RBGrid(&Grid); int threads = GridThread::GetThreads(); std::cout< simd_layout = GridDefaultSimd(Nd, vComplex::Nsimd()); std::vector mpi_layout = GridDefaultMpi(); GridCartesian Grid(latt_size, simd_layout, mpi_layout); - GridRedBlackCartesian RBGrid(latt_size, simd_layout, mpi_layout); + GridRedBlackCartesian RBGrid(&Grid); std::vector seeds({1, 2, 3, 4, 5}); GridSerialRNG sRNG; @@ -149,4 +149,4 @@ JSON } -*/ \ No newline at end of file +*/ diff --git a/tests/solver/Test_dwf_mrhs_cg.cc b/tests/solver/Test_dwf_mrhs_cg.cc new file mode 100644 index 00000000..2d2cfcb1 --- /dev/null +++ b/tests/solver/Test_dwf_mrhs_cg.cc @@ -0,0 +1,195 @@ + /************************************************************************************* + + Grid physics library, www.github.com/paboyle/Grid + + Source file: ./tests/Test_dwf_mrhs_cg.cc + + Copyright (C) 2015 + +Author: Peter Boyle + + 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 +#include + +using namespace std; +using namespace Grid; +using namespace Grid::QCD; + +int main (int argc, char ** argv) +{ + typedef typename DomainWallFermionR::FermionField FermionField; + typedef typename DomainWallFermionR::ComplexField ComplexField; + typename DomainWallFermionR::ImplParams params; + + const int Ls=8; + + Grid_init(&argc,&argv); + + std::vector latt_size = GridDefaultLatt(); + std::vector simd_layout = GridDefaultSimd(Nd,vComplex::Nsimd()); + std::vector mpi_layout = GridDefaultMpi(); + std::vector mpi_split (mpi_layout.size(),1); + + std::cout << "UGrid (world root)"<RankCount() ; + + ///////////////////////////////////////////// + // Split into 1^4 mpi communicators + ///////////////////////////////////////////// + std::cout << "SGrid (world root)"< seeds({1,2,3,4}); + + GridParallelRNG pRNG(UGrid ); pRNG.SeedFixedIntegers(seeds); + GridParallelRNG pRNG5(FGrid); pRNG5.SeedFixedIntegers(seeds); + std::vector src(nrhs,FGrid); + std::vector result(nrhs,FGrid); + + for(int s=0;sThisRank(); + LatticeGaugeField s_Umu(SGrid); + FermionField s_src(SFGrid); + FermionField s_res(SFGrid); + + { + FGrid->Barrier(); + ScidacWriter _ScidacWriter; + _ScidacWriter.open(file); + std::cout << GridLogMessage << "****************************************************************** "<Barrier(); + std::cout << GridLogMessage << "****************************************************************** "<Barrier(); + std::cout << GridLogMessage << "****************************************************************** "<Barrier(); + + std::cout << GridLogMessage << "****************************************************************** "<Barrier(); + } + + + /////////////////////////////////////////////////////////////// + // Set up N-solvers as trivially parallel + /////////////////////////////////////////////////////////////// + + RealD mass=0.01; + RealD M5=1.8; + DomainWallFermionR Ddwf(s_Umu,*SFGrid,*SFrbGrid,*SGrid,*SrbGrid,mass,M5); + + std::cout << GridLogMessage << "****************************************************************** "< HermOp(Ddwf); + ConjugateGradient CG((1.0e-8/(me+1)),10000); + s_res = zero; + CG(HermOp,s_src,s_res); + + /////////////////////////////////////// + // Share the information + /////////////////////////////////////// + std::vector iterations(nrhs,0); + iterations[me] = CG.IterationsToComplete; + + for(int n=0;nGlobalSum(iterations[n]); + } + + ///////////////////////////////////////////////////////////// + // Report how long they all took + ///////////////////////////////////////////////////////////// + for(int r=0;r simd_layout = GridDefaultSimd(Nd,vComplex::Nsimd()); std::vector mpi_layout = GridDefaultMpi(); GridCartesian Grid(latt_size,simd_layout,mpi_layout); - GridRedBlackCartesian RBGrid(latt_size,simd_layout,mpi_layout); + GridRedBlackCartesian RBGrid(&Grid); std::vector seeds({1,2,3,4,5}); GridParallelRNG pRNG(&Grid); pRNG.SeedFixedIntegers(seeds); diff --git a/tests/solver/Test_staggered_block_cg_unprec.cc b/tests/solver/Test_staggered_block_cg_unprec.cc index f54bc3b2..22051ef6 100644 --- a/tests/solver/Test_staggered_block_cg_unprec.cc +++ b/tests/solver/Test_staggered_block_cg_unprec.cc @@ -27,7 +27,6 @@ Author: Peter Boyle *************************************************************************************/ /* END LEGAL */ #include -#include using namespace std; using namespace Grid; diff --git a/tests/solver/Test_staggered_cg_prec.cc b/tests/solver/Test_staggered_cg_prec.cc index 0a803c21..d771f22e 100644 --- a/tests/solver/Test_staggered_cg_prec.cc +++ b/tests/solver/Test_staggered_cg_prec.cc @@ -57,7 +57,7 @@ int main (int argc, char ** argv) std::vector simd_layout = GridDefaultSimd(Nd,vComplex::Nsimd()); std::vector mpi_layout = GridDefaultMpi(); GridCartesian Grid(latt_size,simd_layout,mpi_layout); - GridRedBlackCartesian RBGrid(latt_size,simd_layout,mpi_layout); + GridRedBlackCartesian RBGrid(&Grid); std::vector seeds({1,2,3,4}); GridParallelRNG pRNG(&Grid); pRNG.SeedFixedIntegers(seeds); diff --git a/tests/solver/Test_staggered_cg_unprec.cc b/tests/solver/Test_staggered_cg_unprec.cc index 5e0358d7..eb33c004 100644 --- a/tests/solver/Test_staggered_cg_unprec.cc +++ b/tests/solver/Test_staggered_cg_unprec.cc @@ -57,7 +57,7 @@ int main (int argc, char ** argv) std::vector simd_layout = GridDefaultSimd(Nd,vComplex::Nsimd()); std::vector mpi_layout = GridDefaultMpi(); GridCartesian Grid(latt_size,simd_layout,mpi_layout); - GridRedBlackCartesian RBGrid(latt_size,simd_layout,mpi_layout); + GridRedBlackCartesian RBGrid(&Grid); std::vector seeds({1,2,3,4}); GridParallelRNG pRNG(&Grid); pRNG.SeedFixedIntegers(seeds); diff --git a/tests/solver/Test_wilson_cg_prec.cc b/tests/solver/Test_wilson_cg_prec.cc index 011bc70b..99ddfceb 100644 --- a/tests/solver/Test_wilson_cg_prec.cc +++ b/tests/solver/Test_wilson_cg_prec.cc @@ -52,7 +52,7 @@ int main (int argc, char ** argv) std::vector simd_layout = GridDefaultSimd(Nd,vComplex::Nsimd()); std::vector mpi_layout = GridDefaultMpi(); GridCartesian Grid(latt_size,simd_layout,mpi_layout); - GridRedBlackCartesian RBGrid(latt_size,simd_layout,mpi_layout); + GridRedBlackCartesian RBGrid(&Grid); std::vector seeds({1,2,3,4}); GridParallelRNG pRNG(&Grid); pRNG.SeedFixedIntegers(seeds); diff --git a/tests/solver/Test_wilson_cg_schur.cc b/tests/solver/Test_wilson_cg_schur.cc index 7bbf74d3..13ac0090 100644 --- a/tests/solver/Test_wilson_cg_schur.cc +++ b/tests/solver/Test_wilson_cg_schur.cc @@ -52,7 +52,7 @@ int main (int argc, char ** argv) std::vector simd_layout = GridDefaultSimd(Nd,vComplex::Nsimd()); std::vector mpi_layout = GridDefaultMpi(); GridCartesian Grid(latt_size,simd_layout,mpi_layout); - GridRedBlackCartesian RBGrid(latt_size,simd_layout,mpi_layout); + GridRedBlackCartesian RBGrid(&Grid); std::vector seeds({1,2,3,4}); GridParallelRNG pRNG(&Grid); pRNG.SeedFixedIntegers(seeds); diff --git a/tests/solver/Test_wilson_cg_unprec.cc b/tests/solver/Test_wilson_cg_unprec.cc index 19c5f854..db227ec8 100644 --- a/tests/solver/Test_wilson_cg_unprec.cc +++ b/tests/solver/Test_wilson_cg_unprec.cc @@ -52,7 +52,7 @@ int main (int argc, char ** argv) std::vector simd_layout = GridDefaultSimd(Nd,vComplex::Nsimd()); std::vector mpi_layout = GridDefaultMpi(); GridCartesian Grid(latt_size,simd_layout,mpi_layout); - GridRedBlackCartesian RBGrid(latt_size,simd_layout,mpi_layout); + GridRedBlackCartesian RBGrid(&Grid); std::vector seeds({1,2,3,4}); GridParallelRNG pRNG(&Grid); pRNG.SeedFixedIntegers(seeds); diff --git a/tests/solver/Test_wilson_cr_unprec.cc b/tests/solver/Test_wilson_cr_unprec.cc index 4182c04e..eccd7e74 100644 --- a/tests/solver/Test_wilson_cr_unprec.cc +++ b/tests/solver/Test_wilson_cr_unprec.cc @@ -52,7 +52,7 @@ int main (int argc, char ** argv) std::vector simd_layout = GridDefaultSimd(Nd,vComplex::Nsimd()); std::vector mpi_layout = GridDefaultMpi(); GridCartesian Grid(latt_size,simd_layout,mpi_layout); - GridRedBlackCartesian RBGrid(latt_size,simd_layout,mpi_layout); + GridRedBlackCartesian RBGrid(&Grid); std::vector seeds({1,2,3,4}); GridParallelRNG pRNG(&Grid); pRNG.SeedFixedIntegers(seeds);