diff --git a/README.md b/README.md index faf86faf..c7461368 100644 --- a/README.md +++ b/README.md @@ -20,7 +20,7 @@ License: GPL v2. Last update Nov 2016. -_Please send all pull requests to the `develop` branch._ +_Please do not send pull requests to the `master` branch which is reserved for releases._ ### Bug report @@ -29,7 +29,7 @@ _To help us tracking and solving more efficiently issues with Grid, please repor When you file an issue, please go though the following checklist: 1. Check that the code is pointing to the `HEAD` of `develop` or any commit in `master` which is tagged with a version number. -2. Give a description of the target platform (CPU, network, compiler). +2. Give a description of the target platform (CPU, network, compiler). Please give the full CPU part description, using for example `cat /proc/cpuinfo | grep 'model name' | uniq` (Linux) or `sysctl machdep.cpu.brand_string` (macOS) and the full output the `--version` option of your compiler. 3. Give the exact `configure` command used. 4. Attach `config.log`. 5. Attach `config.summary`. @@ -45,7 +45,7 @@ are provided, similar to HPF and cmfortran, and user control is given over the m array indices to both MPI tasks and SIMD processing elements. * Identically shaped arrays then be processed with perfect data parallelisation. -* Such identically shapped arrays are called conformable arrays. +* Such identically shaped arrays are called conformable arrays. The transformation is based on the observation that Cartesian array processing involves identical processing to be performed on different regions of the Cartesian array. @@ -127,14 +127,15 @@ make -C tests/ tests The following options can be use with the `--enable-simd=` option to target different communication interfaces: -| `` | Description | -| ------------- | -------------------------------------------- | -| `none` | no communications | -| `mpi[-auto]` | MPI communications | -| `mpi3[-auto]` | MPI communications using MPI 3 shared memory | -| `shmem ` | Cray SHMEM communications | +| `` | Description | +| -------------- | ------------------------------------------------------------- | +| `none` | no communications | +| `mpi[-auto]` | MPI communications | +| `mpi3[-auto]` | MPI communications using MPI 3 shared memory | +| `mpi3l[-auto]` | MPI communications using MPI 3 shared memory and leader model | +| `shmem ` | Cray SHMEM communications | -For `mpi` and `mpi3` the optional `-auto` suffix instructs the `configure` scripts to determine all the necessary compilation and linking flags. This is done by extracting the informations from the MPI wrapper specified in the environment variable `MPICXX` (if not specified `configure` will scan though a list of default names). +For the MPI interfaces the optional `-auto` suffix instructs the `configure` scripts to determine all the necessary compilation and linking flags. This is done by extracting the informations from the MPI wrapper specified in the environment variable `MPICXX` (if not specified `configure` will scan though a list of default names). ### Possible SIMD types @@ -160,7 +161,7 @@ Alternatively, some CPU codenames can be directly used: | `BGQ` | Blue Gene/Q | #### Notes: -- We currently support AVX512 only for the Intel compiler. Support for GCC and clang will appear in future versions. +- We currently support AVX512 only for the Intel compiler. Support for GCC and clang will appear in future versions of Grid when the AVX512 support within GCC and clang will be more advanced. - For BG/Q only [bgclang](http://trac.alcf.anl.gov/projects/llvm-bgq) is supported. We do not presently plan to support more compilers for this platform. - BG/Q performances are currently rather poor. This is being investigated for future versions. @@ -171,7 +172,7 @@ The following configuration is recommended for the Intel Knights Landing platfor ``` bash ../configure --enable-precision=double\ --enable-simd=KNL \ - --enable-comms=mpi3-auto \ + --enable-comms=mpi-auto \ --with-gmp= \ --with-mpfr= \ --enable-mkl \ @@ -183,10 +184,9 @@ where `` is the UNIX prefix where GMP and MPFR are installed. If you are w ``` bash ../configure --enable-precision=double\ --enable-simd=KNL \ - --enable-comms=mpi3 \ + --enable-comms=mpi \ --with-gmp= \ --with-mpfr= \ --enable-mkl \ CXX=CC CC=cc -``` - +``` \ No newline at end of file diff --git a/configure.ac b/configure.ac index a8ad40dd..2b759dde 100644 --- a/configure.ac +++ b/configure.ac @@ -4,6 +4,7 @@ AC_CANONICAL_BUILD AC_CANONICAL_HOST AC_CANONICAL_TARGET AM_INIT_AUTOMAKE(subdir-objects) +AM_EXTRA_RECURSIVE_TARGETS([tests]) AC_CONFIG_MACRO_DIR([m4]) AC_CONFIG_SRCDIR([lib/Grid.h]) AC_CONFIG_HEADERS([lib/Config.h]) @@ -253,18 +254,23 @@ AC_ARG_ENABLE([comms],[AC_HELP_STRING([--enable-comms=none|mpi|mpi-auto|mpi3|mpi case ${ac_COMMS} in none) AC_DEFINE([GRID_COMMS_NONE],[1],[GRID_COMMS_NONE] ) + comms_type='none' ;; - mpi|mpi-auto) - AC_DEFINE([GRID_COMMS_MPI],[1],[GRID_COMMS_MPI] ) - ;; - mpi3|mpi3-auto) - AC_DEFINE([GRID_COMMS_MPI3],[1],[GRID_COMMS_MPI3] ) - ;; - mpi3l) + mpi3l*) AC_DEFINE([GRID_COMMS_MPI3L],[1],[GRID_COMMS_MPI3L] ) + comms_type='mpi3l' + ;; + mpi3*) + AC_DEFINE([GRID_COMMS_MPI3],[1],[GRID_COMMS_MPI3] ) + comms_type='mpi3' + ;; + mpi*) + AC_DEFINE([GRID_COMMS_MPI],[1],[GRID_COMMS_MPI] ) + comms_type='mpi' ;; shmem) AC_DEFINE([GRID_COMMS_SHMEM],[1],[GRID_COMMS_SHMEM] ) + comms_type='shmem' ;; *) AC_MSG_ERROR([${ac_COMMS} unsupported --enable-comms option]); @@ -282,11 +288,11 @@ case ${ac_COMMS} in ;; esac -AM_CONDITIONAL(BUILD_COMMS_SHMEM,[ test "X${ac_COMMS}X" == "XshmemX" ]) -AM_CONDITIONAL(BUILD_COMMS_MPI,[ test "X${ac_COMMS}X" == "XmpiX" || test "X${ac_COMMS}X" == "Xmpi-autoX" ]) -AM_CONDITIONAL(BUILD_COMMS_MPI3,[ test "X${ac_COMMS}X" == "Xmpi3X"] ) -AM_CONDITIONAL(BUILD_COMMS_MPI3L,[ test "X${ac_COMMS}X" == "Xmpi3lX"] ) -AM_CONDITIONAL(BUILD_COMMS_NONE,[ test "X${ac_COMMS}X" == "XnoneX" ]) +AM_CONDITIONAL(BUILD_COMMS_SHMEM, [ test "${comms_type}X" == "shmemX" ]) +AM_CONDITIONAL(BUILD_COMMS_MPI, [ test "${comms_type}X" == "mpiX" ]) +AM_CONDITIONAL(BUILD_COMMS_MPI3, [ test "${comms_type}X" == "mpi3X" ] ) +AM_CONDITIONAL(BUILD_COMMS_MPI3L, [ test "${comms_type}X" == "mpi3lX" ] ) +AM_CONDITIONAL(BUILD_COMMS_NONE, [ test "${comms_type}X" == "noneX" ]) ############### RNG selection AC_ARG_ENABLE([rng],[AC_HELP_STRING([--enable-rng=ranlux48|mt19937],\ @@ -379,7 +385,7 @@ compiler version : ${ax_cv_gxx_version} ----- BUILD OPTIONS ----------------------------------- SIMD : ${ac_SIMD} Threading : ${ac_openmp} -Communications type : ${ac_COMMS} +Communications type : ${comms_type} Default precision : ${ac_PRECISION} RNG choice : ${ac_RNG} GMP : `if test "x$have_gmp" = xtrue; then echo yes; else echo no; fi` diff --git a/lib/FFT.h b/lib/FFT.h index f9f91ea6..fda43eb8 100644 --- a/lib/FFT.h +++ b/lib/FFT.h @@ -30,8 +30,10 @@ Author: Peter Boyle #define _GRID_FFT_H_ #ifdef HAVE_FFTW -#include +#include #endif + + namespace Grid { template struct FFTW { }; @@ -98,174 +100,198 @@ namespace Grid { #define FFTW_BACKWARD (+1) #endif - class FFT { + class FFT { private: - + GridCartesian *vgrid; GridCartesian *sgrid; - + int Nd; double flops; double flops_call; uint64_t usec; - + std::vector dimensions; std::vector processors; std::vector processor_coor; - + public: - + static const int forward=FFTW_FORWARD; static const int backward=FFTW_BACKWARD; - + double Flops(void) {return flops;} double MFlops(void) {return flops/usec;} - - FFT ( GridCartesian * grid ) : - vgrid(grid), - Nd(grid->_ndimension), - dimensions(grid->_fdimensions), - processors(grid->_processors), - processor_coor(grid->_processor_coor) + + FFT ( GridCartesian * grid ) : + vgrid(grid), + Nd(grid->_ndimension), + dimensions(grid->_fdimensions), + processors(grid->_processors), + processor_coor(grid->_processor_coor) { flops=0; usec =0; std::vector layout(Nd,1); sgrid = new GridCartesian(dimensions,layout,processors); }; - - ~FFT ( void) { - delete sgrid; + + ~FFT ( void) { + delete sgrid; } template - void FFT_dim(Lattice &result,const Lattice &source,int dim, int inverse){ + void FFT_dim_mask(Lattice &result,const Lattice &source,std::vector mask,int sign){ conformable(result._grid,vgrid); conformable(source._grid,vgrid); + Lattice tmp(vgrid); + tmp = source; + for(int d=0;d + void FFT_all_dim(Lattice &result,const Lattice &source,int sign){ + std::vector mask(Nd,1); + FFT_dim_mask(result,source,mask,sign); + } + + + template + void FFT_dim(Lattice &result,const Lattice &source,int dim, int sign){ +#ifndef HAVE_FFTW + assert(0); +#else + conformable(result._grid,vgrid); + conformable(source._grid,vgrid); int L = vgrid->_ldimensions[dim]; int G = vgrid->_fdimensions[dim]; - + std::vector layout(Nd,1); std::vector pencil_gd(vgrid->_fdimensions); - - pencil_gd[dim] = G*processors[dim]; - + + pencil_gd[dim] = G*processors[dim]; + // Pencil global vol LxLxGxLxL per node GridCartesian pencil_g(pencil_gd,layout,processors); - + // Construct pencils typedef typename vobj::scalar_object sobj; typedef typename sobj::scalar_type scalar; + + Lattice pgbuf(&pencil_g); + - Lattice ssource(vgrid); ssource =source; - Lattice pgsource(&pencil_g); - Lattice pgresult(&pencil_g); pgresult=zero; - -#ifndef HAVE_FFTW - assert(0); -#else typedef typename FFTW::FFTW_scalar FFTW_scalar; typedef typename FFTW::FFTW_plan FFTW_plan; - - { - int Ncomp = sizeof(sobj)/sizeof(scalar); - int Nlow = 1; - for(int d=0;d_ldimensions[d]; - } - - int rank = 1; /* 1d transforms */ - int n[] = {G}; /* 1d transforms of length G */ - int howmany = Ncomp; - int odist,idist,istride,ostride; - idist = odist = 1; /* Distance between consecutive FT's */ - istride = ostride = Ncomp*Nlow; /* distance between two elements in the same FT */ - int *inembed = n, *onembed = n; - - - int sign = FFTW_FORWARD; - if (inverse) sign = FFTW_BACKWARD; - - FFTW_plan p; - { - FFTW_scalar *in = (FFTW_scalar *)&pgsource._odata[0]; - FFTW_scalar *out= (FFTW_scalar *)&pgresult._odata[0]; - p = FFTW::fftw_plan_many_dft(rank,n,howmany, - in,inembed, - istride,idist, - out,onembed, - ostride, odist, - sign,FFTW_ESTIMATE); - } - - std::vector lcoor(Nd), gcoor(Nd); - - // Barrel shift and collect global pencil - for(int p=0;plSites();idx++) { - - - sgrid->LocalIndexToLocalCoor(idx,lcoor); - - sobj s; - - peekLocalSite(s,ssource,lcoor); - - lcoor[dim]+=p*L; - - pokeLocalSite(s,pgsource,lcoor); - } - - ssource = Cshift(ssource,dim,L); - } - - // Loop over orthog coords - int NN=pencil_g.lSites(); - GridStopWatch timer; - timer.Start(); - -//PARALLEL_FOR_LOOP - for(int idx=0;idx::fftw_execute_dft(p,in,out); - } - } - - timer.Stop(); - - double add,mul,fma; - FFTW::fftw_flops(p,&add,&mul,&fma); - flops_call = add+mul+2.0*fma; - usec += timer.useconds(); - flops+= flops_call*NN; - int pc = processor_coor[dim]; - for(int idx=0;idxlSites();idx++) { - sgrid->LocalIndexToLocalCoor(idx,lcoor); - gcoor = lcoor; - // extract the result - sobj s; - gcoor[dim] = lcoor[dim]+L*pc; - peekLocalSite(s,pgresult,gcoor); - pokeLocalSite(s,result,lcoor); - } - - FFTW::fftw_destroy_plan(p); + + int Ncomp = sizeof(sobj)/sizeof(scalar); + int Nlow = 1; + for(int d=0;d_ldimensions[d]; } + + int rank = 1; /* 1d transforms */ + int n[] = {G}; /* 1d transforms of length G */ + int howmany = Ncomp; + int odist,idist,istride,ostride; + idist = odist = 1; /* Distance between consecutive FT's */ + istride = ostride = Ncomp*Nlow; /* distance between two elements in the same FT */ + int *inembed = n, *onembed = n; + + scalar div; + if ( sign == backward ) div = 1.0/G; + else if ( sign == forward ) div = 1.0; + else assert(0); + + FFTW_plan p; + { + FFTW_scalar *in = (FFTW_scalar *)&pgbuf._odata[0]; + FFTW_scalar *out= (FFTW_scalar *)&pgbuf._odata[0]; + p = FFTW::fftw_plan_many_dft(rank,n,howmany, + in,inembed, + istride,idist, + out,onembed, + ostride, odist, + sign,FFTW_ESTIMATE); + } + + // Barrel shift and collect global pencil + std::vector lcoor(Nd), gcoor(Nd); + result = source; + for(int p=0;p cbuf(Nd); + sobj s; + + PARALLEL_FOR_LOOP_INTERN + for(int idx=0;idxlSites();idx++) { + sgrid->LocalIndexToLocalCoor(idx,cbuf); + peekLocalSite(s,result,cbuf); + cbuf[dim]+=p*L; + pokeLocalSite(s,pgbuf,cbuf); + } + } + result = Cshift(result,dim,L); + } + + // Loop over orthog coords + int NN=pencil_g.lSites(); + GridStopWatch timer; + timer.Start(); + PARALLEL_REGION + { + std::vector cbuf(Nd); + + PARALLEL_FOR_LOOP_INTERN + for(int idx=0;idx::fftw_execute_dft(p,in,out); + } + } + } + timer.Stop(); + + // performance counting + double add,mul,fma; + FFTW::fftw_flops(p,&add,&mul,&fma); + flops_call = add+mul+2.0*fma; + usec += timer.useconds(); + flops+= flops_call*NN; + + // writing out result + int pc = processor_coor[dim]; + PARALLEL_REGION + { + std::vector clbuf(Nd), cgbuf(Nd); + sobj s; + + PARALLEL_FOR_LOOP_INTERN + for(int idx=0;idxlSites();idx++) { + sgrid->LocalIndexToLocalCoor(idx,clbuf); + cgbuf = clbuf; + cgbuf[dim] = clbuf[dim]+L*pc; + peekLocalSite(s,pgbuf,cgbuf); + s = s * div; + pokeLocalSite(s,result,clbuf); + } + } + + // destroying plan + FFTW::fftw_destroy_plan(p); #endif - - } - }; - - } #endif diff --git a/lib/Grid.h b/lib/Grid.h index 486ee4d3..0c5983f3 100644 --- a/lib/Grid.h +++ b/lib/Grid.h @@ -77,11 +77,10 @@ Author: paboyle #include #include #include -#include -#include - #include +#include +#include #include #include diff --git a/lib/Init.cc b/lib/Init.cc index d15e4bd1..b5ef0a95 100644 --- a/lib/Init.cc +++ b/lib/Init.cc @@ -44,9 +44,33 @@ Author: paboyle #include #include #include +#include +#include + + +#include +#ifdef __APPLE__ +static int +feenableexcept (unsigned int excepts) +{ + static fenv_t fenv; + unsigned int new_excepts = excepts & FE_ALL_EXCEPT, + old_excepts; // previous masks + + if ( fegetenv (&fenv) ) return -1; + old_excepts = fenv.__control & FE_ALL_EXCEPT; + + // unmask + fenv.__control &= ~new_excepts; + fenv.__mxcsr &= ~(new_excepts << 7); + + return ( fesetenv (&fenv) ? -1 : old_excepts ); +} +#endif namespace Grid { + ////////////////////////////////////////////////////// // Convenience functions to access stadard command line arg // driven parallelism controls @@ -234,7 +258,7 @@ void Grid_init(int *argc,char ***argv) std::cout< -#endif + void Grid_debug_handler_init(void) { struct sigaction sa,osa; @@ -402,9 +425,9 @@ void Grid_debug_handler_init(void) sa.sa_flags = SA_SIGINFO; sigaction(SIGSEGV,&sa,NULL); sigaction(SIGTRAP,&sa,NULL); -#ifdef GRID_FPE + feenableexcept( FE_INVALID|FE_OVERFLOW|FE_DIVBYZERO); + sigaction(SIGFPE,&sa,NULL); -#endif } } diff --git a/lib/Init.h b/lib/Init.h index d68cee37..6b70d42d 100644 --- a/lib/Init.h +++ b/lib/Init.h @@ -54,6 +54,7 @@ namespace Grid { void GridCmdOptionCSL(std::string str,std::vector & vec); void GridCmdOptionIntVector(std::string &str,std::vector & vec); + void GridParseLayout(char **argv,int argc, std::vector &latt, std::vector &simd, diff --git a/lib/Log.cc b/lib/Log.cc index d4ac42ee..733bfc58 100644 --- a/lib/Log.cc +++ b/lib/Log.cc @@ -31,8 +31,23 @@ directory /* END LEGAL */ #include +#include + namespace Grid { + std::string demangle(const char* name) { + + int status = -4; // some arbitrary value to eliminate the compiler warning + + // enable c++11 by passing the flag -std=c++11 to g++ + std::unique_ptr res { + abi::__cxa_demangle(name, NULL, NULL, &status), + std::free + }; + + return (status==0) ? res.get() : name ; + } + GridStopWatch Logger::StopWatch; int Logger::timestamp; std::ostream Logger::devnull(0); diff --git a/lib/Log.h b/lib/Log.h index dd3fe927..d7422ca9 100644 --- a/lib/Log.h +++ b/lib/Log.h @@ -144,6 +144,7 @@ extern GridLogger GridLogIterative ; extern GridLogger GridLogIntegrator ; extern Colours GridLogColours; + std::string demangle(const char* name) ; #define _NBACKTRACE (256) extern void * Grid_backtrace_buffer[_NBACKTRACE]; @@ -162,7 +163,7 @@ std::fclose(fp); \ int symbols = backtrace (Grid_backtrace_buffer,_NBACKTRACE);\ char **strings = backtrace_symbols(Grid_backtrace_buffer,symbols);\ for (int i = 0; i < symbols; i++){\ - std::fprintf (fp,"BackTrace Strings: %d %s\n",i, strings[i]); std::fflush(fp); \ + std::fprintf (fp,"BackTrace Strings: %d %s\n",i, demangle(strings[i]).c_str()); std::fflush(fp); \ }\ } #else diff --git a/lib/Simd.h b/lib/Simd.h index 9568c555..adc2849d 100644 --- a/lib/Simd.h +++ b/lib/Simd.h @@ -237,6 +237,18 @@ namespace Grid { stream<<">"; return stream; } + inline std::ostream& operator<< (std::ostream& stream, const vInteger &o){ + int nn=vInteger::Nsimd(); + std::vector > buf(nn); + vstore(o,&buf[0]); + stream<<"<"; + for(int i=0;i"; + return stream; + } } diff --git a/lib/Threads.h b/lib/Threads.h index 08e5d545..2f270b73 100644 --- a/lib/Threads.h +++ b/lib/Threads.h @@ -38,14 +38,19 @@ Author: paboyle #ifdef GRID_OMP #include #ifdef GRID_NUMA -#define PARALLEL_FOR_LOOP _Pragma("omp parallel for schedule(static)") +#define PARALLEL_FOR_LOOP _Pragma("omp parallel for schedule(static)") +#define PARALLEL_FOR_LOOP_INTERN _Pragma("omp for schedule(static)") #else -#define PARALLEL_FOR_LOOP _Pragma("omp parallel for schedule(runtime)") +#define PARALLEL_FOR_LOOP _Pragma("omp parallel for schedule(runtime)") +#define PARALLEL_FOR_LOOP_INTERN _Pragma("omp for schedule(runtime)") #endif #define PARALLEL_NESTED_LOOP2 _Pragma("omp parallel for collapse(2)") +#define PARALLEL_REGION _Pragma("omp parallel") #else -#define PARALLEL_FOR_LOOP +#define PARALLEL_FOR_LOOP +#define PARALLEL_FOR_LOOP_INTERN #define PARALLEL_NESTED_LOOP2 +#define PARALLEL_REGION #endif namespace Grid { diff --git a/lib/communicator/Communicator_mpi3_leader.cc b/lib/communicator/Communicator_mpi3_leader.cc index b6995a44..71f1a913 100644 --- a/lib/communicator/Communicator_mpi3_leader.cc +++ b/lib/communicator/Communicator_mpi3_leader.cc @@ -39,6 +39,10 @@ Author: Peter Boyle /// //////////////////////////////////////////////////////////////////////////////////////////////////////////////// #include +#include +#include +#include + typedef sem_t *Grid_semaphore; #define SEM_INIT(S) S = sem_open(sem_name,0,0600,0); assert ( S != SEM_FAILED ); @@ -48,7 +52,6 @@ typedef sem_t *Grid_semaphore; #include - namespace Grid { enum { COMMAND_ISEND, COMMAND_IRECV, COMMAND_WAITALL }; @@ -91,18 +94,18 @@ public: void SemInit(void) { sprintf(sem_name,"/Grid_mpi3_sem_head_%d",universe_rank); - printf("SEM_NAME: %s \n",sem_name); + // printf("SEM_NAME: %s \n",sem_name); SEM_INIT(sem_head); sprintf(sem_name,"/Grid_mpi3_sem_tail_%d",universe_rank); - printf("SEM_NAME: %s \n",sem_name); + // printf("SEM_NAME: %s \n",sem_name); SEM_INIT(sem_tail); } void SemInitExcl(void) { sprintf(sem_name,"/Grid_mpi3_sem_head_%d",universe_rank); - printf("SEM_INIT_EXCL: %s \n",sem_name); + // printf("SEM_INIT_EXCL: %s \n",sem_name); SEM_INIT_EXCL(sem_head); sprintf(sem_name,"/Grid_mpi3_sem_tail_%d",universe_rank); - printf("SEM_INIT_EXCL: %s \n",sem_name); + // printf("SEM_INIT_EXCL: %s \n",sem_name); SEM_INIT_EXCL(sem_tail); } void WakeUpDMA(void) { @@ -118,7 +121,7 @@ public: SEM_WAIT(sem_tail); }; void EventLoop (void) { - std::cout<< " Entering event loop "<head = state->tail = state->start = 0; base = (uint64_t)MPIoffloadEngine::VerticalShmBufs[0]; int rank; MPI_Comm_rank(_squadron,&rank); diff --git a/lib/fftw/fftw3.h b/lib/fftw/fftw3.h deleted file mode 100644 index 1bf34396..00000000 --- a/lib/fftw/fftw3.h +++ /dev/null @@ -1,412 +0,0 @@ -/* - * Copyright (c) 2003, 2007-14 Matteo Frigo - * Copyright (c) 2003, 2007-14 Massachusetts Institute of Technology - * - * The following statement of license applies *only* to this header file, - * and *not* to the other files distributed with FFTW or derived therefrom: - * - * Redistribution and use in source and binary forms, with or without - * modification, are permitted provided that the following conditions - * are met: - * - * 1. Redistributions of source code must retain the above copyright - * notice, this list of conditions and the following disclaimer. - * - * 2. Redistributions in binary form must reproduce the above copyright - * notice, this list of conditions and the following disclaimer in the - * documentation and/or other materials provided with the distribution. - * - * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS - * OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED - * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE - * ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY - * DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL - * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE - * GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS - * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, - * WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING - * NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS - * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. - */ - -/***************************** NOTE TO USERS ********************************* - * - * THIS IS A HEADER FILE, NOT A MANUAL - * - * If you want to know how to use FFTW, please read the manual, - * online at http://www.fftw.org/doc/ and also included with FFTW. - * For a quick start, see the manual's tutorial section. - * - * (Reading header files to learn how to use a library is a habit - * stemming from code lacking a proper manual. Arguably, it's a - * *bad* habit in most cases, because header files can contain - * interfaces that are not part of the public, stable API.) - * - ****************************************************************************/ - -#ifndef FFTW3_H -#define FFTW3_H - -#include - -#ifdef __cplusplus -extern "C" -{ -#endif /* __cplusplus */ - -/* If is included, use the C99 complex type. Otherwise - define a type bit-compatible with C99 complex */ -#if !defined(FFTW_NO_Complex) && defined(_Complex_I) && defined(complex) && defined(I) -# define FFTW_DEFINE_COMPLEX(R, C) typedef R _Complex C -#else -# define FFTW_DEFINE_COMPLEX(R, C) typedef R C[2] -#endif - -#define FFTW_CONCAT(prefix, name) prefix ## name -#define FFTW_MANGLE_DOUBLE(name) FFTW_CONCAT(fftw_, name) -#define FFTW_MANGLE_FLOAT(name) FFTW_CONCAT(fftwf_, name) -#define FFTW_MANGLE_LONG_DOUBLE(name) FFTW_CONCAT(fftwl_, name) -#define FFTW_MANGLE_QUAD(name) FFTW_CONCAT(fftwq_, name) - -/* IMPORTANT: for Windows compilers, you should add a line - #define FFTW_DLL - here and in kernel/ifftw.h if you are compiling/using FFTW as a - DLL, in order to do the proper importing/exporting, or - alternatively compile with -DFFTW_DLL or the equivalent - command-line flag. This is not necessary under MinGW/Cygwin, where - libtool does the imports/exports automatically. */ -#if defined(FFTW_DLL) && (defined(_WIN32) || defined(__WIN32__)) - /* annoying Windows syntax for shared-library declarations */ -# if defined(COMPILING_FFTW) /* defined in api.h when compiling FFTW */ -# define FFTW_EXTERN extern __declspec(dllexport) -# else /* user is calling FFTW; import symbol */ -# define FFTW_EXTERN extern __declspec(dllimport) -# endif -#else -# define FFTW_EXTERN extern -#endif - -enum fftw_r2r_kind_do_not_use_me { - FFTW_R2HC=0, FFTW_HC2R=1, FFTW_DHT=2, - FFTW_REDFT00=3, FFTW_REDFT01=4, FFTW_REDFT10=5, FFTW_REDFT11=6, - FFTW_RODFT00=7, FFTW_RODFT01=8, FFTW_RODFT10=9, FFTW_RODFT11=10 -}; - -struct fftw_iodim_do_not_use_me { - int n; /* dimension size */ - int is; /* input stride */ - int os; /* output stride */ -}; - -#include /* for ptrdiff_t */ -struct fftw_iodim64_do_not_use_me { - ptrdiff_t n; /* dimension size */ - ptrdiff_t is; /* input stride */ - ptrdiff_t os; /* output stride */ -}; - -typedef void (*fftw_write_char_func_do_not_use_me)(char c, void *); -typedef int (*fftw_read_char_func_do_not_use_me)(void *); - -/* - huge second-order macro that defines prototypes for all API - functions. We expand this macro for each supported precision - - X: name-mangling macro - R: real data type - C: complex data type -*/ - -#define FFTW_DEFINE_API(X, R, C) \ - \ -FFTW_DEFINE_COMPLEX(R, C); \ - \ -typedef struct X(plan_s) *X(plan); \ - \ -typedef struct fftw_iodim_do_not_use_me X(iodim); \ -typedef struct fftw_iodim64_do_not_use_me X(iodim64); \ - \ -typedef enum fftw_r2r_kind_do_not_use_me X(r2r_kind); \ - \ -typedef fftw_write_char_func_do_not_use_me X(write_char_func); \ -typedef fftw_read_char_func_do_not_use_me X(read_char_func); \ - \ -FFTW_EXTERN void X(execute)(const X(plan) p); \ - \ -FFTW_EXTERN X(plan) X(plan_dft)(int rank, const int *n, \ - C *in, C *out, int sign, unsigned flags); \ - \ -FFTW_EXTERN X(plan) X(plan_dft_1d)(int n, C *in, C *out, int sign, \ - unsigned flags); \ -FFTW_EXTERN X(plan) X(plan_dft_2d)(int n0, int n1, \ - C *in, C *out, int sign, unsigned flags); \ -FFTW_EXTERN X(plan) X(plan_dft_3d)(int n0, int n1, int n2, \ - C *in, C *out, int sign, unsigned flags); \ - \ -FFTW_EXTERN X(plan) X(plan_many_dft)(int rank, const int *n, \ - int howmany, \ - C *in, const int *inembed, \ - int istride, int idist, \ - C *out, const int *onembed, \ - int ostride, int odist, \ - int sign, unsigned flags); \ - \ -FFTW_EXTERN X(plan) X(plan_guru_dft)(int rank, const X(iodim) *dims, \ - int howmany_rank, \ - const X(iodim) *howmany_dims, \ - C *in, C *out, \ - int sign, unsigned flags); \ -FFTW_EXTERN X(plan) X(plan_guru_split_dft)(int rank, const X(iodim) *dims, \ - int howmany_rank, \ - const X(iodim) *howmany_dims, \ - R *ri, R *ii, R *ro, R *io, \ - unsigned flags); \ - \ -FFTW_EXTERN X(plan) X(plan_guru64_dft)(int rank, \ - const X(iodim64) *dims, \ - int howmany_rank, \ - const X(iodim64) *howmany_dims, \ - C *in, C *out, \ - int sign, unsigned flags); \ -FFTW_EXTERN X(plan) X(plan_guru64_split_dft)(int rank, \ - const X(iodim64) *dims, \ - int howmany_rank, \ - const X(iodim64) *howmany_dims, \ - R *ri, R *ii, R *ro, R *io, \ - unsigned flags); \ - \ -FFTW_EXTERN void X(execute_dft)(const X(plan) p, C *in, C *out); \ -FFTW_EXTERN void X(execute_split_dft)(const X(plan) p, R *ri, R *ii, \ - R *ro, R *io); \ - \ -FFTW_EXTERN X(plan) X(plan_many_dft_r2c)(int rank, const int *n, \ - int howmany, \ - R *in, const int *inembed, \ - int istride, int idist, \ - C *out, const int *onembed, \ - int ostride, int odist, \ - unsigned flags); \ - \ -FFTW_EXTERN X(plan) X(plan_dft_r2c)(int rank, const int *n, \ - R *in, C *out, unsigned flags); \ - \ -FFTW_EXTERN X(plan) X(plan_dft_r2c_1d)(int n,R *in,C *out,unsigned flags); \ -FFTW_EXTERN X(plan) X(plan_dft_r2c_2d)(int n0, int n1, \ - R *in, C *out, unsigned flags); \ -FFTW_EXTERN X(plan) X(plan_dft_r2c_3d)(int n0, int n1, \ - int n2, \ - R *in, C *out, unsigned flags); \ - \ - \ -FFTW_EXTERN X(plan) X(plan_many_dft_c2r)(int rank, const int *n, \ - int howmany, \ - C *in, const int *inembed, \ - int istride, int idist, \ - R *out, const int *onembed, \ - int ostride, int odist, \ - unsigned flags); \ - \ -FFTW_EXTERN X(plan) X(plan_dft_c2r)(int rank, const int *n, \ - C *in, R *out, unsigned flags); \ - \ -FFTW_EXTERN X(plan) X(plan_dft_c2r_1d)(int n,C *in,R *out,unsigned flags); \ -FFTW_EXTERN X(plan) X(plan_dft_c2r_2d)(int n0, int n1, \ - C *in, R *out, unsigned flags); \ -FFTW_EXTERN X(plan) X(plan_dft_c2r_3d)(int n0, int n1, \ - int n2, \ - C *in, R *out, unsigned flags); \ - \ -FFTW_EXTERN X(plan) X(plan_guru_dft_r2c)(int rank, const X(iodim) *dims, \ - int howmany_rank, \ - const X(iodim) *howmany_dims, \ - R *in, C *out, \ - unsigned flags); \ -FFTW_EXTERN X(plan) X(plan_guru_dft_c2r)(int rank, const X(iodim) *dims, \ - int howmany_rank, \ - const X(iodim) *howmany_dims, \ - C *in, R *out, \ - unsigned flags); \ - \ -FFTW_EXTERN X(plan) X(plan_guru_split_dft_r2c)( \ - int rank, const X(iodim) *dims, \ - int howmany_rank, \ - const X(iodim) *howmany_dims, \ - R *in, R *ro, R *io, \ - unsigned flags); \ -FFTW_EXTERN X(plan) X(plan_guru_split_dft_c2r)( \ - int rank, const X(iodim) *dims, \ - int howmany_rank, \ - const X(iodim) *howmany_dims, \ - R *ri, R *ii, R *out, \ - unsigned flags); \ - \ -FFTW_EXTERN X(plan) X(plan_guru64_dft_r2c)(int rank, \ - const X(iodim64) *dims, \ - int howmany_rank, \ - const X(iodim64) *howmany_dims, \ - R *in, C *out, \ - unsigned flags); \ -FFTW_EXTERN X(plan) X(plan_guru64_dft_c2r)(int rank, \ - const X(iodim64) *dims, \ - int howmany_rank, \ - const X(iodim64) *howmany_dims, \ - C *in, R *out, \ - unsigned flags); \ - \ -FFTW_EXTERN X(plan) X(plan_guru64_split_dft_r2c)( \ - int rank, const X(iodim64) *dims, \ - int howmany_rank, \ - const X(iodim64) *howmany_dims, \ - R *in, R *ro, R *io, \ - unsigned flags); \ -FFTW_EXTERN X(plan) X(plan_guru64_split_dft_c2r)( \ - int rank, const X(iodim64) *dims, \ - int howmany_rank, \ - const X(iodim64) *howmany_dims, \ - R *ri, R *ii, R *out, \ - unsigned flags); \ - \ -FFTW_EXTERN void X(execute_dft_r2c)(const X(plan) p, R *in, C *out); \ -FFTW_EXTERN void X(execute_dft_c2r)(const X(plan) p, C *in, R *out); \ - \ -FFTW_EXTERN void X(execute_split_dft_r2c)(const X(plan) p, \ - R *in, R *ro, R *io); \ -FFTW_EXTERN void X(execute_split_dft_c2r)(const X(plan) p, \ - R *ri, R *ii, R *out); \ - \ -FFTW_EXTERN X(plan) X(plan_many_r2r)(int rank, const int *n, \ - int howmany, \ - R *in, const int *inembed, \ - int istride, int idist, \ - R *out, const int *onembed, \ - int ostride, int odist, \ - const X(r2r_kind) *kind, unsigned flags); \ - \ -FFTW_EXTERN X(plan) X(plan_r2r)(int rank, const int *n, R *in, R *out, \ - const X(r2r_kind) *kind, unsigned flags); \ - \ -FFTW_EXTERN X(plan) X(plan_r2r_1d)(int n, R *in, R *out, \ - X(r2r_kind) kind, unsigned flags); \ -FFTW_EXTERN X(plan) X(plan_r2r_2d)(int n0, int n1, R *in, R *out, \ - X(r2r_kind) kind0, X(r2r_kind) kind1, \ - unsigned flags); \ -FFTW_EXTERN X(plan) X(plan_r2r_3d)(int n0, int n1, int n2, \ - R *in, R *out, X(r2r_kind) kind0, \ - X(r2r_kind) kind1, X(r2r_kind) kind2, \ - unsigned flags); \ - \ -FFTW_EXTERN X(plan) X(plan_guru_r2r)(int rank, const X(iodim) *dims, \ - int howmany_rank, \ - const X(iodim) *howmany_dims, \ - R *in, R *out, \ - const X(r2r_kind) *kind, unsigned flags); \ - \ -FFTW_EXTERN X(plan) X(plan_guru64_r2r)(int rank, const X(iodim64) *dims, \ - int howmany_rank, \ - const X(iodim64) *howmany_dims, \ - R *in, R *out, \ - const X(r2r_kind) *kind, unsigned flags); \ - \ -FFTW_EXTERN void X(execute_r2r)(const X(plan) p, R *in, R *out); \ - \ -FFTW_EXTERN void X(destroy_plan)(X(plan) p); \ -FFTW_EXTERN void X(forget_wisdom)(void); \ -FFTW_EXTERN void X(cleanup)(void); \ - \ -FFTW_EXTERN void X(set_timelimit)(double t); \ - \ -FFTW_EXTERN void X(plan_with_nthreads)(int nthreads); \ -FFTW_EXTERN int X(init_threads)(void); \ -FFTW_EXTERN void X(cleanup_threads)(void); \ - \ -FFTW_EXTERN int X(export_wisdom_to_filename)(const char *filename); \ -FFTW_EXTERN void X(export_wisdom_to_file)(FILE *output_file); \ -FFTW_EXTERN char *X(export_wisdom_to_string)(void); \ -FFTW_EXTERN void X(export_wisdom)(X(write_char_func) write_char, \ - void *data); \ -FFTW_EXTERN int X(import_system_wisdom)(void); \ -FFTW_EXTERN int X(import_wisdom_from_filename)(const char *filename); \ -FFTW_EXTERN int X(import_wisdom_from_file)(FILE *input_file); \ -FFTW_EXTERN int X(import_wisdom_from_string)(const char *input_string); \ -FFTW_EXTERN int X(import_wisdom)(X(read_char_func) read_char, void *data); \ - \ -FFTW_EXTERN void X(fprint_plan)(const X(plan) p, FILE *output_file); \ -FFTW_EXTERN void X(print_plan)(const X(plan) p); \ -FFTW_EXTERN char *X(sprint_plan)(const X(plan) p); \ - \ -FFTW_EXTERN void *X(malloc)(size_t n); \ -FFTW_EXTERN R *X(alloc_real)(size_t n); \ -FFTW_EXTERN C *X(alloc_complex)(size_t n); \ -FFTW_EXTERN void X(free)(void *p); \ - \ -FFTW_EXTERN void X(flops)(const X(plan) p, \ - double *add, double *mul, double *fmas); \ -FFTW_EXTERN double X(estimate_cost)(const X(plan) p); \ -FFTW_EXTERN double X(cost)(const X(plan) p); \ - \ -FFTW_EXTERN int X(alignment_of)(R *p); \ -FFTW_EXTERN const char X(version)[]; \ -FFTW_EXTERN const char X(cc)[]; \ -FFTW_EXTERN const char X(codelet_optim)[]; - - -/* end of FFTW_DEFINE_API macro */ - -FFTW_DEFINE_API(FFTW_MANGLE_DOUBLE, double, fftw_complex) -FFTW_DEFINE_API(FFTW_MANGLE_FLOAT, float, fftwf_complex) -FFTW_DEFINE_API(FFTW_MANGLE_LONG_DOUBLE, long double, fftwl_complex) - -/* __float128 (quad precision) is a gcc extension on i386, x86_64, and ia64 - for gcc >= 4.6 (compiled in FFTW with --enable-quad-precision) */ -#if (__GNUC__ > 4 || (__GNUC__ == 4 && __GNUC_MINOR__ >= 6)) \ - && !(defined(__ICC) || defined(__INTEL_COMPILER)) \ - && (defined(__i386__) || defined(__x86_64__) || defined(__ia64__)) -# if !defined(FFTW_NO_Complex) && defined(_Complex_I) && defined(complex) && defined(I) -/* note: __float128 is a typedef, which is not supported with the _Complex - keyword in gcc, so instead we use this ugly __attribute__ version. - However, we can't simply pass the __attribute__ version to - FFTW_DEFINE_API because the __attribute__ confuses gcc in pointer - types. Hence redefining FFTW_DEFINE_COMPLEX. Ugh. */ -# undef FFTW_DEFINE_COMPLEX -# define FFTW_DEFINE_COMPLEX(R, C) typedef _Complex float __attribute__((mode(TC))) C -# endif -FFTW_DEFINE_API(FFTW_MANGLE_QUAD, __float128, fftwq_complex) -#endif - -#define FFTW_FORWARD (-1) -#define FFTW_BACKWARD (+1) - -#define FFTW_NO_TIMELIMIT (-1.0) - -/* documented flags */ -#define FFTW_MEASURE (0U) -#define FFTW_DESTROY_INPUT (1U << 0) -#define FFTW_UNALIGNED (1U << 1) -#define FFTW_CONSERVE_MEMORY (1U << 2) -#define FFTW_EXHAUSTIVE (1U << 3) /* NO_EXHAUSTIVE is default */ -#define FFTW_PRESERVE_INPUT (1U << 4) /* cancels FFTW_DESTROY_INPUT */ -#define FFTW_PATIENT (1U << 5) /* IMPATIENT is default */ -#define FFTW_ESTIMATE (1U << 6) -#define FFTW_WISDOM_ONLY (1U << 21) - -/* undocumented beyond-guru flags */ -#define FFTW_ESTIMATE_PATIENT (1U << 7) -#define FFTW_BELIEVE_PCOST (1U << 8) -#define FFTW_NO_DFT_R2HC (1U << 9) -#define FFTW_NO_NONTHREADED (1U << 10) -#define FFTW_NO_BUFFERING (1U << 11) -#define FFTW_NO_INDIRECT_OP (1U << 12) -#define FFTW_ALLOW_LARGE_GENERIC (1U << 13) /* NO_LARGE_GENERIC is default */ -#define FFTW_NO_RANK_SPLITS (1U << 14) -#define FFTW_NO_VRANK_SPLITS (1U << 15) -#define FFTW_NO_VRECURSE (1U << 16) -#define FFTW_NO_SIMD (1U << 17) -#define FFTW_NO_SLOW (1U << 18) -#define FFTW_NO_FIXED_RADIX_LARGE_N (1U << 19) -#define FFTW_ALLOW_PRUNING (1U << 20) - -#ifdef __cplusplus -} /* extern "C" */ -#endif /* __cplusplus */ - -#endif /* FFTW3_H */ diff --git a/lib/lattice/Lattice_ET.h b/lib/lattice/Lattice_ET.h index 7ebac99d..1bb83901 100644 --- a/lib/lattice/Lattice_ET.h +++ b/lib/lattice/Lattice_ET.h @@ -261,6 +261,7 @@ GridUnopClass(UnaryExp, exp(a)); GridBinOpClass(BinaryAdd, lhs + rhs); GridBinOpClass(BinarySub, lhs - rhs); GridBinOpClass(BinaryMul, lhs *rhs); +GridBinOpClass(BinaryDiv, lhs /rhs); GridBinOpClass(BinaryAnd, lhs &rhs); GridBinOpClass(BinaryOr, lhs | rhs); @@ -385,6 +386,7 @@ GRID_DEF_UNOP(exp, UnaryExp); GRID_DEF_BINOP(operator+, BinaryAdd); GRID_DEF_BINOP(operator-, BinarySub); GRID_DEF_BINOP(operator*, BinaryMul); +GRID_DEF_BINOP(operator/, BinaryDiv); GRID_DEF_BINOP(operator&, BinaryAnd); GRID_DEF_BINOP(operator|, BinaryOr); diff --git a/lib/lattice/Lattice_base.h b/lib/lattice/Lattice_base.h index e04c1443..e4dc1ca8 100644 --- a/lib/lattice/Lattice_base.h +++ b/lib/lattice/Lattice_base.h @@ -300,17 +300,6 @@ PARALLEL_FOR_LOOP *this = (*this)+r; return *this; } - - strong_inline friend Lattice operator / (const Lattice &lhs,const Lattice &rhs){ - conformable(lhs,rhs); - Lattice ret(lhs._grid); -PARALLEL_FOR_LOOP - for(int ss=0;ssoSites();ss++){ - ret._odata[ss] = lhs._odata[ss]*pow(rhs._odata[ss],-1.0); - } - return ret; - }; - }; // class Lattice template std::ostream& operator<< (std::ostream& stream, const Lattice &o){ diff --git a/lib/lattice/Lattice_rng.h b/lib/lattice/Lattice_rng.h index 2df5e2d1..51cc16ec 100644 --- a/lib/lattice/Lattice_rng.h +++ b/lib/lattice/Lattice_rng.h @@ -294,7 +294,7 @@ namespace Grid { int rank,o_idx,i_idx; _grid->GlobalIndexToGlobalCoor(gidx,gcoor); _grid->GlobalCoorToRankIndex(rank,o_idx,i_idx,gcoor); - + int l_idx=generator_idx(o_idx,i_idx); const int num_rand_seed=16; diff --git a/lib/parallelIO/BinaryIO.h b/lib/parallelIO/BinaryIO.h index 5eddb57d..e2af0545 100644 --- a/lib/parallelIO/BinaryIO.h +++ b/lib/parallelIO/BinaryIO.h @@ -457,7 +457,7 @@ class BinaryIO { // available (how short sighted is that?) ////////////////////////////////////////////////////////// Umu = zero; - static uint32_t csum=0; + static uint32_t csum; csum=0; fobj fileObj; static sobj siteObj; // Static to place in symmetric region for SHMEM diff --git a/lib/qcd/action/fermion/CayleyFermion5D.cc b/lib/qcd/action/fermion/CayleyFermion5D.cc index c4196043..57b047d4 100644 --- a/lib/qcd/action/fermion/CayleyFermion5D.cc +++ b/lib/qcd/action/fermion/CayleyFermion5D.cc @@ -50,6 +50,30 @@ namespace QCD { mass(_mass) { } +template +void CayleyFermion5D::Dminus(const FermionField &psi, FermionField &chi) +{ + int Ls=this->Ls; + FermionField tmp(psi._grid); + + this->DW(psi,tmp,DaggerNo); + + for(int s=0;s +void CayleyFermion5D::DminusDag(const FermionField &psi, FermionField &chi) +{ + int Ls=this->Ls; + FermionField tmp(psi._grid); + + this->DW(psi,tmp,DaggerYes); + + for(int s=0;s void CayleyFermion5D::M5D (const FermionField &psi, FermionField &chi) { diff --git a/lib/qcd/action/fermion/CayleyFermion5D.h b/lib/qcd/action/fermion/CayleyFermion5D.h index 53a93a05..1d8c2b95 100644 --- a/lib/qcd/action/fermion/CayleyFermion5D.h +++ b/lib/qcd/action/fermion/CayleyFermion5D.h @@ -56,6 +56,9 @@ namespace Grid { virtual void M5D (const FermionField &psi, FermionField &chi); virtual void M5Ddag(const FermionField &psi, FermionField &chi); + virtual void Dminus(const FermionField &psi, FermionField &chi); + virtual void DminusDag(const FermionField &psi, FermionField &chi); + ///////////////////////////////////////////////////// // Instantiate different versions depending on Impl ///////////////////////////////////////////////////// @@ -117,6 +120,7 @@ namespace Grid { GridRedBlackCartesian &FourDimRedBlackGrid, RealD _mass,RealD _M5,const ImplParams &p= ImplParams()); + protected: void SetCoefficientsZolotarev(RealD zolohi,Approx::zolotarev_data *zdata,RealD b,RealD c); void SetCoefficientsTanh(Approx::zolotarev_data *zdata,RealD b,RealD c); diff --git a/lib/qcd/action/fermion/DomainWallFermion.h b/lib/qcd/action/fermion/DomainWallFermion.h index b72ca8ec..c0b6b6aa 100644 --- a/lib/qcd/action/fermion/DomainWallFermion.h +++ b/lib/qcd/action/fermion/DomainWallFermion.h @@ -42,6 +42,10 @@ namespace Grid { INHERIT_IMPL_TYPES(Impl); public: + void MomentumSpacePropagator(FermionField &out,const FermionField &in,RealD _m) { + this->MomentumSpacePropagatorHt(out,in,_m); + }; + virtual void Instantiatable(void) {}; // Constructors DomainWallFermion(GaugeField &_Umu, @@ -51,6 +55,7 @@ namespace Grid { GridRedBlackCartesian &FourDimRedBlackGrid, RealD _mass,RealD _M5,const ImplParams &p= ImplParams()) : + CayleyFermion5D(_Umu, FiveDimGrid, FiveDimRedBlackGrid, diff --git a/lib/qcd/action/fermion/FermionOperator.h b/lib/qcd/action/fermion/FermionOperator.h index ea5583eb..742c6e08 100644 --- a/lib/qcd/action/fermion/FermionOperator.h +++ b/lib/qcd/action/fermion/FermionOperator.h @@ -91,6 +91,20 @@ namespace Grid { virtual void Mdiag (const FermionField &in, FermionField &out) { Mooee(in,out);}; // Same as Mooee applied to both CB's virtual void Mdir (const FermionField &in, FermionField &out,int dir,int disp)=0; // case by case Wilson, Clover, Cayley, ContFrac, PartFrac + + virtual void MomentumSpacePropagator(FermionField &out,const FermionField &in,RealD _m) { assert(0);}; + + virtual void FreePropagator(const FermionField &in,FermionField &out,RealD mass) { + FFT theFFT((GridCartesian *) in._grid); + + FermionField in_k(in._grid); + FermionField prop_k(in._grid); + + theFFT.FFT_all_dim(in_k,in,FFT::forward); + this->MomentumSpacePropagator(prop_k,in_k,mass); + theFFT.FFT_all_dim(out,prop_k,FFT::backward); + }; + /////////////////////////////////////////////// // Updates gauge field during HMC /////////////////////////////////////////////// diff --git a/lib/qcd/action/fermion/OverlapWilsonCayleyTanhFermion.h b/lib/qcd/action/fermion/OverlapWilsonCayleyTanhFermion.h index cd23ddd5..9cab0e22 100644 --- a/lib/qcd/action/fermion/OverlapWilsonCayleyTanhFermion.h +++ b/lib/qcd/action/fermion/OverlapWilsonCayleyTanhFermion.h @@ -42,7 +42,11 @@ namespace Grid { INHERIT_IMPL_TYPES(Impl); public: - // Constructors + void MomentumSpacePropagator(FermionField &out,const FermionField &in,RealD _m) { + this->MomentumSpacePropagatorHw(out,in,_m); + }; + + // Constructors OverlapWilsonCayleyTanhFermion(GaugeField &_Umu, GridCartesian &FiveDimGrid, GridRedBlackCartesian &FiveDimRedBlackGrid, diff --git a/lib/qcd/action/fermion/WilsonFermion.cc b/lib/qcd/action/fermion/WilsonFermion.cc index 4bc28bc7..99baa8a0 100644 --- a/lib/qcd/action/fermion/WilsonFermion.cc +++ b/lib/qcd/action/fermion/WilsonFermion.cc @@ -101,6 +101,7 @@ void WilsonFermion::Meooe(const FermionField &in, FermionField &out) { DhopOE(in, out, DaggerNo); } } + template void WilsonFermion::MeooeDag(const FermionField &in, FermionField &out) { if (in.checkerboard == Odd) { @@ -109,32 +110,87 @@ void WilsonFermion::MeooeDag(const FermionField &in, FermionField &out) { DhopOE(in, out, DaggerYes); } } + + template + void WilsonFermion::Mooee(const FermionField &in, FermionField &out) { + out.checkerboard = in.checkerboard; + typename FermionField::scalar_type scal(4.0 + mass); + out = scal * in; + } -template -void WilsonFermion::Mooee(const FermionField &in, FermionField &out) { - out.checkerboard = in.checkerboard; - typename FermionField::scalar_type scal(4.0 + mass); - out = scal * in; -} + template + void WilsonFermion::MooeeDag(const FermionField &in, FermionField &out) { + out.checkerboard = in.checkerboard; + Mooee(in, out); + } -template -void WilsonFermion::MooeeDag(const FermionField &in, FermionField &out) { - out.checkerboard = in.checkerboard; - Mooee(in, out); -} + template + void WilsonFermion::MooeeInv(const FermionField &in, FermionField &out) { + out.checkerboard = in.checkerboard; + out = (1.0/(4.0+mass))*in; + } + + template + void WilsonFermion::MooeeInvDag(const FermionField &in, FermionField &out) { + out.checkerboard = in.checkerboard; + MooeeInv(in,out); + } -template -void WilsonFermion::MooeeInv(const FermionField &in, FermionField &out) { - out.checkerboard = in.checkerboard; - out = (1.0 / (4.0 + mass)) * in; -} + template + void WilsonFermion::MomentumSpacePropagator(FermionField &out, const FermionField &in,RealD _m) { -template -void WilsonFermion::MooeeInvDag(const FermionField &in, - FermionField &out) { - out.checkerboard = in.checkerboard; - MooeeInv(in, out); -} + // what type LatticeComplex + conformable(_grid,out._grid); + + typedef typename FermionField::vector_type vector_type; + typedef typename FermionField::scalar_type ScalComplex; + + typedef Lattice > LatComplex; + + Gamma::GammaMatrix Gmu [] = { + Gamma::GammaX, + Gamma::GammaY, + Gamma::GammaZ, + Gamma::GammaT + }; + + std::vector latt_size = _grid->_fdimensions; + + FermionField num (_grid); num = zero; + LatComplex wilson(_grid); wilson= zero; + LatComplex one (_grid); one = ScalComplex(1.0,0.0); + + LatComplex denom(_grid); denom= zero; + LatComplex kmu(_grid); + ScalComplex ci(0.0,1.0); + // momphase = n * 2pi / L + for(int mu=0;mu, public WilsonFermionStatic { virtual void MooeeInv(const FermionField &in, FermionField &out); virtual void MooeeInvDag(const FermionField &in, FermionField &out); + virtual void MomentumSpacePropagator(FermionField &out,const FermionField &in,RealD _mass) ; + //////////////////////// // Derivative interface //////////////////////// // Interface calls an internal routine - void DhopDeriv(GaugeField &mat, const FermionField &U, const FermionField &V, - int dag); - void DhopDerivOE(GaugeField &mat, const FermionField &U, - const FermionField &V, int dag); - void DhopDerivEO(GaugeField &mat, const FermionField &U, - const FermionField &V, int dag); + void DhopDeriv(GaugeField &mat,const FermionField &U,const FermionField &V,int dag); + void DhopDerivOE(GaugeField &mat,const FermionField &U,const FermionField &V,int dag); + void DhopDerivEO(GaugeField &mat,const FermionField &U,const FermionField &V,int dag); /////////////////////////////////////////////////////////////// // non-hermitian hopping term; half cb or both diff --git a/lib/qcd/action/fermion/WilsonFermion5D.cc b/lib/qcd/action/fermion/WilsonFermion5D.cc index 8c059667..4c2d24bf 100644 --- a/lib/qcd/action/fermion/WilsonFermion5D.cc +++ b/lib/qcd/action/fermion/WilsonFermion5D.cc @@ -482,6 +482,148 @@ void WilsonFermion5D::DW(const FermionField &in, FermionField &out,int dag axpy(out,4.0-M5,in,out); } +template +void WilsonFermion5D::MomentumSpacePropagatorHt(FermionField &out,const FermionField &in, RealD mass) +{ + // what type LatticeComplex + GridBase *_grid = _FourDimGrid; + conformable(_grid,out._grid); + + typedef typename FermionField::vector_type vector_type; + typedef typename FermionField::scalar_type ScalComplex; + typedef iSinglet Tcomplex; + typedef Lattice > LatComplex; + + Gamma::GammaMatrix Gmu [] = { + Gamma::GammaX, + Gamma::GammaY, + Gamma::GammaZ, + Gamma::GammaT + }; + + std::vector latt_size = _grid->_fdimensions; + + + FermionField num (_grid); num = zero; + + LatComplex sk(_grid); sk = zero; + LatComplex sk2(_grid); sk2= zero; + LatComplex W(_grid); W= zero; + LatComplex a(_grid); a= zero; + LatComplex one (_grid); one = ScalComplex(1.0,0.0); + LatComplex denom(_grid); denom= zero; + LatComplex cosha(_grid); + LatComplex kmu(_grid); + LatComplex Wea(_grid); + LatComplex Wema(_grid); + + ScalComplex ci(0.0,1.0); + + for(int mu=0;mu alpha + //////////////////////////////////////////// + cosha = (one + W*W + sk) / (W*2.0); + + // FIXME Need a Lattice acosh + for(int idx=0;idx<_grid->lSites();idx++){ + std::vector lcoor(Nd); + Tcomplex cc; + RealD sgn; + _grid->LocalIndexToLocalCoor(idx,lcoor); + peekLocalSite(cc,cosha,lcoor); + assert((double)real(cc)>=1.0); + assert(fabs((double)imag(cc))<=1.0e-15); + cc = ScalComplex(::acosh(real(cc)),0.0); + pokeLocalSite(cc,a,lcoor); + } + + Wea = ( exp( a) * W ); + Wema= ( exp(-a) * W ); + + num = num + ( one - Wema ) * mass * in; + denom= ( Wea - one ) + mass*mass * (one - Wema); + out = num/denom; +} + +template +void WilsonFermion5D::MomentumSpacePropagatorHw(FermionField &out,const FermionField &in,RealD mass) +{ + Gamma::GammaMatrix Gmu [] = { + Gamma::GammaX, + Gamma::GammaY, + Gamma::GammaZ, + Gamma::GammaT + }; + + GridBase *_grid = _FourDimGrid; + conformable(_grid,out._grid); + + typedef typename FermionField::vector_type vector_type; + typedef typename FermionField::scalar_type ScalComplex; + + typedef Lattice > LatComplex; + + + std::vector latt_size = _grid->_fdimensions; + + LatComplex sk(_grid); sk = zero; + LatComplex sk2(_grid); sk2= zero; + + LatComplex w_k(_grid); w_k= zero; + LatComplex b_k(_grid); b_k= zero; + + LatComplex one (_grid); one = ScalComplex(1.0,0.0); + + FermionField num (_grid); num = zero; + LatComplex denom(_grid); denom= zero; + LatComplex kmu(_grid); + ScalComplex ci(0.0,1.0); + + for(int mu=0;mu directions; - static const std::vector displacements; - const int npoint = 8; - }; - - template - class WilsonFermion5D : public WilsonKernels, public WilsonFermion5DStatic - { - public: - INHERIT_IMPL_TYPES(Impl); - typedef WilsonKernels Kernels; - PmuStat stat; - - void Report(void); - void ZeroCounters(void); - double DhopCalls; - double DhopCommTime; - double DhopComputeTime; - - double DerivCalls; - double DerivCommTime; - double DerivComputeTime; - double DerivDhopComputeTime; - - /////////////////////////////////////////////////////////////// - // Implement the abstract base - /////////////////////////////////////////////////////////////// - GridBase *GaugeGrid(void) { return _FourDimGrid ;} - GridBase *GaugeRedBlackGrid(void) { return _FourDimRedBlackGrid ;} - GridBase *FermionGrid(void) { return _FiveDimGrid;} - GridBase *FermionRedBlackGrid(void) { return _FiveDimRedBlackGrid;} - - // full checkerboard operations; leave unimplemented as abstract for now - virtual RealD M (const FermionField &in, FermionField &out){assert(0); return 0.0;}; - virtual RealD Mdag (const FermionField &in, FermionField &out){assert(0); return 0.0;}; - - // half checkerboard operations; leave unimplemented as abstract for now - virtual void Meooe (const FermionField &in, FermionField &out){assert(0);}; - virtual void Mooee (const FermionField &in, FermionField &out){assert(0);}; - virtual void MooeeInv (const FermionField &in, FermionField &out){assert(0);}; - - virtual void MeooeDag (const FermionField &in, FermionField &out){assert(0);}; - virtual void MooeeDag (const FermionField &in, FermionField &out){assert(0);}; - virtual void MooeeInvDag (const FermionField &in, FermionField &out){assert(0);}; - virtual void Mdir (const FermionField &in, FermionField &out,int dir,int disp){assert(0);}; // case by case Wilson, Clover, Cayley, ContFrac, PartFrac - - // These can be overridden by fancy 5d chiral action - virtual void DhopDeriv (GaugeField &mat,const FermionField &U,const FermionField &V,int dag); - virtual void DhopDerivEO(GaugeField &mat,const FermionField &U,const FermionField &V,int dag); - virtual void DhopDerivOE(GaugeField &mat,const FermionField &U,const FermionField &V,int dag); - - // Implement hopping term non-hermitian hopping term; half cb or both - // Implement s-diagonal DW - void DW (const FermionField &in, FermionField &out,int dag); - void Dhop (const FermionField &in, FermionField &out,int dag); - void DhopOE(const FermionField &in, FermionField &out,int dag); - void DhopEO(const FermionField &in, FermionField &out,int dag); - - // add a DhopComm + //////////////////////////////////////////////////////////////////////////////// + // This is the 4d red black case appropriate to support + // + // parity = (x+y+z+t)|2; + // generalised five dim fermions like mobius, zolotarev etc.. + // + // i.e. even even contains fifth dim hopping term. + // + // [DIFFERS from original CPS red black implementation parity = (x+y+z+t+s)|2 ] + //////////////////////////////////////////////////////////////////////////////// + + class WilsonFermion5DStatic { + public: + // S-direction is INNERMOST and takes no part in the parity. + static const std::vector directions; + static const std::vector displacements; + const int npoint = 8; + }; + + template + class WilsonFermion5D : public WilsonKernels, public WilsonFermion5DStatic + { + public: + INHERIT_IMPL_TYPES(Impl); + typedef WilsonKernels Kernels; + PmuStat stat; + + void Report(void); + void ZeroCounters(void); + double DhopCalls; + double DhopCommTime; + double DhopComputeTime; + + double DerivCalls; + double DerivCommTime; + double DerivComputeTime; + double DerivDhopComputeTime; + + /////////////////////////////////////////////////////////////// + // Implement the abstract base + /////////////////////////////////////////////////////////////// + GridBase *GaugeGrid(void) { return _FourDimGrid ;} + GridBase *GaugeRedBlackGrid(void) { return _FourDimRedBlackGrid ;} + GridBase *FermionGrid(void) { return _FiveDimGrid;} + GridBase *FermionRedBlackGrid(void) { return _FiveDimRedBlackGrid;} + + // full checkerboard operations; leave unimplemented as abstract for now + virtual RealD M (const FermionField &in, FermionField &out){assert(0); return 0.0;}; + virtual RealD Mdag (const FermionField &in, FermionField &out){assert(0); return 0.0;}; + + // half checkerboard operations; leave unimplemented as abstract for now + virtual void Meooe (const FermionField &in, FermionField &out){assert(0);}; + virtual void Mooee (const FermionField &in, FermionField &out){assert(0);}; + virtual void MooeeInv (const FermionField &in, FermionField &out){assert(0);}; + + virtual void MeooeDag (const FermionField &in, FermionField &out){assert(0);}; + virtual void MooeeDag (const FermionField &in, FermionField &out){assert(0);}; + virtual void MooeeInvDag (const FermionField &in, FermionField &out){assert(0);}; + virtual void Mdir (const FermionField &in, FermionField &out,int dir,int disp){assert(0);}; // case by case Wilson, Clover, Cayley, ContFrac, PartFrac + + // These can be overridden by fancy 5d chiral action + virtual void DhopDeriv (GaugeField &mat,const FermionField &U,const FermionField &V,int dag); + virtual void DhopDerivEO(GaugeField &mat,const FermionField &U,const FermionField &V,int dag); + virtual void DhopDerivOE(GaugeField &mat,const FermionField &U,const FermionField &V,int dag); + + void MomentumSpacePropagatorHt(FermionField &out,const FermionField &in,RealD mass) ; + void MomentumSpacePropagatorHw(FermionField &out,const FermionField &in,RealD mass) ; + + // Implement hopping term non-hermitian hopping term; half cb or both + // Implement s-diagonal DW + void DW (const FermionField &in, FermionField &out,int dag); + void Dhop (const FermionField &in, FermionField &out,int dag); + void DhopOE(const FermionField &in, FermionField &out,int dag); + void DhopEO(const FermionField &in, FermionField &out,int dag); + + // add a DhopComm // -- suboptimal interface will presently trigger multiple comms. void DhopDir(const FermionField &in, FermionField &out,int dir,int disp); diff --git a/lib/qcd/action/fermion/WilsonKernels.h b/lib/qcd/action/fermion/WilsonKernels.h index 919f7540..47da2b14 100644 --- a/lib/qcd/action/fermion/WilsonKernels.h +++ b/lib/qcd/action/fermion/WilsonKernels.h @@ -61,14 +61,8 @@ public: switch(Opt) { #ifdef AVX512 case OptInlineAsm: - for (int site = 0; site < Ns; site++) { - for (int s = 0; s < Ls; s++) { - WilsonKernels::DiracOptAsmDhopSite(st,lo,U,buf,sF,sU,Ls,Ns,in,out); - sF++; - } - sU++; - } - break; + WilsonKernels::DiracOptAsmDhopSite(st,lo,U,buf,sF,sU,Ls,Ns,in,out); + break; #endif case OptHandUnroll: for (int site = 0; site < Ns; site++) { @@ -115,13 +109,7 @@ public: switch(Opt) { #ifdef AVX512 case OptInlineAsm: - for (int site = 0; site < Ns; site++) { - for (int s = 0; s < Ls; s++) { - WilsonKernels::DiracOptAsmDhopSiteDag(st,lo,U,buf,sF,sU,Ls,Ns,in,out); - sF++; - } - sU++; - } + WilsonKernels::DiracOptAsmDhopSiteDag(st,lo,U,buf,sF,sU,Ls,Ns,in,out); break; #endif case OptHandUnroll: diff --git a/lib/qcd/utils/LinalgUtils.h b/lib/qcd/utils/LinalgUtils.h index c2cb95ee..e7e6a794 100644 --- a/lib/qcd/utils/LinalgUtils.h +++ b/lib/qcd/utils/LinalgUtils.h @@ -39,8 +39,8 @@ namespace QCD{ //on the 5d (rb4d) checkerboarded lattices //////////////////////////////////////////////////////////////////////// -template -void axpibg5x(Lattice &z,const Lattice &x,RealD a,RealD b) +template +void axpibg5x(Lattice &z,const Lattice &x,Coeff a,Coeff b) { z.checkerboard = x.checkerboard; conformable(x,z); @@ -57,8 +57,8 @@ PARALLEL_FOR_LOOP } } -template -void axpby_ssp(Lattice &z, RealD a,const Lattice &x,RealD b,const Lattice &y,int s,int sp) +template +void axpby_ssp(Lattice &z, Coeff a,const Lattice &x,Coeff b,const Lattice &y,int s,int sp) { z.checkerboard = x.checkerboard; conformable(x,y); @@ -72,8 +72,8 @@ PARALLEL_FOR_LOOP } } -template -void ag5xpby_ssp(Lattice &z,RealD a,const Lattice &x,RealD b,const Lattice &y,int s,int sp) +template +void ag5xpby_ssp(Lattice &z,Coeff a,const Lattice &x,Coeff b,const Lattice &y,int s,int sp) { z.checkerboard = x.checkerboard; conformable(x,y); @@ -90,8 +90,8 @@ PARALLEL_FOR_LOOP } } -template -void axpbg5y_ssp(Lattice &z,RealD a,const Lattice &x,RealD b,const Lattice &y,int s,int sp) +template +void axpbg5y_ssp(Lattice &z,Coeff a,const Lattice &x,Coeff b,const Lattice &y,int s,int sp) { z.checkerboard = x.checkerboard; conformable(x,y); @@ -108,8 +108,8 @@ PARALLEL_FOR_LOOP } } -template -void ag5xpbg5y_ssp(Lattice &z,RealD a,const Lattice &x,RealD b,const Lattice &y,int s,int sp) +template +void ag5xpbg5y_ssp(Lattice &z,Coeff a,const Lattice &x,Coeff b,const Lattice &y,int s,int sp) { z.checkerboard = x.checkerboard; conformable(x,y); @@ -127,8 +127,8 @@ PARALLEL_FOR_LOOP } } -template -void axpby_ssp_pminus(Lattice &z,RealD a,const Lattice &x,RealD b,const Lattice &y,int s,int sp) +template +void axpby_ssp_pminus(Lattice &z,Coeff a,const Lattice &x,Coeff b,const Lattice &y,int s,int sp) { z.checkerboard = x.checkerboard; conformable(x,y); @@ -144,8 +144,8 @@ PARALLEL_FOR_LOOP } } -template -void axpby_ssp_pplus(Lattice &z,RealD a,const Lattice &x,RealD b,const Lattice &y,int s,int sp) +template +void axpby_ssp_pplus(Lattice &z,Coeff a,const Lattice &x,Coeff b,const Lattice &y,int s,int sp) { z.checkerboard = x.checkerboard; conformable(x,y); diff --git a/lib/qcd/utils/SUn.h b/lib/qcd/utils/SUn.h index 20cd1889..914c3bad 100644 --- a/lib/qcd/utils/SUn.h +++ b/lib/qcd/utils/SUn.h @@ -674,6 +674,37 @@ class SU { out += la; } } +/* + add GaugeTrans +*/ + +template + static void GaugeTransform( GaugeField &Umu, GaugeMat &g){ + GridBase *grid = Umu._grid; + conformable(grid,g._grid); + + GaugeMat U(grid); + GaugeMat ag(grid); ag = adj(g); + + for(int mu=0;mu(Umu,mu); + U = g*U*Cshift(ag, mu, 1); + PokeIndex(Umu,U,mu); + } + } + template + static void GaugeTransform( std::vector &U, GaugeMat &g){ + GridBase *grid = g._grid; + GaugeMat ag(grid); ag = adj(g); + for(int mu=0;mu + static void RandomGaugeTransform(GridParallelRNG &pRNG, GaugeField &Umu, GaugeMat &g){ + LieRandomize(pRNG,g,1.0); + GaugeTransform(Umu,g); + } // Projects the algebra components a lattice matrix (of dimension ncol*ncol -1 ) // inverse operation: FundamentalLieAlgebraMatrix @@ -702,23 +733,33 @@ class SU { PokeIndex(out, Umu, mu); } } - static void TepidConfiguration(GridParallelRNG &pRNG, - LatticeGaugeField &out) { - LatticeMatrix Umu(out._grid); - for (int mu = 0; mu < Nd; mu++) { - LieRandomize(pRNG, Umu, 0.01); - PokeIndex(out, Umu, mu); + template + static void TepidConfiguration(GridParallelRNG &pRNG,GaugeField &out){ + typedef typename GaugeField::vector_type vector_type; + typedef iSUnMatrix vMatrixType; + typedef Lattice LatticeMatrixType; + + LatticeMatrixType Umu(out._grid); + for(int mu=0;mu(out,Umu,mu); } } - static void ColdConfiguration(GridParallelRNG &pRNG, LatticeGaugeField &out) { - LatticeMatrix Umu(out._grid); - Umu = 1.0; - for (int mu = 0; mu < Nd; mu++) { - PokeIndex(out, Umu, mu); + template + static void ColdConfiguration(GridParallelRNG &pRNG,GaugeField &out){ + typedef typename GaugeField::vector_type vector_type; + typedef iSUnMatrix vMatrixType; + typedef Lattice LatticeMatrixType; + + LatticeMatrixType Umu(out._grid); + Umu=1.0; + for(int mu=0;mu(out,Umu,mu); } } - static void taProj(const LatticeMatrix &in, LatticeMatrix &out) { + template + static void taProj( const LatticeMatrixType &in, LatticeMatrixType &out){ out = Ta(in); } template diff --git a/lib/qcd/utils/WilsonLoops.h b/lib/qcd/utils/WilsonLoops.h index 10022f50..03d45c07 100644 --- a/lib/qcd/utils/WilsonLoops.h +++ b/lib/qcd/utils/WilsonLoops.h @@ -522,4 +522,4 @@ typedef WilsonLoops SU3WilsonLoops; } } -#endif \ No newline at end of file +#endif diff --git a/lib/simd/Grid_avx.h b/lib/simd/Grid_avx.h index 05f2f1ab..f50eae2b 100644 --- a/lib/simd/Grid_avx.h +++ b/lib/simd/Grid_avx.h @@ -365,6 +365,18 @@ namespace Optimization { } }; + struct Div{ + // Real float + inline __m256 operator()(__m256 a, __m256 b){ + return _mm256_div_ps(a,b); + } + // Real double + inline __m256d operator()(__m256d a, __m256d b){ + return _mm256_div_pd(a,b); + } + }; + + struct Conj{ // Complex single inline __m256 operator()(__m256 in){ @@ -437,14 +449,13 @@ namespace Optimization { }; -#if defined (AVX2) || defined (AVXFMA4) -#define _mm256_alignr_epi32(ret,a,b,n) ret=(__m256) _mm256_alignr_epi8((__m256i)a,(__m256i)b,(n*4)%16) -#define _mm256_alignr_epi64(ret,a,b,n) ret=(__m256d) _mm256_alignr_epi8((__m256i)a,(__m256i)b,(n*8)%16) +#if defined (AVX2) +#define _mm256_alignr_epi32_grid(ret,a,b,n) ret=(__m256) _mm256_alignr_epi8((__m256i)a,(__m256i)b,(n*4)%16) +#define _mm256_alignr_epi64_grid(ret,a,b,n) ret=(__m256d) _mm256_alignr_epi8((__m256i)a,(__m256i)b,(n*8)%16) #endif -#if defined (AVX1) || defined (AVXFMA) - -#define _mm256_alignr_epi32(ret,a,b,n) { \ +#if defined (AVX1) || defined (AVXFMA) +#define _mm256_alignr_epi32_grid(ret,a,b,n) { \ __m128 aa, bb; \ \ aa = _mm256_extractf128_ps(a,1); \ @@ -458,7 +469,7 @@ namespace Optimization { ret = _mm256_insertf128_ps(ret,aa,0); \ } -#define _mm256_alignr_epi64(ret,a,b,n) { \ +#define _mm256_alignr_epi64_grid(ret,a,b,n) { \ __m128d aa, bb; \ \ aa = _mm256_extractf128_pd(a,1); \ @@ -474,19 +485,6 @@ namespace Optimization { #endif - inline std::ostream & operator << (std::ostream& stream, const __m256 a) - { - const float *p=(const float *)&a; - stream<< "{"< 3 ) { - _mm256_alignr_epi32(ret,in,tmp,n); + _mm256_alignr_epi32_grid(ret,in,tmp,n); } else { - _mm256_alignr_epi32(ret,tmp,in,n); + _mm256_alignr_epi32_grid(ret,tmp,in,n); } - // std::cout << " align epi32 n=" < "<< ret < 1 ) { - _mm256_alignr_epi64(ret,in,tmp,n); + _mm256_alignr_epi64_grid(ret,in,tmp,n); } else { - _mm256_alignr_epi64(ret,tmp,in,n); + _mm256_alignr_epi64_grid(ret,tmp,in,n); } - // std::cout << " align epi64 n=" < "<< ret < inline Grid::ComplexF Reduce::operator()(__m256 in){ @@ -631,6 +625,7 @@ namespace Optimization { // Arithmetic operations typedef Optimization::Sum SumSIMD; typedef Optimization::Sub SubSIMD; + typedef Optimization::Div DivSIMD; typedef Optimization::Mult MultSIMD; typedef Optimization::MultComplex MultComplexSIMD; typedef Optimization::Conj ConjSIMD; diff --git a/lib/simd/Grid_avx512.h b/lib/simd/Grid_avx512.h index 136c940e..93e4c9ad 100644 --- a/lib/simd/Grid_avx512.h +++ b/lib/simd/Grid_avx512.h @@ -240,6 +240,17 @@ namespace Optimization { } }; + struct Div{ + // Real float + inline __m512 operator()(__m512 a, __m512 b){ + return _mm512_div_ps(a,b); + } + // Real double + inline __m512d operator()(__m512d a, __m512d b){ + return _mm512_div_pd(a,b); + } + }; + struct Conj{ // Complex single @@ -497,6 +508,7 @@ namespace Optimization { typedef Optimization::Sum SumSIMD; typedef Optimization::Sub SubSIMD; typedef Optimization::Mult MultSIMD; + typedef Optimization::Div DivSIMD; typedef Optimization::MultComplex MultComplexSIMD; typedef Optimization::Conj ConjSIMD; typedef Optimization::TimesMinusI TimesMinusISIMD; diff --git a/lib/simd/Grid_imci.h b/lib/simd/Grid_imci.h index 119a31a3..48cff96d 100644 --- a/lib/simd/Grid_imci.h +++ b/lib/simd/Grid_imci.h @@ -244,6 +244,17 @@ namespace Optimization { } }; + struct Div{ + // Real float + inline __m512 operator()(__m512 a, __m512 b){ + return _mm512_div_ps(a,b); + } + // Real double + inline __m512d operator()(__m512d a, __m512d b){ + return _mm512_div_pd(a,b); + } + }; + struct Conj{ // Complex single @@ -437,6 +448,7 @@ namespace Optimization { // Arithmetic operations typedef Optimization::Sum SumSIMD; typedef Optimization::Sub SubSIMD; + typedef Optimization::Div DivSIMD; typedef Optimization::Mult MultSIMD; typedef Optimization::MultComplex MultComplexSIMD; typedef Optimization::Conj ConjSIMD; diff --git a/lib/simd/Grid_sse4.h b/lib/simd/Grid_sse4.h index a4002373..560eda11 100644 --- a/lib/simd/Grid_sse4.h +++ b/lib/simd/Grid_sse4.h @@ -224,6 +224,18 @@ namespace Optimization { } }; + struct Div{ + // Real float + inline __m128 operator()(__m128 a, __m128 b){ + return _mm_div_ps(a,b); + } + // Real double + inline __m128d operator()(__m128d a, __m128d b){ + return _mm_div_pd(a,b); + } + }; + + struct Conj{ // Complex single inline __m128 operator()(__m128 in){ @@ -372,6 +384,8 @@ namespace Optimization { } } + + ////////////////////////////////////////////////////////////////////////////////////// // Here assign types @@ -398,6 +412,7 @@ namespace Optimization { // Arithmetic operations typedef Optimization::Sum SumSIMD; typedef Optimization::Sub SubSIMD; + typedef Optimization::Div DivSIMD; typedef Optimization::Mult MultSIMD; typedef Optimization::MultComplex MultComplexSIMD; typedef Optimization::Conj ConjSIMD; diff --git a/lib/simd/Grid_vector_types.h b/lib/simd/Grid_vector_types.h index 30576fb6..e1592eef 100644 --- a/lib/simd/Grid_vector_types.h +++ b/lib/simd/Grid_vector_types.h @@ -77,38 +77,24 @@ struct RealPart > { ////////////////////////////////////// // demote a vector to real type ////////////////////////////////////// - // type alias used to simplify the syntax of std::enable_if -template -using Invoke = typename T::type; -template -using EnableIf = Invoke >; -template -using NotEnableIf = Invoke >; +template using Invoke = typename T::type; +template using EnableIf = Invoke >; +template using NotEnableIf = Invoke >; //////////////////////////////////////////////////////// // Check for complexity with type traits -template -struct is_complex : public std::false_type {}; -template <> -struct is_complex > : public std::true_type {}; -template <> -struct is_complex > : public std::true_type {}; +template struct is_complex : public std::false_type {}; +template <> struct is_complex > : public std::true_type {}; +template <> struct is_complex > : public std::true_type {}; -template -using IfReal = Invoke::value, int> >; -template -using IfComplex = Invoke::value, int> >; -template -using IfInteger = Invoke::value, int> >; +template using IfReal = Invoke::value, int> >; +template using IfComplex = Invoke::value, int> >; +template using IfInteger = Invoke::value, int> >; -template -using IfNotReal = - Invoke::value, int> >; -template -using IfNotComplex = Invoke::value, int> >; -template -using IfNotInteger = Invoke::value, int> >; +template using IfNotReal = Invoke::value, int> >; +template using IfNotComplex = Invoke::value, int> >; +template using IfNotInteger = Invoke::value, int> >; //////////////////////////////////////////////////////// // Define the operation templates functors @@ -285,6 +271,20 @@ class Grid_simd { return a * b; } + ////////////////////////////////// + // Divides + ////////////////////////////////// + friend inline Grid_simd operator/(const Scalar_type &a, Grid_simd b) { + Grid_simd va; + vsplat(va, a); + return va / b; + } + friend inline Grid_simd operator/(Grid_simd b, const Scalar_type &a) { + Grid_simd va; + vsplat(va, a); + return b / a; + } + /////////////////////// // Unary negation /////////////////////// @@ -428,7 +428,6 @@ inline void rotate(Grid_simd &ret,Grid_simd b,int nrot) ret.v = Optimization::Rotate::rotate(b.v,2*nrot); } - template inline void vbroadcast(Grid_simd &ret,const Grid_simd &src,int lane){ S* typepun =(S*) &src; @@ -512,7 +511,6 @@ template = 0> inline void vfalse(Grid_simd &ret) { vsplat(ret, 0); } - template inline void zeroit(Grid_simd &z) { vzero(z); @@ -530,7 +528,6 @@ inline void vstream(Grid_simd &out, const Grid_simd &in) { typedef typename S::value_type T; binary((T *)&out.v, in.v, VstreamSIMD()); } - template = 0> inline void vstream(Grid_simd &out, const Grid_simd &in) { out = in; @@ -569,6 +566,34 @@ inline Grid_simd operator*(Grid_simd a, Grid_simd b) { return ret; }; +// Distinguish between complex types and others +template = 0> +inline Grid_simd operator/(Grid_simd a, Grid_simd b) { + typedef Grid_simd simd; + + simd ret; + simd den; + typename simd::conv_t conv; + + ret = a * conjugate(b) ; + den = b * conjugate(b) ; + + + auto real_den = toReal(den); + + ret.v=binary(ret.v, real_den.v, DivSIMD()); + + return ret; +}; + +// Real/Integer types +template = 0> +inline Grid_simd operator/(Grid_simd a, Grid_simd b) { + Grid_simd ret; + ret.v = binary(a.v, b.v, DivSIMD()); + return ret; +}; + /////////////////////// // Conjugate /////////////////////// @@ -582,7 +607,6 @@ template = 0> inline Grid_simd conjugate(const Grid_simd &in) { return in; // for real objects } - // Suppress adj for integer types... // odd; why conjugate above but not adj?? template = 0> inline Grid_simd adj(const Grid_simd &in) { @@ -596,14 +620,12 @@ template = 0> inline void timesMinusI(Grid_simd &ret, const Grid_simd &in) { ret.v = binary(in.v, ret.v, TimesMinusISIMD()); } - template = 0> inline Grid_simd timesMinusI(const Grid_simd &in) { Grid_simd ret; timesMinusI(ret, in); return ret; } - template = 0> inline Grid_simd timesMinusI(const Grid_simd &in) { return in; @@ -616,14 +638,12 @@ template = 0> inline void timesI(Grid_simd &ret, const Grid_simd &in) { ret.v = binary(in.v, ret.v, TimesISIMD()); } - template = 0> inline Grid_simd timesI(const Grid_simd &in) { Grid_simd ret; timesI(ret, in); return ret; } - template = 0> inline Grid_simd timesI(const Grid_simd &in) { return in; diff --git a/lib/tensors/Tensor_arith_mul.h b/lib/tensors/Tensor_arith_mul.h index 917f7373..c24853b7 100644 --- a/lib/tensors/Tensor_arith_mul.h +++ b/lib/tensors/Tensor_arith_mul.h @@ -126,6 +126,36 @@ iVector operator * (const iVector& lhs,const iScalar& r mult(&ret,&lhs,&rhs); return ret; } + +////////////////////////////////////////////////////////////////// +// Divide by scalar +////////////////////////////////////////////////////////////////// +template strong_inline +iScalar operator / (const iScalar& lhs,const iScalar& rhs) +{ + iScalar ret; + ret._internal = lhs._internal/rhs._internal; + return ret; +} +template strong_inline +iVector operator / (const iVector& lhs,const iScalar& rhs) +{ + iVector ret; + for(int i=0;i strong_inline +iMatrix operator / (const iMatrix& lhs,const iScalar& rhs) +{ + iMatrix ret; + for(int i=0;i Make.inc + echo "tests-local: ${TESTLIST}" > Make.inc echo ${PREF}_PROGRAMS = ${TESTLIST} >> Make.inc echo >> Make.inc for f in $TESTS; do diff --git a/tests/core/Test_fft.cc b/tests/core/Test_fft.cc index 2fdaeed2..03c8bc52 100644 --- a/tests/core/Test_fft.cc +++ b/tests/core/Test_fft.cc @@ -1,6 +1,6 @@ /************************************************************************************* - Grid physics library, www.github.com/paboyle/Grid + grid` physics library, www.github.com/paboyle/Grid Source file: ./tests/Test_cshift.cc @@ -46,57 +46,79 @@ int main (int argc, char ** argv) for(int d=0;d p({1,2,3,2}); + std::vector p({1,3,2,3}); one = ComplexD(1.0,0.0); zz = ComplexD(0.0,0.0); ComplexD ci(0.0,1.0); + + std::cout<<"*************************************************"< seeds({1,2,3,4}); + GridSerialRNG sRNG; sRNG.SeedFixedIntegers(seeds); // naughty seeding + GridParallelRNG pRNG(&GRID); + pRNG.SeedFixedIntegers(seeds); + + LatticeGaugeFieldD Umu(&GRID); + + SU3::ColdConfiguration(pRNG,Umu); // Unit gauge + // Umu=zero; + //////////////////////////////////////////////////// + // Wilson test + //////////////////////////////////////////////////// + { + LatticeFermionD src(&GRID); gaussian(pRNG,src); + LatticeFermionD tmp(&GRID); + LatticeFermionD ref(&GRID); + + RealD mass=0.01; + WilsonFermionD Dw(Umu,GRID,RBGRID,mass); + + Dw.M(src,tmp); + + std::cout << "Dw src = " < 1/2 gmu (eip - emip) = i sinp gmu + Kinetic = Kinetic + sin(kmu)*ci*(Gamma(Gmu[mu])*src5_p); + + } + + // NB implicit sum over mu + // + // 1-1/2 Dw = 1 - 1/2 ( eip+emip) + // = - 1/2 (ei - 2 + emi) + // = - 1/4 2 (eih - eimh)(eih - eimh) + // = 2 sink/2 ink/2 = sk2 + + W = one - M5 + sk2; + Kinetic = Kinetic + W * src5_p; + + LatticeCoordinate(scoor,sdir); + + tmp5 = Cshift(src5_p,sdir,+1); + tmp5 = (tmp5 - G5*tmp5)*0.5; + tmp5 = where(scoor==Integer(Ls-1),mass*tmp5,-tmp5); + Kinetic = Kinetic + tmp5; + + tmp5 = Cshift(src5_p,sdir,-1); + tmp5 = (tmp5 + G5*tmp5)*0.5; + tmp5 = where(scoor==Integer(0),mass*tmp5,-tmp5); + Kinetic = Kinetic + tmp5; + + std::cout<<"Momentum space Ddwf "<< norm2(Kinetic)< point(4,0); + src=zero; + SpinColourVectorD ferm; gaussian(sRNG,ferm); + pokeSite(ferm,src,point); + + const int Ls=32; + GridCartesian * FGrid = SpaceTimeGrid::makeFiveDimGrid(Ls,&GRID); + GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,&GRID); + + RealD mass=0.01; + RealD M5 =0.8; + DomainWallFermionD Ddwf(Umu,*FGrid,*FrbGrid,GRID,RBGRID,mass,M5); + + // Momentum space prop + std::cout << " Solving by FFT and Feynman rules" < HermOp(Ddwf); + ConjugateGradient CG(1.0e-16,10000); + CG(HermOp,src5,result5); + + //////////////////////////////////////////////////////////////////////// + // Domain wall physical field propagator + //////////////////////////////////////////////////////////////////////// + /* + psi = chiralProjectMinus(psi_5[0]); + psi += chiralProjectPlus(psi_5[Ls-1]); + */ + ExtractSlice(tmp,result5,0 ,sdir); result4 = (tmp-G5*tmp)*0.5; + ExtractSlice(tmp,result5,Ls-1,sdir); result4 = result4+(tmp+G5*tmp)*0.5; + + std::cout << " Taking difference" < point(4,0); + src=zero; + SpinColourVectorD ferm; gaussian(sRNG,ferm); + pokeSite(ferm,src,point); + + const int Ls=48; + GridCartesian * FGrid = SpaceTimeGrid::makeFiveDimGrid(Ls,&GRID); + GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,&GRID); + + RealD mass=0.01; + RealD M5 =0.8; + + OverlapWilsonCayleyTanhFermionD Dov(Umu,*FGrid,*FrbGrid,GRID,RBGRID,mass,M5,1.0); + + // Momentum space prop + std::cout << " Solving by FFT and Feynman rules" < HermOp(Dov); + ConjugateGradient CG(1.0e-16,10000); + CG(HermOp,src5,result5); + + //////////////////////////////////////////////////////////////////////// + // Domain wall physical field propagator + //////////////////////////////////////////////////////////////////////// + /* + psi = chiralProjectMinus(psi_5[0]); + psi += chiralProjectPlus(psi_5[Ls-1]); + */ + ExtractSlice(tmp,result5,0 ,sdir); result4 = (tmp-G5*tmp)*0.5; + ExtractSlice(tmp,result5,Ls-1,sdir); result4 = result4+(tmp+G5*tmp)*0.5; + + std::cout << " Taking difference" < QEDGimplTypesD; + typedef Photon QEDGaction; + + QEDGaction Maxwell(QEDGaction::FEYNMAN_L); + QEDGaction::GaugeField Prop(&GRID);Prop=zero; + QEDGaction::GaugeField Source(&GRID);Source=zero; + + Maxwell.FreePropagator (Source,Prop); + std::cout << " MaxwellFree propagator\n"; + */ + } Grid_finalize(); } diff --git a/tests/core/Test_fft_gfix.cc b/tests/core/Test_fft_gfix.cc new file mode 100644 index 00000000..6a2868b0 --- /dev/null +++ b/tests/core/Test_fft_gfix.cc @@ -0,0 +1,300 @@ + /************************************************************************************* + + grid` physics library, www.github.com/paboyle/Grid + + Source file: ./tests/Test_cshift.cc + + Copyright (C) 2015 + +Author: Azusa Yamaguchi +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 + +using namespace Grid; +using namespace Grid::QCD; + +template +class FourierAcceleratedGaugeFixer : public Gimpl { + public: + INHERIT_GIMPL_TYPES(Gimpl); + + typedef typename Gimpl::GaugeLinkField GaugeMat; + typedef typename Gimpl::GaugeField GaugeLorentz; + + static void GaugeLinkToLieAlgebraField(const std::vector &U,std::vector &A) { + for(int mu=0;mu &A,GaugeMat &dmuAmu) { + dmuAmu=zero; + for(int mu=0;mu::avgPlaquette(Umu); + RealD org_link_trace=WilsonLoops::linkTrace(Umu); + RealD old_trace = org_link_trace; + RealD trG; + + std::vector U(Nd,grid); + GaugeMat dmuAmu(grid); + + for(int i=0;i(Umu,mu); + //trG = SteepestDescentStep(U,alpha,dmuAmu); + trG = FourierAccelSteepestDescentStep(U,alpha,dmuAmu); + for(int mu=0;mu(Umu,U[mu],mu); + // Monitor progress and convergence test + // infrequently to minimise cost overhead + if ( i %20 == 0 ) { + RealD plaq =WilsonLoops::avgPlaquette(Umu); + RealD link_trace=WilsonLoops::linkTrace(Umu); + + std::cout << GridLogMessage << " Iteration "< &U,RealD & alpha, GaugeMat & dmuAmu) { + GridBase *grid = U[0]._grid; + + std::vector A(Nd,grid); + GaugeMat g(grid); + + GaugeLinkToLieAlgebraField(U,A); + ExpiAlphaDmuAmu(A,g,alpha,dmuAmu); + + + RealD vol = grid->gSites(); + RealD trG = TensorRemove(sum(trace(g))).real()/vol/Nc; + + SU::GaugeTransform(U,g); + + return trG; + } + + static RealD FourierAccelSteepestDescentStep(std::vector &U,RealD & alpha, GaugeMat & dmuAmu) { + + GridBase *grid = U[0]._grid; + + RealD vol = grid->gSites(); + + FFT theFFT((GridCartesian *)grid); + + LatticeComplex Fp(grid); + LatticeComplex psq(grid); psq=zero; + LatticeComplex pmu(grid); + LatticeComplex one(grid); one = ComplexD(1.0,0.0); + + GaugeMat g(grid); + GaugeMat dmuAmu_p(grid); + std::vector A(Nd,grid); + + GaugeLinkToLieAlgebraField(U,A); + + DmuAmu(A,dmuAmu); + + theFFT.FFT_all_dim(dmuAmu_p,dmuAmu,FFT::forward); + + ////////////////////////////////// + // Work out Fp = psq_max/ psq... + ////////////////////////////////// + std::vector latt_size = grid->GlobalDimensions(); + std::vector coor(grid->_ndimension,0); + for(int mu=0;mu::taExp(ciadmam,g); + + RealD trG = TensorRemove(sum(trace(g))).real()/vol/Nc; + + SU::GaugeTransform(U,g); + + return trG; + } + + static void ExpiAlphaDmuAmu(const std::vector &A,GaugeMat &g,RealD & alpha, GaugeMat &dmuAmu) { + GridBase *grid = g._grid; + ComplexD cialpha(0.0,-alpha); + GaugeMat ciadmam(grid); + DmuAmu(A,dmuAmu); + ciadmam = dmuAmu*cialpha; + SU::taExp(ciadmam,g); + } +/* + //////////////////////////////////////////////////////////////// + // NB The FT for fields living on links has an extra phase in it + // Could add these to the FFT class as a later task since this code + // might be reused elsewhere ???? + //////////////////////////////////////////////////////////////// + static void InverseFourierTransformAmu(FFT &theFFT,const std::vector &Ap,std::vector &Ax) { + GridBase * grid = theFFT.Grid(); + std::vector latt_size = grid->GlobalDimensions(); + + ComplexField pmu(grid); + ComplexField pha(grid); + GaugeMat Apha(grid); + + ComplexD ci(0.0,1.0); + + for(int mu=0;mu &Ax,std::vector &Ap) { + GridBase * grid = theFFT.Grid(); + std::vector latt_size = grid->GlobalDimensions(); + + ComplexField pmu(grid); + ComplexField pha(grid); + ComplexD ci(0.0,1.0); + + // Sign convention for FFTW calls: + // A(x)= Sum_p e^ipx A(p) / V + // A(p)= Sum_p e^-ipx A(x) + + for(int mu=0;mu seeds({1,2,3,4}); + + Grid_init(&argc,&argv); + + int threads = GridThread::GetThreads(); + + std::vector latt_size = GridDefaultLatt(); + std::vector simd_layout( { vComplexD::Nsimd(),1,1,1}); + std::vector mpi_layout = GridDefaultMpi(); + + int vol = 1; + for(int d=0;d::avgPlaquette(Umu); + std::cout << " Initial plaquette "<::SteepestDescentGaugeFix(Umu,alpha,10000,1.0e-10, 1.0e-10); + + + plaq=WilsonLoops::avgPlaquette(Umu); + std::cout << " Final plaquette "< +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 + +using namespace Grid; +using namespace Grid::QCD; + +int main (int argc, char ** argv) +{ + Grid_init(&argc,&argv); + + int threads = GridThread::GetThreads(); + std::cout< latt_size ({N,N}); + std::vector simd_layout({vComplexD::Nsimd(),1}); + std::vector mpi_layout ({1,1}); + + int vol = 1; + int nd = latt_size.size(); + for(int d=0;dInteger(N2+W),zz,Charge); + } + + // std::cout << Charge < k(4,&GRID); + LatticeComplexD ksq(&GRID); + + ksq=zero; + for(int mu=0;mu (L/2), k[mu]-L, k[mu]); + + // std::cout << k[mu]< zero_mode(nd,0); + TComplexD Tone = ComplexD(1.0,0.0); + pokeSite(Tone,ksq,zero_mode); + + // std::cout << "Charge\n" << Charge <(mom,mommu,mu); // fourth order exponential approx - parallel_for(auto i=mom.begin();i