diff --git a/Grid/build-lime-mpfr.sh b/Grid/build-lime-mpfr.sh new file mode 100644 index 000000000..427755bc9 --- /dev/null +++ b/Grid/build-lime-mpfr.sh @@ -0,0 +1,29 @@ +#! /bin/bash + +# Builds two packages needed by Grid + +################## +# prerequisites: Lime, MPFR place in $HOME/prefix +################## +export prefix=$HOME/prefix +export grid=$HOME/GridCompile +mkdir $prefix +################## +#LIME +################## +cd $prefix +wget http://usqcd-software.github.io/downloads/c-lime/lime-1.3.2.tar.gz +tar xvzf lime-1.3.2.tar.gz +cd lime-1.3.2 +./configure --prefix $prefix +make all install +################## +#MPFR - summit is badly configured +################## +cd $prefix +wget https://www.mpfr.org/mpfr-current/mpfr-4.0.2.tar.gz +tar xvzf mpfr-4.0.2.tar.gz +cd mpfr-4.0.2 +./configure --prefix $prefix +make all install + diff --git a/Grid/build.sh b/Grid/build.sh index 31ab0dc11..c0f244cf1 100755 --- a/Grid/build.sh +++ b/Grid/build.sh @@ -3,7 +3,9 @@ ARCH=$1 PK_CC=$2 PK_CXX=$3 -GIT_BRANCH=feature/staggered-comms-compute +#GIT_REPO=https://github.com/paboyle/Grid +GIT_REPO=https://github.com/milc-qcd/Grid +GIT_BRANCH=CGinfo if [ -z ${PK_CXX} ] then @@ -12,11 +14,11 @@ then fi case ${ARCH} in - scalar|avx512|avx2) + scalar|avx512|avx2|gpu-cuda) ;; *) echo "Unsupported ARCH" - echo "Usage $0 " + echo "Usage $0 " exit 1 esac @@ -27,12 +29,12 @@ SRCDIR=${TOPDIR}/Grid BUILDDIR=${TOPDIR}/build-${ARCH} INSTALLDIR=${TOPDIR}/install-${ARCH} -MAKE="make -j4" +MAKE="make -j4 V=1" if [ ! -d ${SRCDIR} ] then echo "Fetching ${GIT_BRANCH} branch of Grid package from github" - git clone https://github.com/paboyle/Grid -b ${GIT_BRANCH} + git clone ${GIT_REPO} -b ${GIT_BRANCH} fi # Fetch Eigen package, set up Make.inc files and create Grid configure @@ -51,17 +53,16 @@ then scalar) - module swap craype-mic-knl craype-haswell - ${SRCDIR}/configure \ --prefix=${INSTALLDIR} \ --enable-precision=double \ --enable-simd=GEN \ --enable-comms=none \ - --with-lime=${HOME}/scidac/install/qio-cori-omp-knl-icc \ + --with-lime=${HOME}/scidac/install/qio-single \ + --with-fftw=${HOME}/fftw/build-gcc \ --with-openssl=/global/common/cori/software/openssl/1.1.0a/hsw \ CXX="${PK_CXX}" \ - CXXFLAGS="-std=c++11" \ + CXXFLAGS="-std=c++14" \ # --with-hdf5=/opt/cray/pe/hdf5/1.10.0/INTEL/15.0 \ @@ -72,20 +73,20 @@ then ${SRCDIR}/configure \ --prefix=${INSTALLDIR} \ + --enable-mkl=yes \ --enable-precision=double \ --enable-simd=GEN \ --enable-comms=mpi \ - --with-lime=${HOME}/scidac/install/qio-cori-omp-knl-icc \ + --with-lime=${HOME}/scidac/install/qio-skx \ --with-openssl=/global/common/cori/software/openssl/1.1.0a/hsw \ CXX="${PK_CXX}" CC="${PK_CC}" \ CXXFLAGS="-std=c++11 -xCORE-AVX2" \ # --with-hdf5=/opt/cray/pe/hdf5/1.10.0/INTEL/15.0 \ -# --enable-mkl=yes \ status=$? ;; - avx512) + avx512-knl) INCMKL="-I/opt/intel/compilers_and_libraries_2018.1.163/linux/mkl/include" LIBMKL="-L/opt/intel/compilers_and_libraries_2018.1.163/linux/mkl/lib/intel64_lin" @@ -96,16 +97,50 @@ then --enable-simd=KNL \ --enable-comms=mpi \ --host=x86_64-unknown-linux-gnu \ - --with-lime=${HOME}/scidac/install/qio-cori-omp-knl-icc \ - --with-openssl=/global/common/cori/software/openssl/1.1.0a/hsw \ + --with-lime=${HOME}/scidac/install/qio-impi-knl \ CXX="${PK_CXX}" CC="${PK_CC}" \ - CXXFLAGS="-std=c++11 -xMIC-AVX512" \ + CXXFLAGS="-std=c++17 -xMIC-AVX512 -O2 -g -vec -simd -qopenmp" \ # --with-hdf5=/opt/cray/pe/hdf5/1.10.0.3/INTEL/16.0 \ + # --with-openssl=/global/common/cori/software/openssl/1.1.0a/hsw \ status=$? echo "Configure exit status $status" - ;; + ;; + avx512-skx) + + INCMKL="-I/opt/intel/compilers_and_libraries_2018.1.163/linux/mkl/include" + LIBMKL="-L/opt/intel/compilers_and_libraries_2018.1.163/linux/mkl/lib/intel64_lin" + + ${SRCDIR}/configure \ + --prefix=${INSTALLDIR} \ + --enable-precision=double \ + --enable-simd=KNL \ + --enable-comms=mpi \ + --host=x86_64-unknown-linux-gnu \ + --with-lime=${HOME}/scidac/install/qio-impi-knl \ + CXX="${PK_CXX}" CC="${PK_CC}" \ + CXXFLAGS="-std=c++17 -xCORE-AVX512 -O2 -g -vec -simd -qopenmp" \ + + # --with-hdf5=/opt/cray/pe/hdf5/1.10.0.3/INTEL/16.0 \ + # --with-openssl=/global/common/cori/software/openssl/1.1.0a/hsw \ + + status=$? + echo "Configure exit status $status" + ;; + gpu-cuda) + ${SRCDIR}/configure \ + --prefix ${INSTALLDIR} \ + --enable-precision=double \ + --enable-simd=GEN \ + --enable-comms=mpi \ + --host=x86_64-unknown-linux-gnu \ + CXX=nvcc \ + LDFLAGS=-L$HOME/prefix/lib/ \ + CXXFLAGS="-ccbin ${PK_CXX} -gencode arch=compute_70,code=sm_70 -I$HOME/prefix/include/ -std=c++11" + status=$? + echo "Configure exit status $status" + ;; *) echo "Unsupported ARCH ${ARCH}" exit 1; diff --git a/file_utilities/Make_template b/file_utilities/Make_template index 7c4cd5b37..e906a6a1b 100644 --- a/file_utilities/Make_template +++ b/file_utilities/Make_template @@ -133,7 +133,7 @@ check_gauge:: "DEFINES= -DONLY_GLUON_FILES" \ "EXTRA_OBJECTS= ${G_OBJECTS} check_gauge.o nersc_cksum.o \ gauge_utilities.o \ - check_unitarity.o gauge_info_dummy.o d_plaq4.o io_lat4.o io_scidac.o" + check_unitarity.o gauge_info_dummy.o d_plaq4.o io_lat4.o" diff_colormatrix:: ${MAKE} -f ${MAKEFILE} target "MYTARGET= diff_colormatrix" \ diff --git a/generic/gridMap.cc b/generic/gridMap.cc index dfbc763a1..732f4e8af 100644 --- a/generic/gridMap.cc +++ b/generic/gridMap.cc @@ -17,20 +17,22 @@ extern int even_sites_on_node; /* number of even sites on this node */ extern int this_node; using namespace Grid; -using namespace Grid::QCD; +//using namespace Grid::QCD; using namespace std; -extern std::vector squaresize; +extern Coordinate squaresize; static void -indexToCoords(uint64_t idx, std::vector &x){ +indexToCoords(uint64_t idx, Coordinate &x){ + + int r[4]; // Gets the lattice coordinates from the MILC index - get_coords(x.data(), this_node, idx); + get_coords(r, this_node, idx); // For Grid, we need the coordinates within the sublattice hypercube for the current MPI rank // NOTE: Would be better to provide a get_subl_coords() in MILC layout_*.c for(int i = 0; i < 4; i++) - x[i] %= squaresize[i]; + x[i] = r[i] % squaresize[i]; //printf("Converted %d to %d %d %d %d\n", idx, x[0], x[1], x[2], x[3]); fflush(stdout); } @@ -41,7 +43,6 @@ indexToCoords(uint64_t idx, std::vector &x){ template static struct GRID_ColorVector_struct * create_V(int milc_parity, GridCartesian *CGrid, GridRedBlackCartesian *RBGrid){ - struct GRID_ColorVector_struct *out; out = (struct GRID_ColorVector_struct *) @@ -54,13 +55,13 @@ create_V(int milc_parity, GridCartesian *CGrid, GridRedBlackCartesian *RBGrid){ case EVEN: out->cv = new typename ImprovedStaggeredFermion::FermionField(RBGrid); GRID_ASSERT(out->cv != NULL, GRID_MEM_ERROR); - out->cv->checkerboard = Even; + out->cv->Checkerboard() = Even; break; case ODD: out->cv = new typename ImprovedStaggeredFermion::FermionField(RBGrid); GRID_ASSERT(out->cv != NULL, GRID_MEM_ERROR); - out->cv->checkerboard = Odd; + out->cv->Checkerboard() = Odd; break; case EVENANDODD: @@ -92,7 +93,7 @@ create_nV(int n, int milc_parity, std::cout << "Constructing 5D field\n" << std::flush; out->cv = new typename ImprovedStaggeredFermion5D::FermionField(FRBGrid); GRID_ASSERT(out->cv != NULL, GRID_MEM_ERROR); - out->cv->checkerboard = milc_parity == EVEN ? Even : Odd ; + out->cv->Checkerboard() = milc_parity == EVEN ? Even : Odd ; } else { out->cv = new typename ImprovedStaggeredFermion5D::FermionField(FCGrid); GRID_ASSERT(out->cv != NULL, GRID_MEM_ERROR); @@ -142,10 +143,10 @@ create_V_from_vec( su3_vector *src, int milc_parity, int loopstart=((milc_parity)==ODD ? even_sites_on_node : 0 ); auto start = std::chrono::system_clock::now(); - PARALLEL_FOR_LOOP + #pragma omp parallel for for( uint64_t idx = loopstart; idx < loopend; idx++){ - std::vector x(4); + Coordinate x(4); indexToCoords(idx,x); ColourVector cVec; @@ -184,17 +185,20 @@ create_nV_from_vecs( su3_vector *src[], int n, int milc_parity, std::cout << "create_nv_from_vecs: ColourVector size = " << sizeof(ColourVector) - << " ColourVectorField size =" << sizeof(*(out->cv)) << "\n" << std::flush; + << " ColourVectorField size = " << sizeof(*(out->cv)) << "\n" << std::flush; auto start = std::chrono::system_clock::now(); - PARALLEL_FOR_LOOP + #pragma omp parallel for for( uint64_t idx = loopstart; idx < loopend; idx++){ - std::vector x(4); + Coordinate x(4); indexToCoords(idx,x); - std::vector x5(1,0); +// Coordinate x5(1,0); +// for( int d = 0; d < 4; d++ ) +// x5.push_back(x[d]); + Coordinate x5(5); for( int d = 0; d < 4; d++ ) - x5.push_back(x[d]); - + x5[d+1] = x[d]; + for( int j = 0; j < n; j++ ){ x5[0] = j; ColourVector cVec; @@ -222,7 +226,7 @@ static void extract_V_to_vec( su3_vector *dest, FORSOMEFIELDPARITY_OMP(idx, milc_parity, ) { - std::vector x(4); + Coordinate x(4); indexToCoords(idx, x); ColourVector cVec; @@ -247,9 +251,9 @@ static void extract_nV_to_vecs( su3_vector *dest[], int n, FORSOMEFIELDPARITY_OMP(idx, milc_parity, ) { - std::vector x(4); + Coordinate x(4); indexToCoords(idx, x); - std::vector x5(1,0); + Coordinate x5(1,0); for( int d = 0; d < 4; d++ ) x5.push_back(x[d]); @@ -289,13 +293,13 @@ static void milcGaugeFieldToGrid(su3_matrix *in, LatticeGaugeField &out){ typedef typename LatticeGaugeField::vector_object vobj; typedef typename vobj::scalar_object sobj; - GridBase *grid = out._grid; + GridBase *grid = out.Grid(); int lsites = grid->lSites(); std::vector scalardata(lsites); - PARALLEL_FOR_LOOP + #pragma omp parallel for for (uint64_t milc_idx = 0; milc_idx < sites_on_node; milc_idx++){ - std::vector x(4); + Coordinate x(4); indexToCoords(milc_idx, x); int grid_idx; Lexicographic::IndexFromCoor(x, grid_idx, grid->LocalDimensions()); @@ -372,10 +376,10 @@ asqtad_destroy_L( struct GRID_FermionLinksAsqtad_struct *Link // Create a 4D, full-grid wrapper GRID_4Dgrid * GRID_create_grid(void){ - std::vector latt_size = GridDefaultLatt(); - std::vector simd_layoutF = GridDefaultSimd(Nd,vComplexF::Nsimd()); - std::vector simd_layoutD = GridDefaultSimd(Nd,vComplexD::Nsimd()); - std::vector mpi_layout = GridDefaultMpi(); + Coordinate latt_size = GridDefaultLatt(); + Coordinate simd_layoutF = GridDefaultSimd(Nd,vComplexF::Nsimd()); + Coordinate simd_layoutD = GridDefaultSimd(Nd,vComplexD::Nsimd()); + Coordinate mpi_layout = GridDefaultMpi(); GridCartesian *CGridF = new GridCartesian(latt_size,simd_layoutF,mpi_layout); GRID_ASSERT(CGridF != NULL, GRID_MEM_ERROR); diff --git a/generic/milc_to_grid_utilities.cc b/generic/milc_to_grid_utilities.cc index ea828bdb0..6e6029443 100644 --- a/generic/milc_to_grid_utilities.cc +++ b/generic/milc_to_grid_utilities.cc @@ -23,12 +23,12 @@ extern "C" { } using namespace std; using namespace Grid; -using namespace Grid::QCD; +//using namespace Grid::QCD; GRID_4Dgrid *grid_full; GRID_4DRBgrid *grid_rb; -std::vector squaresize; +Coordinate squaresize; static int grid_is_initialized = 0; @@ -120,8 +120,8 @@ initialize_grid(void){ grid_full = GRID_create_grid(); grid_rb = GRID_create_RBgrid(grid_full); - std::vector latt_size = GridDefaultLatt(); - std::vector mpi_layout = GridDefaultMpi(); + Coordinate latt_size = GridDefaultLatt(); + Coordinate mpi_layout = GridDefaultMpi(); if(mynode()==0){ printf("milc_to_grid_utilities: Initialized Grid with args\n%s %s\n%s %s\n%s %s\n", @@ -149,7 +149,7 @@ void setup_grid_communicator(int peGrid[]){ terminate(1); } - std::vector processors; + Coordinate processors; for(int i=0;i<4;i++) processors.push_back(peGrid[i]); grid_cart = new Grid::CartesianCommunicator(processors); } @@ -205,7 +205,8 @@ int grid_rank_from_processor_coor(int x, int y, int z, int t){ terminate(1); } - std::vector coor = {x, y, z, t}; + Coordinate coor(4); + coor[0] = x; coor[1] = y; coor[2] = z; coor[3] = t; int worldrank = grid_cart->RankFromProcessorCoor(coor); return worldrank; } @@ -221,7 +222,7 @@ void grid_coor_from_processor_rank(int coords[], int worldrank){ terminate(1); } - std::vector coor; + Coordinate coor; grid_cart->ProcessorCoorFromRank(worldrank, coor); for(int i = 0; i < 4; i++) diff --git a/generic_ks/d_congrad5_fn_gpu.c b/generic_ks/d_congrad5_fn_gpu.c index 58d54eaa5..dd7c4ddc8 100644 --- a/generic_ks/d_congrad5_fn_gpu.c +++ b/generic_ks/d_congrad5_fn_gpu.c @@ -142,7 +142,6 @@ int ks_congrad_parity_gpu(su3_vector *t_src, su3_vector *t_dest, &residual, &relative_residual, &num_iters); - qic->final_rsq = residual*residual; qic->final_relrsq = relative_residual*relative_residual; @@ -249,13 +248,12 @@ int ks_congrad_block_parity_gpu(int nsrc, su3_vector **t_src, su3_vector **t_des const int quda_precision = qic->prec; double residual, relative_residual; - int num_iters; + int num_iters = 0; // for newer versions of QUDA we need to invalidate the gauge field if the links are new - static imp_ferm_links_t *fn_last = NULL; - if ( fn != fn_last || fresh_fn_links(fn) ){ + if ( fn != get_fn_last() || fresh_fn_links(fn) ){ cancel_quda_notification(fn); - fn_last = fn; + set_fn_last(fn); num_iters = -1; node0_printf("%s: fn, notify: Signal QUDA to refresh links\n", myname); } @@ -299,8 +297,8 @@ int ks_congrad_block_parity_gpu(int nsrc, su3_vector **t_src, su3_vector **t_des #ifdef CGTIME if(this_node==0){ - printf("CONGRAD5: time = %e (fn_QUDA %s) masses = 1 iters = %d mflops = %e\n", - dtimec, prec_label[quda_precision-1], qic->final_iters, + printf("CONGRAD5: time = %e (fn_QUDA %s) masses = 1 srcs = %d iters = %d mflops = %e\n", + dtimec, prec_label[quda_precision-1], nsrc,qic->final_iters, (double)(nflop*volume*qic->final_iters/(1.0e6*dtimec*numnodes())) ); fflush(stdout);} #endif diff --git a/generic_ks/gridStaggInvert.cc b/generic_ks/gridStaggInvert.cc index 168e5c01a..ba3610324 100644 --- a/generic_ks/gridStaggInvert.cc +++ b/generic_ks/gridStaggInvert.cc @@ -13,7 +13,7 @@ #include "../include/macros.h" using namespace Grid; -using namespace Grid::QCD; +//using namespace Grid::QCD; template static void @@ -52,7 +52,7 @@ asqtadInvert (GRID_info_t *info, struct GRID_FermionLinksAsqtad_structparity == GRID_EVENODD) || (inv_arg->parity == GRID_EVEN) || (inv_arg->parity == GRID_ODD), GRID_FAIL); // In and out fields must be on the same lattice - GRID_ASSERT(in->cv->_grid == out->cv->_grid, GRID_FAIL); + GRID_ASSERT(in->cv->Grid() == out->cv->Grid(), GRID_FAIL); auto start = std::chrono::system_clock::now(); // Call using c1 = c2 = 2. and u0 = 1. to neutralize link rescaling -- probably ignored anyway. @@ -70,7 +70,7 @@ asqtadInvert (GRID_info_t *info, struct GRID_FermionLinksAsqtad_structcv->_grid == CGrid) && (out->cv->_grid == CGrid), GRID_FAIL); + GRID_ASSERT((in->cv->Grid() == CGrid) && (out->cv->Grid() == CGrid), GRID_FAIL); std::cout << "WARNING: inversion with EVENODD is untested.\n"; MdagMLinearOperator HermOp(Ds); @@ -90,8 +90,8 @@ asqtadInvert (GRID_info_t *info, struct GRID_FermionLinksAsqtad_structcv->_grid == RBGrid) && (out->cv->_grid == RBGrid) - && (in->cv->checkerboard == out->cv->checkerboard), GRID_FAIL); + GRID_ASSERT((in->cv->Grid() == RBGrid) && (out->cv->Grid() == RBGrid) + && (in->cv->Checkerboard() == out->cv->Checkerboard()), GRID_FAIL); SchurStaggeredOperator HermOp(Ds); @@ -125,7 +125,7 @@ asqtadInvertMulti (GRID_info_t *info, struct GRID_FermionLinksAsqtad_structcv->_grid == out[i]->cv->_grid, GRID_FAIL); + GRID_ASSERT(in->cv->Grid() == out[i]->cv->Grid(), GRID_FAIL); } typedef typename ImprovedStaggeredFermion::FermionField FermionField; @@ -161,7 +161,7 @@ asqtadInvertMulti (GRID_info_t *info, struct GRID_FermionLinksAsqtad_structcv->_grid == CGrid), GRID_FAIL); + GRID_ASSERT((in->cv->Grid() == CGrid), GRID_FAIL); std::vector outvec(nmass, CGrid); MdagMLinearOperator HermOp(Ds); @@ -171,9 +171,9 @@ asqtadInvertMulti (GRID_info_t *info, struct GRID_FermionLinksAsqtad_structfinal_iter = MSCG.IterationsToComplete[i]; - iters += MSCG.IterationsToComplete[i]; - res_arg[i]->final_rsq = MSCG.TrueResidual[i]*MSCG.TrueResidual[i]; + res_arg[i]->final_iter = MSCG.IterationsToCompleteShift[i]; + iters += MSCG.IterationsToCompleteShift[i]; + res_arg[i]->final_rsq = MSCG.TrueResidualShift[i]*MSCG.TrueResidualShift[i]; } auto elapsed = end - start; info->final_sec = std::chrono::duration_cast(elapsed).count()/1000.; @@ -189,10 +189,10 @@ asqtadInvertMulti (GRID_info_t *info, struct GRID_FermionLinksAsqtad_structcv->_grid == RBGrid, GRID_FAIL); + GRID_ASSERT(in->cv->Grid() == RBGrid, GRID_FAIL); std::vector outvec(nmass, RBGrid); for(int i = 0; i < nmass; i++) - outvec[i].checkerboard = inv_arg->parity == EVEN ? Even : Odd ; + outvec[i].Checkerboard() = inv_arg->parity == EVEN ? Even : Odd ; SchurStaggeredOperator HermOp(Ds); // Do the solve @@ -206,9 +206,9 @@ asqtadInvertMulti (GRID_info_t *info, struct GRID_FermionLinksAsqtad_structfinal_iter = MSCG.IterationsToComplete[i]; - iters += MSCG.IterationsToComplete[i]; - res_arg[i]->final_rsq = MSCG.TrueResidual[i]*MSCG.TrueResidual[i]; + res_arg[i]->final_iter = MSCG.IterationsToCompleteShift[i]; + iters += MSCG.IterationsToCompleteShift[i]; + res_arg[i]->final_rsq = MSCG.TrueResidualShift[i]*MSCG.TrueResidualShift[i]; } std::cout << "iters summed over masses = " << iters << "\n" << std::flush; @@ -238,7 +238,7 @@ asqtadInvertBlock (GRID_info_t *info, GridCartesian *CGrid, GridRedBlackCartesian *RBGrid) { // In and out fields must be on the same lattice - GRID_ASSERT(in->cv->_grid == out->cv->_grid, GRID_FAIL); + GRID_ASSERT(in->cv->Grid() == out->cv->Grid(), GRID_FAIL); typedef typename ImprovedStaggeredFermion5D::FermionField FermionField; typedef typename ImprovedStaggeredFermion5D::ComplexField ComplexField; diff --git a/generic_ks/io_helpers_ks_eigen.c b/generic_ks/io_helpers_ks_eigen.c index d834bd417..c353062ec 100644 --- a/generic_ks/io_helpers_ks_eigen.c +++ b/generic_ks/io_helpers_ks_eigen.c @@ -209,7 +209,7 @@ int reload_ks_eigen(int flag, char *eigfile, int *Nvecs, double *eigVal, } } else { /* Read using QUDA format */ - qio_status = read_quda_ks_eigenvectors(infile, eigVec, eigVal, Nvecs, EVEN, fn); + qio_status = read_quda_ks_eigenvectors(infile, eigVec, eigVal, Nvecs, EVEN); if(qio_status != QIO_SUCCESS){ if(qio_status == QIO_EOF){ node0_printf("WARNING: Premature EOF at eigenvectors\n"); @@ -222,6 +222,10 @@ int reload_ks_eigen(int flag, char *eigfile, int *Nvecs, double *eigVal, } } close_ks_eigen_infile(infile); + + /* Reconstruct eigenvalues */ + reset_eigenvalues( eigVec, eigVal, *Nvecs, EVEN, fn); + break; default: node0_printf("%s: Unrecognized reload flag.\n", myname); @@ -247,7 +251,7 @@ int reload_ks_eigen(int flag, char *eigfile, int *Nvecs, double *eigVal, >1 for seek, read error, or missing data error */ int reload_ks_eigen(int flag, char *eigfile, int *Nvecs, double *eigVal, - su3_vector **eigVec, int timing){ + su3_vector **eigVec, imp_ferm_links_t *fn, int timing){ register int i, j; int status = 0; diff --git a/generic_ks/io_scidac_ks_eigen.c b/generic_ks/io_scidac_ks_eigen.c index 36eb0162b..1d8a8c5d0 100644 --- a/generic_ks/io_scidac_ks_eigen.c +++ b/generic_ks/io_scidac_ks_eigen.c @@ -360,7 +360,7 @@ read_ks_eigenvector(QIO_Reader *infile, int packed, su3_vector *eigVec, double * int read_quda_ks_eigenvectors(QIO_Reader *infile, su3_vector *eigVec[], double *eigVal, int *Nvecs, - int parity, imp_ferm_links_t *fn){ + int parity){ char myname[] = "read_quda_ks_eigenvectors"; int status; char *xml; @@ -422,9 +422,6 @@ read_quda_ks_eigenvectors(QIO_Reader *infile, su3_vector *eigVec[], double *eigV } } END_LOOP_OMP; - /* Reconstruct eigenvalues */ - reset_eigenvalues( eigVec, eigVal, *Nvecs, parity, fn); - free(eigVecs); return status; diff --git a/generic_ks/mat_invert.c b/generic_ks/mat_invert.c index 4b94ba20f..309743780 100644 --- a/generic_ks/mat_invert.c +++ b/generic_ks/mat_invert.c @@ -534,13 +534,11 @@ void check_invert_field( su3_vector *src, su3_vector *dest, Real mass, tmp = create_v_field(); /* Compute tmp = M src */ - node0_printf("check_invert_field: calling ks_dirac_op\n"); fflush(stdout); ks_dirac_op( src, tmp, mass, parity, fn); sum2=sum=0.0; dmaxerr=0; flag = 0; - node0_printf("check_invert_field: checking diffs\n"); fflush(stdout); FORSOMEFIELDPARITY(i,parity){ for(k=0;k<3;k++){ r_diff = dest[i].c[k].real - tmp[i].c[k].real; diff --git a/include/generic_quda.h b/include/generic_quda.h index e7cc42abb..4fad8b369 100644 --- a/include/generic_quda.h +++ b/include/generic_quda.h @@ -116,16 +116,4 @@ static void destroy_M_quda(anti_hermitmat *momentum) { qudaFreePinned(momentum); } -/* - Return the most recent fermion link field passed to QUDA - (defined in generic_ks/ks_multicg_offset_gpu.c) -*/ -imp_ferm_links_t* get_fn_last(); - -/* - Update the fermion link field passed to QUDA - (defined in generic_ks/ks_multicg_offset_gpu.c) -*/ -void set_fn_last(imp_ferm_links_t *fn_last_new); - #endif /* GENERIC_QUDA_H */ diff --git a/include/imp_ferm_links.h b/include/imp_ferm_links.h index 7b609ee31..cd23a7cec 100644 --- a/include/imp_ferm_links.h +++ b/include/imp_ferm_links.h @@ -281,6 +281,18 @@ int mat_invert_multi( imp_ferm_links_t *fn[] /* Storage for fat and Naik links */ ); +/* + Return the most recent fermion link field passed to QUDA + (defined in generic_ks/ks_multicg_offset_gpu.c) +*/ +imp_ferm_links_t* get_fn_last(); + +/* + Update the fermion link field passed to QUDA + (defined in generic_ks/ks_multicg_offset_gpu.c) +*/ +void set_fn_last(imp_ferm_links_t *fn_last_new); + /* eigen_stuff*.c */ typedef struct { diff --git a/include/io_scidac_ks.h b/include/io_scidac_ks.h index 84afa5e54..2e1f72cf1 100644 --- a/include/io_scidac_ks.h +++ b/include/io_scidac_ks.h @@ -56,7 +56,7 @@ void close_ks_eigen_outfile(QIO_Writer *outfile); QIO_Reader *open_ks_eigen_infile(char *filename, int *Nvecs, int *packed, int *file_type, int serpar); int read_ks_eigenvector(QIO_Reader *infile, int packed, su3_vector *eigVec, double *eigVal); int read_quda_ks_eigenvectors(QIO_Reader *infile, su3_vector *eigVec[], double *eigVal, int *Nvecs, - int parity, imp_ferm_links_t *fn); + int parity); void close_ks_eigen_infile(QIO_Reader *infile); /**********************************************************************/ diff --git a/include/mGrid/mGrid_internal.h b/include/mGrid/mGrid_internal.h index 4dbbff9f7..6544468e0 100644 --- a/include/mGrid/mGrid_internal.h +++ b/include/mGrid/mGrid_internal.h @@ -4,7 +4,7 @@ #include using namespace Grid; -using namespace Grid::QCD; +//using namespace Grid::QCD; // Containers for grid definitions diff --git a/include/milc_datatypes.h b/include/milc_datatypes.h index 906ef0897..1d9a87c71 100644 --- a/include/milc_datatypes.h +++ b/include/milc_datatypes.h @@ -33,7 +33,7 @@ typedef struct { #define complex dcomplex #endif -#ifdef WANT_QUDA +#ifdef HAVE_QUDA // When using QUDA, we need to set the site struct member arrays alignment to a multiple of 16 bytes #define ALIGNAS(n) __attribute__((aligned(n))) #define ALIGNMENT ALIGNAS(16) diff --git a/ks_spectrum/Make_template b/ks_spectrum/Make_template index 16b56792c..2b4b120f4 100644 --- a/ks_spectrum/Make_template +++ b/ks_spectrum/Make_template @@ -72,9 +72,12 @@ DEFLATE_OBJECTS = \ eigen_stuff_helpers.o \ io_helpers_ks_eigen.o \ io_ks_eigen.o \ - io_scidac_ks_eigen.o \ jacobi.o +ifeq ($(strip ${HAVEQIO}),true) + DEFLATE_OBJECTS += io_scidac_ks_eigen.o +endif + OBJECTS = \ ${MY_OBJECTS} \ ${GAUGE_OBJECTS} \ @@ -165,14 +168,14 @@ quark_action.h: ${QUARKIMP}/${QUARK} ks_spectrum_asqtad:: ${MAKE} -f ${MAKEFILE} target "MYTARGET= $@" \ "QUARK = asqtad_action.h" \ - "DEFINES= -DFN -DHAVE_KS -DEIGEN_QIO" \ + "DEFINES= -DFN -DHAVE_KS" \ "EXTRA_OBJECTS += ${FN_OBJECTS}" ks_spectrum_eigcg_asqtad:: ${MAKE} -f ${MAKEFILE} target "MYTARGET= $@" \ "QUARK = asqtad_action.h" \ "DEFINES= -DFN -DHAVE_KS" \ - "DEFINES += -DEIGMODE=EIGCG -DEIGEN_QIO" \ + "DEFINES += -DEIGMODE=EIGCG" \ "EXTRA_OBJECTS += ${FN_OBJECTS} ${EIGCG_OBJECTS}" ks_spectrum_asqtad_eos:: @@ -221,12 +224,12 @@ HISQ_OPTIONS = "DEFINES= -DFN -DHAVE_KS \ ks_spectrum_hisq:: ${MAKE} -f ${MAKEFILE} target "MYTARGET= $@" \ - ${HISQ_OPTIONS} "DEFINES += -DEIGEN_QIO" + ${HISQ_OPTIONS} ks_spectrum_eigcg_hisq:: ${MAKE} -f ${MAKEFILE} target "MYTARGET= $@" \ - ${HISQ_OPTIONS} "DEFINES += -DEIGMODE=EIGCG -DEIGEN_QIO" \ + ${HISQ_OPTIONS} "DEFINES += -DEIGMODE=EIGCG" \ "EXTRA_OBJECTS += ${EIGCG_OBJECTS}" ks_spectrum_hisq_eos:: diff --git a/ks_spectrum/control.c b/ks_spectrum/control.c index 96b95d4ff..1328e892b 100644 --- a/ks_spectrum/control.c +++ b/ks_spectrum/control.c @@ -342,9 +342,9 @@ int main(int argc, char *argv[]) for(k=0; k 1: + requires.append(self.load[1]) + pass + if len(self.save) > 1: + produces.append(self.save[1]) + pass + return { 'name': wname, 'depends': depends, + 'requires': requires, 'produces': produces } + pass + +class PointSource: + """Point base source.""" + _Template = """ + #== source ${id}: ${_classType} == + point + field_type ${field_type} + subset ${subset} + origin #echo ' '.join(map(str,$origin))# + #if $scaleFactor is not None: + scale_factor ${scaleFactor} + #end if + source_label ${label} + #echo ' '.join($save)#""" + def __init__(self,origin,field_type,subset,scaleFactor,label,save): + self._classType = self.__class__.__name__ + self._objectID = self._classType+'_'+base36(id(self)) + self.origin = origin + self.field_type = field_type + self.subset = subset + self.scaleFactor = scaleFactor + self.label = label + self.save = save + self._template = Template(source=textwrap.dedent(self._Template),searchList=vars(self)) + pass + def generate(self,ostream): + print(self._template, file=ostream) + return + def dataflow(self): + wname = self._objectID + depends = list() + requires = list() + produces = list() + return { 'name': wname, 'depends': depends, + 'requires': requires, 'produces': produces } + pass + +class RandomColorWallSource: + """Random color wall base source.""" + _Template = """ + #== source ${id}: ${_classType} == + random_color_wall + field_type ${field_type} + subset ${subset} + t0 ${tsrc} + ncolor ${ncolor} + momentum #echo ' '.join(map(str,$momentum))# + #if $scaleFactor is not None: + scale_factor ${scaleFactor} + #end if + source_label ${label} + #echo ' '.join($save)#""" + def __init__(self,tsrc,ncolor,field_type,subset,momentum,scaleFactor,label,save): + self._classType = self.__class__.__name__ + self._objectID = self._classType+'_'+base36(id(self)) + self.tsrc = tsrc + self.ncolor = ncolor + self.field_type = field_type + self.subset = subset + self.momentum = momentum + self.scaleFactor = scaleFactor + self.label = label + self.save = save + self._template = Template(source=textwrap.dedent(self._Template),searchList=vars(self)) + return + def generate(self,ostream): + print(self._template, file=ostream) + return + def dataflow(self): + wname = self._objectID + depends = list() + requires = list() + produces = list() + if len(self.save) > 1: + produces.append(self.save[1]) + pass + return { 'name': wname, 'depends': depends, + 'requires': requires, 'produces': produces } + pass + +class VectorFieldSource: + """Color vector field base source input from a file.""" + _Template = """ + #== source ${id}: ${_classType} == + vector_field + field_type ${field_type} + subset ${subset} + origin #echo ' '.join(map(str,$origin))# + #echo ' '.join($load)# + ncolor ${ncolor} + momentum #echo ' '.join(map(str,$momentum))# + #if $scaleFactor is not None: + scale_factor ${scaleFactor} + #end if + source_label ${label} + #echo ' '.join($save)#""" + def __init__(self,load,origin,ncolor,field_type,subset,momentum,scaleFactor,label,save): + self._classType = self.__class__.__name__ + self._objectID = self._classType+'_'+base36(id(self)) + self.load = load + self.origin = origin + self.ncolor = ncolor + self.field_type = field_type + self.subset = subset + self.momentum = momentum + self.scaleFactor = scaleFactor + self.label = label + self.save = save + self._template = Template(source=textwrap.dedent(self._Template),searchList=vars(self)) + return + def generate(self,ostream): + print(self._template, file=ostream) + return + def dataflow(self): + wname = self._objectID + depends = list() + requires = list() + produces = list() + if len(self.load) > 1: + requires.append(self.load[1]) + pass + if len(self.save) > 1: + produces.append(self.save[1]) + pass + return { 'name': wname, 'depends': depends, + 'requires': requires, 'produces': produces } + pass + +class DiracFieldSource: + """Dirac (spin and color) base source.""" + _Template = """ + #== source ${id}: ${_classType} == + dirac_field + field_type clover + subset ${subset} + origin #echo ' '.join(map(str,$origin))# + #echo ' '.join($load)# + nsource #echo 4 * ${ncolor}# + momentum #echo ' '.join(map(str,$momentum))# + #if $scaleFactor is not None: + scale_factor ${scaleFactor} + #end if + source_label ${label} + #echo ' '.join($save)#""" + def __init__(self,load,origin,ncolor,subset,momentum,scaleFactor,label,save): + self._classType = self.__class__.__name__ + self._objectID = self._classType+'_'+base36(id(self)) + self.load = load + self.origin = origin + self.ncolor = ncolor + self.subset = subset + self.momentum = momentum + self.scaleFactor = scaleFactor + self.label = label + self.save = save + self._template = Template(source=textwrap.dedent(self._Template),searchList=vars(self)) + return + def generate(self,ostream): + print(self._template, file=ostream) + return + def dataflow(self): + wname = self._objectID + depends = list() + requires = list() + produces = list() + if len(self.load) > 1: + requires.append(self.load[1]) + pass + if len(self.save) > 1: + produces.append(self.save[1]) + pass + return { 'name': wname, 'depends': depends, + 'requires': requires, 'produces': produces } + pass + +class RadialWavefunction: + """Radial wavefunction source from a file. Either as a base or a derived source type.""" + _Template = """ + #== source ${id}: ${_classType} == + #if $startSource is not None: + source ${startSource.id} + #end if + wavefunction + #echo ' '.join($load)# + #if $stride is not None: + stride ${stride} + #end if + a ${afm} + op_label ${label} + #echo ' '.join($save)#""" + def __init__(self,label,afm,stride,load,save,startSource=None): + self._classType = self.__class__.__name__ + self._objectID = self._classType+'_'+base36(id(self)) + self.label = label + self.afm = afm + self.stride = stride + self.load = load + self.save = save + self.startSource = startSource + self._template = Template(source=textwrap.dedent(self._Template),searchList=vars(self)) + return + def generate(self,ostream): + print(self._template, file=ostream) + return + def dataflow(self): + wname = self._objectID + depends = list() + requires = list() + produces = list() + if len(self.load) > 1: + requires.append(self.load[1]) + pass + if len(self.save) > 1: + produces.append(self.save[1]) + pass + if self.startSource is not None: + depends.append(self.startSource._objectID) + pass + return { 'name': wname, 'depends': depends, + 'requires': requires, 'produces': produces } + pass + +class FermilabRotation: + """Apply a Fermilab rotation to a source.""" + _Template = """ + #== source ${id}: ${_classType} == + #if $startSource is not None: + source ${startSource.id} + #end if + rotate_3D + d1 ${d1} + op_label ${label} + #echo ' '.join($save)#""" + def __init__(self,label,d1,save,startSource=None): + self._classType = self.__class__.__name__ + self._objectID = self._classType+'_'+base36(id(self)) + self.label = label + self.d1 = d1 + self.save = save + self.startSource = startSource + self._template = Template(source=textwrap.dedent(self._Template),searchList=vars(self)) + return + def generate(self,ostream): + print(self._template, file=ostream) + return + def dataflow(self): + wname = self._objectID + depends = list() + requires = list() + produces = list() + if len(self.load) > 1: + requires.append(self.load[1]) + pass + if len(self.save) > 1: + produces.append(self.save[1]) + pass + if self.startSource is not None: + depends.append(self.startSource._objectID) + pass + return { 'name': wname, 'depends': depends, + 'requires': requires, 'produces': produces } + pass + +class FatCovariantGaussian: + """Covariant Gaussian source from APE smeared links. Either as a base or a derived source type.""" + _Template = """ + #== source ${id}: ${_classType} == + #if $startSource is not None: + source ${startSource.id} + #end if + fat_covariant_gaussian + stride ${gparams.stride} + r0 ${gparams.r0} + source_iters ${gparams.iters} + op_label ${label} + #echo ' '.join($save)#""" + def __init__(self,gparams,label,save,startSource=None): + self._classType = self.__class__.__name__ + self._objectID = self._classType+'_'+base36(id(self)) + self.gparams = gparams + self.label = label + self.save = save + self.startSource = startSource + self._template = Template(source=textwrap.dedent(self._Template),searchList=vars(self)) + return + def generate(self,ostream): + print(self._template, file=ostream) + return + def dataflow(self): + wname = self._objectID + depends = list() + requires = list() + produces = list() + if len(self.save) > 1: + produces.append(self.save[1]) + pass + if self.startSource is not None: + depends.append(self.startSource._objectID) + pass + return { 'name': wname, 'depends': depends, + 'requires': requires, 'produces': produces } + pass + +class Eigen: + """Read eigenpairs from file.""" + _Template = """ + #== ${_classType} == + max_number_of_eigenpairs ${Nvecs} + #if $Nvecs > 0: + #echo ' '.join($load)# + #echo ' '.join($save)# + #end if""" + def __init__(self,load, Nvecs, save): + self._classType = self.__class__.__name__ + self._objectID = self._classType+'_'+base36(id(self)) + self.Nvecs = Nvecs + self.load = load + self.save =save + self._template = Template(source=textwrap.dedent(self._Template),searchList=vars(self)) + return + def generate(self,ostream): + print(self._template, file=ostream) + return + def dataflow(self): + wname = self._objectID + depends = list() + requires = list() + produces = list() + if len(self.load) > 1: + requires.append(self.load[1]) + pass + if len(self.save) > 1: + produces.append(self.save[1]) + pass + depends.append(self.source._objectID) + return { 'name': wname, 'depends': depends, + 'requires': requires, 'produces': produces } + pass + + +class SolveKS: + """su3_clov KS propagator solve.""" + _Template = """ + #== propagator ${id}: ${_classType} == + propagator_type KS + mass ${mass} + #if $naik_epsilon is not None: + naik_term_epsilon ${naik_epsilon} + #end if + check ${check} + error_for_propagator ${residual.L2} + rel_error_for_propagator ${residual.R2} + precision ${precision} + momentum_twist #echo ' '.join(map(str,$twist))# + source ${source.id} + #echo ' '.join($load)# + #echo ' '.join($save)#""" + def __init__(self,mass,naik_epsilon,source,twist,load,save,residual,precision,check): + self._classType = self.__class__.__name__ + self._objectID = self._classType+'_'+base36(id(self)) + self.mass = mass + self.naik_epsilon = naik_epsilon + self.source = source + self.twist = twist + self.load = load + self.save = save + self.residual = residual + self.precision = precision + self.check = check + self._template = Template(source=textwrap.dedent(self._Template),searchList=vars(self)) + return + def generate(self,ostream): + print(self._template, file=ostream) + return + def dataflow(self): + wname = self._objectID + depends = list() + requires = list() + produces = list() + if len(self.load) > 1: + requires.append(self.load[1]) + pass + if len(self.save) > 1: + produces.append(self.save[1]) + pass + depends.append(self.source._objectID) + return { 'name': wname, 'depends': depends, + 'requires': requires, 'produces': produces } + pass + +class SolveClover: + """su3_clov Clover propagator solve.""" + _Template = """ + #== propagator ${id}: ${_classType} == + propagator_type clover + kappa ${kappa} + clov_c ${cSW} + check ${check} + error_for_propagator ${residual.L2} + rel_error_for_propagator ${residual.R2} + precision ${precision} + momentum_twist #echo ' '.join(map(str,$twist))# + source ${source.id} + #echo ' '.join($load)# + #echo ' '.join($save)#""" + def __init__(self,kappa,cSW,source,twist,load,save,residual,precision,check): + self._classType = self.__class__.__name__ + self._objectID = self._classType+'_'+base36(id(self)) + self.kappa = kappa + self.cSW = cSW + self.source = source + self.twist = twist + self.load = load + self.save =save + self.residual = residual + self.precision = precision + self.check = check + self._template = Template(source=textwrap.dedent(self._Template),searchList=vars(self)) + return + def generate(self,ostream): + print(self._template, file=ostream) + return + def dataflow(self): + wname = self._objectID + depends = list() + requires = list() + produces = list() + if len(self.load) > 1: + requires.append(self.load[1]) + pass + if len(self.save) > 1: + produces.append(self.save[1]) + pass + depends.append(self.source._objectID) + return { 'name': wname, 'depends': depends, + 'requires': requires, 'produces': produces } + pass + +class QuarkIdentitySink: + """NOP fermion sink smearing.""" + _Template = """ + #== quark ${id}: ${_classType} == + ${prop.type} ${prop.id} + identity + op_label ${label} + #echo ' '.join($save)#""" + def __init__(self,prop,label,save): + self._classType = self.__class__.__name__ + self._objectID = self._classType+'_'+base36(id(self)) + self.prop = prop + self.label = label + self.save = save + self._template = Template(source=textwrap.dedent(self._Template),searchList=vars(self)) + return + def generate(self,ostream): + print(self._template, file=ostream) + return + def dataflow(self): + wname = self._objectID + depends = list() + requires = list() + produces = list() + depends.append(self.prop._objectID) + if len(self.save) > 1: + produces.append(self.save[1]) + pass + return { 'name': wname, 'depends': depends, + 'requires': requires, 'produces': produces } + pass + +class RadialWavefunctionSink: + """Smear fermion sink with a radial wavefunction input from file.""" + _Template = """ + #== quark ${id}: ${_classType} == + ${prop.type} ${prop.id} + wavefunction + #echo ' '.join($load)# + #if $stride is not None: + stride ${stride} + #end if + a ${afm} + op_label ${label} + #echo ' '.join($save)#""" + def __init__(self,prop,label,afm,stride,load,save): + self._classType = self.__class__.__name__ + self._objectID = self._classType+'_'+base36(id(self)) + self.prop = prop + self.label = label + self.afm = afm + self.stride = stride + self.load = load + self.save = save + self._template = Template(source=textwrap.dedent(self._Template),searchList=vars(self)) + return + def generate(self,ostream): + print(self._template, file=ostream) + return + def dataflow(self): + wname = self._objectID + depends = list() + requires = list() + produces = list() + depends.append(self.prop._objectID) + if len(self.load) > 1: + requires.append(self.load[1]) + pass + if len(self.save) > 1: + produces.append(self.save[1]) + pass + return { 'name': wname, 'depends': depends, + 'requires': requires, 'produces': produces } + pass + +class FermilabRotateSink: + """Apply a Fermilab rotation to a propagator.""" + _Template = """ + #== quark ${id}: ${_classType} == + ${prop.type} ${prop.id} + rotate_3D + d1 ${d1} + op_label ${label} + #echo ' '.join($save)#""" + def __init__(self,prop,label,d1,save): + self._classType = self.__class__.__name__ + self._objectID = self._classType+'_'+base36(id(self)) + self.prop = prop + self.label = label + self.d1 = d1 + self.save = save + self._template = Template(source=textwrap.dedent(self._Template),searchList=vars(self)) + return + def generate(self,ostream): + print(self._template, file=ostream) + return + def dataflow(self): + wname = self._objectID + depends = list() + requires = list() + produces = list() + depends.append(self.prop._objectID) + if len(self.load) > 1: + requires.append(self.load[1]) + pass + if len(self.save) > 1: + produces.append(self.save[1]) + pass + return { 'name': wname, 'depends': depends, + 'requires': requires, 'produces': produces } + pass + +class DiracExtSrcSink: + """Restrict Dirac fermion propagator to time slice.""" + _Template = """ + #== quark ${id}: ${_classType} == + ${prop.type} ${prop.id} + ext_src_dirac + gamma ${gamma} + momentum #echo ' '.join(map(str,$momentum))# + t0 ${time_slice} + op_label ${label} + #echo ' '.join($save)#""" + def __init__(self,prop,gamma,momentum,time_slice,label,save): + self._classType = self.__class__.__name__ + self._objectID = self._classType+'_'+base36(id(self)) + self.prop = prop + self.label = label + self.gamma = gamma + self.momentum = momentum + self.time_slice = time_slice + self.save = save + self._template = Template(source=textwrap.dedent(self._Template),searchList=vars(self)) + return + def generate(self,ostream): + print(self._template, file=ostream) + return + def dataflow(self): + wname = self._objectID + depends = list() + requires = list() + produces = list() + depends.append(self.prop._objectID) + if len(self.save) > 1: + produces.append(self.save[1]) + pass + return { 'name': wname, 'depends': depends, + 'requires': requires, 'produces': produces } + pass + +class KSExtSrcSink: + """Restrict KS fermion propagator to time slice.""" + _Template = """ + #== quark ${id}: ${_classType} == + ${prop.type} ${prop.id} + ext_src_ks + spin_taste_extend ${spin_taste_op} + momentum #echo ' '.join(map(str,$momentum))# + t0 ${time_slice} + op_label ${label} + #echo ' '.join($save)#""" + def __init__(self,prop,spin_taste_op,momentum,time_slice,label,save): + self._classType = self.__class__.__name__ + self._objectID = self._classType+'_'+base36(id(self)) + self.prop = prop + self.label = label + self.spin_taste_op = spin_taste_op + self.momentum = momentum + self.time_slice = time_slice + self.save = save + self._template = Template(source=textwrap.dedent(self._Template),searchList=vars(self)) + return + def generate(self,ostream): + print(self._template, file=ostream) + return + def dataflow(self): + wname = self._objectID + depends = list() + requires = list() + produces = list() + depends.append(self.prop._objectID) + if len(self.save) > 1: + produces.append(self.save[1]) + pass + return { 'name': wname, 'depends': depends, + 'requires': requires, 'produces': produces } + pass + +class KSInverseSink: + """Extended KS inverse sink operator.""" + _Template = """ + #== quark ${id}: ${_classType} == + ${prop.type} ${prop.id} + ks_inverse + mass ${mass} + #if $naik_epsilon is not None: + naik_term_epsilon ${naik_epsilon} + #end if + u0 ${u0} + max_cg_iterations ${maxCG.iters} + max_cg_restarts ${maxCG.restarts} + deflate ${deflate} + error_for_propagator ${residual.L2} + rel_error_for_propagator ${residual.R2} + precision ${precision} + momentum_twist #echo ' '.join(map(str,$twist))# + op_label ${label} + #echo ' '.join($save)#""" + def __init__(self,prop,mass,naik_epsilon,u0,maxCG,deflate,residual,precision,twist,label,save): + self._classType = self.__class__.__name__ + self._objectID = self._classType+'_'+base36(id(self)) + self.prop = prop + self.label = label + self.mass = mass + self.naik_epsilon = naik_epsilon + self.u0 = u0 + self.maxCG = maxCG + self.deflate = deflate + self.residual = residual + self.precision = precision + self.twist = twist + self.save = save + self._template = Template(source=textwrap.dedent(self._Template),searchList=vars(self)) + return + def generate(self,ostream): + print(self._template, file=ostream) + return + def dataflow(self): + wname = self._objectID + depends = list() + requires = list() + produces = list() + depends.append(self.prop._objectID) + if len(self.save) > 1: + produces.append(self.save[1]) + pass + return { 'name': wname, 'depends': depends, + 'requires': requires, 'produces': produces } + pass + +class DiracInverseSink: + """Extended Dirac inverse sink operator.""" + _Template = """ + #== quark ${id}: ${_classType} == + ${prop.type} ${prop.id} + dirac_inverse + kappa ${kappa} + clov_c ${clov_c} + u0 ${u0} + max_cg_iterations ${maxCG.iters} + max_cg_restarts ${maxCG.restarts} + error_for_propagator ${residual.L2} + rel_error_for_propagator ${residual.R2} + precision ${precision} + momentum_twist #echo ' '.join(map(str,$twist))# + op_label ${label} + #echo ' '.join($save)#""" + def __init__(self,prop,kappa,clov_c,u0,maxCG,residual,precision,twist,label,save): + self._classType = self.__class__.__name__ + self._objectID = self._classType+'_'+base36(id(self)) + self.prop = prop + self.label = label + self.kappa = kappa + self.clov_c = clov_c + self.u0 = u0 + self.maxCG = maxCG + self.residual = residual + self.precision = precision + self.twist = twist + self.save = save + self._template = Template(source=textwrap.dedent(self._Template),searchList=vars(self)) + return + def generate(self,ostream): + print(self._template, file=ostream) + return + def dataflow(self): + wname = self._objectID + depends = list() + requires = list() + produces = list() + depends.append(self.prop._objectID) + if len(self.save) > 1: + produces.append(self.save[1]) + pass + return { 'name': wname, 'depends': depends, + 'requires': requires, 'produces': produces } + pass + +class MesonSpectrum: + """Meson spectrum specification.""" + _Template = """ + #== ${_classType} == + pair ${antiQuark.id} ${quark.id} + spectrum_request meson + #echo ' '.join($save)# + r_offset #echo ' '.join(map(str,$relOffset))# + number_of_correlators #echo len($npts)""" + def __init__(self,antiQuark,quark,relOffset,npts,save): + self._classType = self.__class__.__name__ + self._objectID = self._classType+'_'+base36(id(self)) + self.antiQuark = antiQuark + self.quark = quark + self.relOffset = relOffset + self.npts = npts + self.save = save + self._template = Template(source=textwrap.dedent(self._Template),searchList=vars(self)) + return + def generate(self,ostream): + print(self._template, file=ostream) + for c in self.npts: + c.generate(ostream) + return + def dataflow(self): + wname = self._objectID + depends = list() + requires = list() + produces = list() + depends.append(self.antiQuark._objectID) + depends.append(self.quark._objectID) + if len(self.save) > 1: + produces.append(self.save[1]) + pass + return { 'name': wname, 'depends': depends, + 'requires': requires, 'produces': produces } + pass + +class MesonNpt: + """A meson n-point function.""" + _Template = """correlator ${prefix} ${postfix} #echo ' '.join(map(str,$norm))# #echo ' '.join($gamma)# #echo ' '.join(map(str,$momentum))# #echo ' '.join($parity)""" + def __init__(self,prefix,postfix,norm,gamma,momentum,parity): + self._classType = self.__class__.__name__ + self._objectID = self._classType+'_'+base36(id(self)) + self.prefix = prefix + self.postfix = postfix + self.norm = norm + self.gamma = gamma + self.momentum = momentum + self.parity = parity + self._template = Template(source=textwrap.dedent(self._Template),searchList=vars(self)) + return + def generate(self,ostream): + print(self._template, file=ostream) + return + +class _Cycle_su3_clov: + def __init__(self): + self.gauge = None + self.maxIters = None + self.maxRestarts = None + self.bsource = list() + self.msource = list() + self.propagator = list() + self.quark = list() + self.spectrum = list() + return + def generate(self,ostream): + self.gauge.generate(ostream) + print('max_cg_iterations', self.maxIters, file=ostream) + print('max_cg_restarts', self.maxRestarts, file=ostream) + print('########### base sources ###############', file=ostream) + print("",file=ostream) + print('number_of_base_sources', len(self.bsource), file=ostream) + for x in self.bsource: + x.generate(ostream) + pass + print('########### modified sources ###############', file=ostream) + print("",file=ostream) + print('number_of_modified_sources', len(self.msource), file=ostream) + for x in self.msource: + x.generate(ostream) + pass + print("",file=ostream) + print('########### propagators ###############', file=ostream) + print("",file=ostream) + print('number_of_propagators', len(self.propagator), file=ostream) + for x in self.propagator: + x.generate(ostream) + print("",file=ostream) + print('########### quarks ###############', file=ostream) + print("",file=ostream) + print('number_of_quarks', len(self.quark), file=ostream) + for x in self.quark: + x.generate(ostream) + pass + print("",file=ostream) + print('########### mesons ###############', file=ostream) + print("",file=ostream) + print('number_of_pairings', len(self.spectrum), file=ostream) + for x in self.spectrum: + x.generate(ostream) + return + def dataflow(self): + dinfo = list() + dinfo.append(self.gauge.dataflow()) + for x in self.bsource: + dinfo.append(x.dataflow()) + pass + for x in self.msource: + dinfo.append(x.dataflow()) + pass + for x in self.propagator: + dinfo.append(x.dataflow()) + pass + for x in self.quark: + dinfo.append(x.dataflow()) + pass + for x in self.spectrum: + dinfo.append(x.dataflow()) + pass + return dinfo + pass + +class su3_clov: + """ + The MILC/su3_clov application. + + :Parameters: + - `participantName`: A unique label for the workflow. + - `dim`: List of lattice dimensions, [nX, nY, nZ, nT]. + - `seed`: Random number generator seed, use seed='None' to skip. + - `jobID`: A unique identifier string used for data provenance, usually set to PBS_JOBID. + - `layout`: The SciDAC node and io layout = { node: [Lx, Ly, Lz, Lt], io: [Ix, Iy, Iz, It] }. Use layout = 'None' to skip. + - `prompt`: Echo MILC prompts to stdout (=0), no echo (=1) or check input validity (=2). + + The 'su3_clov' application is a pipline: + + - initialize lattice size, rng.seed + - loop: + - initialize gauge field + - gauge fix + - fat link APE smear + - define base sources + - define derived sources + - define propagator solves + - define quarks (propagator sink treatments) + - spectroscopy: n-point tie-ups + + """ + def __init__(self,participantName,dim,seed,jobID,layout=None,prompt=0): + self._classType = self.__class__.__name__ + self._objectID = self._classType+'_'+base36(id(self)) + self.pName = participantName + self.geometry = Geometry(dim,seed,jobID,prompt,layout) + self.requires = list() + self.produces = list() + self.cycle = [ ] + self.cycnt = -1 + return + def newGauge(self,uspec): + """Add the gauge definition.""" + self.cycle.append(_Cycle_su3_clov()) + self.cycnt += 1 + self.cycle[self.cycnt].gauge = uspec + return uspec + def setCGparams(self,maxIters,maxRestarts): + """Specify CG maximum restarts and iterations between restarts.""" + self.cycle[self.cycnt].maxIters = maxIters + self.cycle[self.cycnt].maxRestarts = maxRestarts + return + def addBaseSource(self,src): + """Add a base source specification""" + self.cycle[self.cycnt].bsource.append(src) + return src + def addModSource(self,src): + """Add a source derived from one of the base sources""" + self.cycle[self.cycnt].msource.append(src) + return src + def addProp(self,prop): + """Add a propagator solve.""" + prop.type = 'propagator' + self.cycle[self.cycnt].propagator.append(prop) + return prop + def addQuark(self,quark): + """Add a sink treatment to a propagator.""" + quark.type = 'quark' + self.cycle[self.cycnt].quark.append(quark) + return quark + def addSpectrum(self,spect): + """Add a spectrum specification that depends upon pairs of quarks.""" + self.cycle[self.cycnt].spectrum.append(spect) + return spect + def _bind_indices(self): + """Bind indices to sources, propagators and quarks.""" + for cy in self.cycle: + nbs = len(cy.bsource) + for idx, s in zip(range(nbs),cy.bsource): + s.id = idx + pass + nms = len(cy.msource) + for idx, s in zip(range(nbs,nbs+nms),cy.msource): + s.id = idx + pass + for idx, p in zip(range(len(cy.propagator)),cy.propagator): + p.id = idx + pass + for idx, q in zip(range(len(cy.quark)),cy.quark): + q.id = idx + pass + return + def generate(self,ostream=sys.stdout): + """Write MILC prompts to output stream 'ostream'.""" + self._bind_indices() + self.geometry.generate(ostream) + for x in self.cycle: + x.generate(ostream) + pass + return + def dataflow(self): + self._bind_indices() + w = list() + for x in self.cycle: + w.append(x.dataflow()) + pass + dinfo = { 'title': self.pName, 'workflow': w } + return dinfo + pass + +#---- + +class KSsolveSet: + """A set of KS solves that have a common source specification, momentum twist, and precision.""" + _Template = """ + #== ${_classType} == + max_cg_iterations ${maxCG.iters} + max_cg_restarts ${maxCG.restarts} + check ${check} + momentum_twist #echo ' '.join(map(str,$twist))# + precision ${precision} + source ${source.id}""" + def __init__(self,source,twist,check,maxCG,precision): + self._classType = self.__class__.__name__ + self._objectID = self._classType+'_'+base36(id(self)) + self.source = source + self.twist = twist + self.check = check + self.maxCG = maxCG + self.precision = precision + self.propagator = list() + self._template = Template(source=textwrap.dedent(self._Template),searchList=vars(self)) + return + # container behavior + def __len__(self): + return len(self.propagator) + def __iter__(self): + return self.propagator.__iter__() + def __add__(self,other): + return [ x for x in self.propagator+other.propagator ] + def addPropagator(self,prop): + """Add a KSsolveElement object to the set of solves.""" + prop.parent = self + prop.type = 'propagator' + self.propagator.append(prop) + return prop + def generate(self,ostream): + print(self._template, file=ostream) + print('number_of_propagators', len(self.propagator), file=ostream) + for p in self.propagator: + p.generate(ostream) + pass + return + def dataflow(self): + df = list() + wname = self._objectID + depends = list() + requires = list() + produces = list() + depends.append(self.source._objectID) + df.append( { 'name': wname, 'depends': depends, + 'requires': requires, 'produces': produces } ) + for p in self.propagator: + df.append(p.dataflow()) + pass + return df + pass + +class KSsolveElement: + """Specification of a single KS solve in a KSsolveSet.""" + _Template = """ + #== propagator ${id}: ${_classType} == + mass ${mass} + #if $naik is not None + naik_term_epsilon ${naik} + #end if + #if $deflate is not None: + deflate ${deflate} + #end if + error_for_propagator ${residual.L2} + rel_error_for_propagator ${residual.R2} + #echo ' '.join($load)# + #echo ' '.join($save)# + """ + def __init__(self,mass,naik,load,save,deflate,residual): + self._classType = self.__class__.__name__ + self._objectID = self._classType+'_'+base36(id(self)) + self.id = None + self.parent = None # the KSsolveSet + self.mass = mass + self.naik = naik + self.load = load + self.save = save + self.deflate = deflate + self.residual = residual + self._template = Template(source=textwrap.dedent(self._Template),searchList=vars(self)) + return + def generate(self,ostream): + print(self._template, file=ostream) + return + def dataflow(self): + wname = self._objectID + depends = list() + requires = list() + produces = list() + depends.append(self.parent._objectID) + if len(self.load) > 1: + requires.append(self.load[1]) + pass + if len(self.save) > 1: + produces.append(self.save[1]) + pass + return { 'name': wname, 'depends': depends, + 'requires': requires, 'produces': produces } + pass + +class ks_spectrum: + """ + MILC application MILC/ks_spectrum application. + The milc applications are designed as pipelines. + The ks_spectrum pipline: + + - initialize lattice size, rng + - loop: + - init gauge field + - gauge fix + - link smear + - eigenpairs?? + - pbp masses?? + - define sources + - define KS propsets + - define quarks (sink treatments) + - meson spectroscopy + - baryon spectroscopy + + """ + def __init__(self,participantName,dim,seed,jobID,layout=None,prompt=0): + self._classType = self.__class__.__name__ + self._objectID = self._classType+'_'+base36(id(self)) + self.pName = participantName + self.geometry = Geometry(dim,seed,jobID,prompt,layout) + self.requires = list() + self.produces = list() + self.cycle = [ ] + self.cycnt = -1 + return + def newGauge(self,uspec): + """Add a gauge field specification.""" + self.cycle.append(_Cycle_ks_spectrum()) + self.cycnt += 1 + self.cycle[self.cycnt].gauge = uspec + return uspec + def newEigen(self,eigspec): + """Add an eigensolve specification.""" + self.cycle[self.cycnt].eigen = eigspec + return eigspec + def addBaseSource(self,src): + """Add a base source specification.""" + self.cycle[self.cycnt].bsource.append(src) + return src + def addModSource(self,src): + """Add a modified source dependent upon one of the base sources.""" + self.cycle[self.cycnt].msource.append(src) + return src + def addPropSet(self,pset): + """Add a KSsolveSet object.""" + self.cycle[self.cycnt].pset.append(pset) + return pset + def addQuark(self,quark): + """Add a propagator sink treatment (a quark) to the workflow. The quark depends upon an KSsolveElement object.""" + quark.type = 'quark' + self.cycle[self.cycnt].quark.append(quark) + return quark + def addMeson(self,meson): + """Add meson spectroscopy.""" + self.cycle[self.cycnt].meson.append(meson) + return meson + def addBaryon(self,baryon): + """Add baryon spectroscopy (*NOTE* unimplemented).""" + raise + def _bind_indices(self): + """Bind indices to sources, propagators and quarks.""" + for cy in self.cycle: + nbs = len(cy.bsource) + for idx, s in zip(range(nbs),cy.bsource): + s.id = idx + pass + nms = len(cy.msource) + for idx, s in zip(range(nbs,nbs+nms),cy.msource): + s.id = idx + pass + idx = 0 + for ps in cy.pset: + for p in ps: + p.id = idx + idx += 1 + pass + pass + for idx, q in zip(range(len(cy.quark)),cy.quark): + q.id = idx + pass + return + def generate(self,ostream=sys.stdout): + """Write MILC prompts to output stream 'ostream'.""" + self._bind_indices() + self.geometry.generate(ostream) + for x in self.cycle: + x.generate(ostream) + pass + return + def dataflow(self): + self._bind_indices() + w = list() + for x in self.cycle: + w.append(x.dataflow()) + pass + dinfo = { 'title': self.pName, 'workflow': w } + return dinfo + pass + +class _Cycle_ks_spectrum: + def __init__(self): + self.gauge = None + self.eigen = None + self.bsource = list() + self.msource = list() + self.nextpropidx = 0 + self.pset = list() + self.quark = list() + self.meson = list() + self.baryon = list() + return + def generate(self,ostream): + self.gauge.generate(ostream) + self.eigen.generate(ostream) + print('#== PBP Masses ==', file=ostream) + print("",file=ostream) + print('number_of_pbp_masses 0', file=ostream) + print("",file=ostream) + print('#== Base Sources ==', file=ostream) + print("",file=ostream) + print('number_of_base_sources', len(self.bsource), file=ostream) + for x in self.bsource: + x.generate(ostream) + pass + print("",file=ostream) + print('#== Modified Sources ==', file=ostream) + print("",file=ostream) + print('number_of_modified_sources', len(self.msource), file=ostream) + for x in self.msource: + x.generate(ostream) + pass + print("",file=ostream) + print('#== KSsolveSets ==', file=ostream) + print("",file=ostream) + print('number_of_sets', len(self.pset), file=ostream) + for x in self.pset: + x.generate(ostream) + print("",file=ostream) + print('#== Quarks ==', file=ostream) + print("",file=ostream) + print('number_of_quarks', len(self.quark), file=ostream) + for x in self.quark: + x.generate(ostream) + pass + print("",file=ostream) + print('#== Mesons ==', file=ostream) + print("",file=ostream) + print('number_of_mesons', len(self.meson), file=ostream) + for x in self.meson: + x.generate(ostream) + print("",file=ostream) + print('#== Baryons ==', file=ostream) + print("",file=ostream) + print('number_of_baryons', len(self.baryon), file=ostream) + for x in self.baryon: + x.generate(ostream) + return + def dataflow(self): + dinfo = list() + dinfo.append(self.gauge.dataflow()) + for x in self.bsource: + dinfo.append(x.dataflow()) + pass + for x in self.msource: + dinfo.append(x.dataflow()) + pass + for x in self.pset: + for y in x.dataflow(): + dinfo.append(y) + pass + pass + for x in self.quark: + dinfo.append(x.dataflow()) + pass + for x in self.meson: + dinfo.append(x.dataflow()) + pass + for x in self.baryon: + dinfo.append(x.dataflow()) + pass + return dinfo + pass