portelli/Grid
1
0
Fork 0
mirror of https://github.com/paboyle/Grid.git synced 2025-06-12 20:27:06 +01:00
Code Releases Activity

Merge branch 'master' of https://github.com/paboyle/Grid

Browse Source
This commit is contained in:
Azusa Yamaguchi
2015-06-05 10:04:46 +01:00
parent 94ea84d83f b1b412f63c
commit ac504bea6c
150 changed files with 5950 additions and 2966 deletions
Download Patch File Download Diff File Expand all files Collapse all files

0
lib/.dirstamp Normal file
View File

0
lib/Grid_algorithms.h → lib/Algorithms.h
View File

0
lib/Grid_aligned_allocator.h → lib/AlignedAllocator.h
View File

8
lib/Cartesian.h Normal file
View File

@ -0,0 +1,8 @@
#ifndef GRID_CARTESIAN_H
#define GRID_CARTESIAN_H
#include <cartesian/Cartesian_base.h>
#include <cartesian/Cartesian_full.h>
#include <cartesian/Cartesian_red_black.h>
#endif

2
lib/Grid_communicator.h → lib/Communicator.h
View File

@ -1,6 +1,6 @@
#ifndef GRID_COMMUNICATOR_H
#define GRID_COMMUNICATOR_H
#include <communicator/Grid_communicator_base.h>
#include <communicator/Communicator_base.h>
#endif

4
lib/Grid_comparison.h → lib/Comparison.h
View File

@ -141,7 +141,7 @@ namespace Grid {
}
}
#include <lattice/Grid_lattice_comparison.h>
#include <lattice/Grid_lattice_where.h>
#include <lattice/Lattice_comparison.h>
#include <lattice/Lattice_where.h>
#endif

6
lib/Grid_cshift.h → lib/Cshift.h
View File

@ -1,13 +1,13 @@
#ifndef _GRID_CSHIFT_H_
#define _GRID_CSHIFT_H_
#include <cshift/Grid_cshift_common.h>
#include <cshift/Cshift_common.h>
#ifdef GRID_COMMS_NONE
#include <cshift/Grid_cshift_none.h>
#include <cshift/Cshift_none.h>
#endif
#ifdef GRID_COMMS_MPI
#include <cshift/Grid_cshift_mpi.h>
#include <cshift/Cshift_mpi.h>
#endif
#endif

26
lib/Grid.h
View File

@ -33,7 +33,7 @@
#define strong_inline __attribute__((always_inline)) inline
#include <Grid_config.h>
#include <GridConfig.h>
////////////////////////////////////////////////////////////
// Tunable header includes
@ -46,22 +46,22 @@
#include <malloc.h>
#endif
#include <Grid_aligned_allocator.h>
#include <AlignedAllocator.h>
#include <Grid_simd.h>
#include <Grid_threads.h>
#include <Simd.h>
#include <Threads.h>
#include <Grid_cartesian.h> // subdir aggregate
#include <Grid_math.h> // subdir aggregate
#include <Grid_lattice.h> // subdir aggregate
#include <Grid_comparison.h>
#include <Grid_cshift.h> // subdir aggregate
#include <Grid_stencil.h> // subdir aggregate
#include <Grid_algorithms.h>// subdir aggregate
#include <Communicator.h> // subdir aggregate
#include <Cartesian.h> // subdir aggregate
#include <Tensors.h> // subdir aggregate
#include <Lattice.h> // subdir aggregate
#include <Comparison.h>
#include <Cshift.h> // subdir aggregate
#include <Stencil.h> // subdir aggregate
#include <Algorithms.h>// subdir aggregate
#include <qcd/QCD.h>
#include <parallelIO/GridNerscIO.h>
#include <parallelIO/NerscIO.h>
namespace Grid {

37
lib/Grid_config.h → lib/GridConfig.h
View File

@ -1,5 +1,5 @@
/* lib/Grid_config.h. Generated from Grid_config.h.in by configure. */
/* lib/Grid_config.h.in. Generated from configure.ac by autoheader. */
/* lib/GridConfig.h. Generated from GridConfig.h.in by configure. */
/* lib/GridConfig.h.in. Generated from configure.ac by autoheader. */
/* AVX */
/* #undef AVX1 */
@ -16,6 +16,15 @@
/* GRID_COMMS_NONE */
#define GRID_COMMS_NONE 1
/* Support Altivec instructions */
/* #undef HAVE_ALTIVEC */
/* Support AVX (Advanced Vector Extensions) instructions */
/* #undef HAVE_AVX */
/* Support AVX2 (Advanced Vector Extensions 2) instructions */
/* #undef HAVE_AVX2 */
/* define if the compiler supports basic C++11 syntax */
/* #undef HAVE_CXX11 */
@ -30,6 +39,9 @@
/* Define to 1 if you have the <endian.h> header file. */
#define HAVE_ENDIAN_H 1
/* Support FMA3 (Fused Multiply-Add) instructions */
/* #undef HAVE_FMA */
/* Define to 1 if you have the `gettimeofday' function. */
#define HAVE_GETTIMEOFDAY 1
@ -54,9 +66,30 @@
/* Define to 1 if you have the <memory.h> header file. */
#define HAVE_MEMORY_H 1
/* Support mmx instructions */
#define HAVE_MMX /**/
/* Define to 1 if you have the <mm_malloc.h> header file. */
#define HAVE_MM_MALLOC_H 1
/* Support SSE (Streaming SIMD Extensions) instructions */
#define HAVE_SSE /**/
/* Support SSE2 (Streaming SIMD Extensions 2) instructions */
#define HAVE_SSE2 /**/
/* Support SSE3 (Streaming SIMD Extensions 3) instructions */
#define HAVE_SSE3 /**/
/* Support SSSE4.1 (Streaming SIMD Extensions 4.1) instructions */
#define HAVE_SSE4_1 /**/
/* Support SSSE4.2 (Streaming SIMD Extensions 4.2) instructions */
#define HAVE_SSE4_2 /**/
/* Support SSSE3 (Supplemental Streaming SIMD Extensions 3) instructions */
#define HAVE_SSSE3 /**/
/* Define to 1 if you have the <stdint.h> header file. */
#define HAVE_STDINT_H 1

35
lib/Grid_config.h.in → lib/GridConfig.h.in
View File

@ -1,4 +1,4 @@
/* lib/Grid_config.h.in. Generated from configure.ac by autoheader. */
/* lib/GridConfig.h.in. Generated from configure.ac by autoheader. */
/* AVX */
#undef AVX1
@ -15,6 +15,15 @@
/* GRID_COMMS_NONE */
#undef GRID_COMMS_NONE
/* Support Altivec instructions */
#undef HAVE_ALTIVEC
/* Support AVX (Advanced Vector Extensions) instructions */
#undef HAVE_AVX
/* Support AVX2 (Advanced Vector Extensions 2) instructions */
#undef HAVE_AVX2
/* define if the compiler supports basic C++11 syntax */
#undef HAVE_CXX11
@ -29,6 +38,9 @@
/* Define to 1 if you have the <endian.h> header file. */
#undef HAVE_ENDIAN_H
/* Support FMA3 (Fused Multiply-Add) instructions */
#undef HAVE_FMA
/* Define to 1 if you have the `gettimeofday' function. */
#undef HAVE_GETTIMEOFDAY
@ -53,9 +65,30 @@
/* Define to 1 if you have the <memory.h> header file. */
#undef HAVE_MEMORY_H
/* Support mmx instructions */
#undef HAVE_MMX
/* Define to 1 if you have the <mm_malloc.h> header file. */
#undef HAVE_MM_MALLOC_H
/* Support SSE (Streaming SIMD Extensions) instructions */
#undef HAVE_SSE
/* Support SSE2 (Streaming SIMD Extensions 2) instructions */
#undef HAVE_SSE2
/* Support SSE3 (Streaming SIMD Extensions 3) instructions */
#undef HAVE_SSE3
/* Support SSSE4.1 (Streaming SIMD Extensions 4.1) instructions */
#undef HAVE_SSE4_1
/* Support SSSE4.2 (Streaming SIMD Extensions 4.2) instructions */
#undef HAVE_SSE4_2
/* Support SSSE3 (Supplemental Streaming SIMD Extensions 3) instructions */
#undef HAVE_SSSE3
/* Define to 1 if you have the <stdint.h> header file. */
#undef HAVE_STDINT_H

2
lib/Grid_init.cc → lib/GridInit.cc
View File

@ -143,7 +143,7 @@ void Grid_init(int *argc,char ***argv)
}
if( GridCmdOptionExists(*argv,*argv+*argc,"--dslash-opt") ){
WilsonFermion::HandOptDslash=1;
FiveDimWilsonFermion::HandOptDslash=1;
WilsonFermion5D::HandOptDslash=1;
}
if( GridCmdOptionExists(*argv,*argv+*argc,"--lebesgue") ){
LebesgueOrder::UseLebesgueOrder=1;

8
lib/Grid_cartesian.h
View File

@ -1,8 +0,0 @@
#ifndef GRID_CARTESIAN_H
#define GRID_CARTESIAN_H
#include <cartesian/Grid_cartesian_base.h>
#include <cartesian/Grid_cartesian_full.h>
#include <cartesian/Grid_cartesian_red_black.h>
#endif

16
lib/Grid_math.h
View File

@ -1,16 +0,0 @@
#ifndef GRID_MATH_H
#define GRID_MATH_H
#include <math/Grid_math_traits.h>
#include <math/Grid_math_tensors.h>
#include <math/Grid_math_arith.h>
#include <math/Grid_math_inner.h>
#include <math/Grid_math_outer.h>
#include <math/Grid_math_transpose.h>
#include <math/Grid_math_trace.h>
#include <math/Grid_math_peek.h>
#include <math/Grid_math_poke.h>
#include <math/Grid_math_reality.h>
#endif

2
lib/Grid_lattice.h → lib/Lattice.h
View File

@ -1,6 +1,6 @@
#ifndef GRID_LATTICE_H
#define GRID_LATTICE_H
#include <lattice/Grid_lattice_base.h>
#include <lattice/Lattice_base.h>
#endif

4
lib/Make.inc Normal file
View File

@ -0,0 +1,4 @@
HFILES=./algorithms/approx/bigfloat.h ./algorithms/approx/bigfloat_double.h ./algorithms/approx/Chebyshev.h ./algorithms/approx/Remez.h ./algorithms/approx/Zolotarev.h ./algorithms/iterative/ConjugateGradient.h ./algorithms/iterative/NormalEquations.h ./algorithms/iterative/SchurRedBlack.h ./algorithms/LinearOperator.h ./algorithms/SparseMatrix.h ./Algorithms.h ./AlignedAllocator.h ./cartesian/Cartesian_base.h ./cartesian/Cartesian_full.h ./cartesian/Cartesian_red_black.h ./Cartesian.h ./communicator/Communicator_base.h ./Communicator.h ./Comparison.h ./cshift/Cshift_common.h ./cshift/Cshift_mpi.h ./cshift/Cshift_none.h ./Cshift.h ./Grid.h ./GridConfig.h ./lattice/Lattice_arith.h ./lattice/Lattice_base.h ./lattice/Lattice_comparison.h ./lattice/Lattice_conformable.h ./lattice/Lattice_coordinate.h ./lattice/Lattice_ET.h ./lattice/Lattice_local.h ./lattice/Lattice_overload.h ./lattice/Lattice_peekpoke.h ./lattice/Lattice_reality.h ./lattice/Lattice_reduction.h ./lattice/Lattice_rng.h ./lattice/Lattice_trace.h ./lattice/Lattice_transfer.h ./lattice/Lattice_transpose.h ./lattice/Lattice_where.h ./Lattice.h ./parallelIO/NerscIO.h ./qcd/action/Actions.h ./qcd/action/DiffAction.h ./qcd/action/fermion/CayleyFermion5D.h ./qcd/action/fermion/ContinuedFractionFermion5D.h ./qcd/action/fermion/DomainWallFermion.h ./qcd/action/fermion/FermionOperator.h ./qcd/action/fermion/MobiusFermion.h ./qcd/action/fermion/MobiusZolotarevFermion.h ./qcd/action/fermion/OverlapWilsonCayleyTanhFermion.h ./qcd/action/fermion/OverlapWilsonCayleyZolotarevFermion.h ./qcd/action/fermion/OverlapWilsonContfracTanhFermion.h ./qcd/action/fermion/OverlapWilsonContfracZolotarevFermion.h ./qcd/action/fermion/OverlapWilsonPartialFractionTanhFermion.h ./qcd/action/fermion/OverlapWilsonPartialFractionZolotarevFermion.h ./qcd/action/fermion/PartialFractionFermion5D.h ./qcd/action/fermion/ScaledShamirFermion.h ./qcd/action/fermion/ShamirZolotarevFermion.h ./qcd/action/fermion/WilsonCompressor.h ./qcd/action/fermion/WilsonFermion.h ./qcd/action/fermion/WilsonFermion5D.h ./qcd/action/fermion/WilsonKernels.h ./qcd/Dirac.h ./qcd/LinalgUtils.h ./qcd/QCD.h ./qcd/SpaceTimeGrid.h ./qcd/TwoSpinor.h ./simd/Grid_avx.h ./simd/Grid_avx512.h ./simd/Grid_qpx.h ./simd/Grid_sse4.h ./simd/Grid_vector_types.h ./simd/Old/Grid_vComplexD.h ./simd/Old/Grid_vComplexF.h ./simd/Old/Grid_vInteger.h ./simd/Old/Grid_vRealD.h ./simd/Old/Grid_vRealF.h ./Simd.h ./stencil/Lebesgue.h ./Stencil.h ./tensors/Tensor_arith.h ./tensors/Tensor_arith_add.h ./tensors/Tensor_arith_mac.h ./tensors/Tensor_arith_mul.h ./tensors/Tensor_arith_scalar.h ./tensors/Tensor_arith_sub.h ./tensors/Tensor_class.h ./tensors/Tensor_extract_merge.h ./tensors/Tensor_inner.h ./tensors/Tensor_outer.h ./tensors/Tensor_peek.h ./tensors/Tensor_poke.h ./tensors/Tensor_reality.h ./tensors/Tensor_Ta.h ./tensors/Tensor_trace.h ./tensors/Tensor_traits.h ./tensors/Tensor_transpose.h ./Tensors.h ./Threads.h
CCFILES=./algorithms/approx/Remez.cc ./algorithms/approx/Zolotarev.cc ./GridInit.cc ./qcd/action/fermion/CayleyFermion5D.cc ./qcd/action/fermion/ContinuedFractionFermion5D.cc ./qcd/action/fermion/PartialFractionFermion5D.cc ./qcd/action/fermion/WilsonFermion.cc ./qcd/action/fermion/WilsonFermion5D.cc ./qcd/action/fermion/WilsonKernels.cc ./qcd/action/fermion/WilsonKernelsHand.cc ./qcd/Dirac.cc ./qcd/SpaceTimeGrid.cc ./stencil/Lebesgue.cc ./stencil/Stencil_common.cc

109
lib/Makefile.am
View File

@ -3,115 +3,26 @@ AM_CXXFLAGS = -I$(top_srcdir)/
extra_sources=
if BUILD_COMMS_MPI
extra_sources+=communicator/Grid_communicator_mpi.cc
extra_sources+=communicator/Communicator_mpi.cc
endif
if BUILD_COMMS_NONE
extra_sources+=communicator/Grid_communicator_none.cc
extra_sources+=communicator/Communicator_none.cc
endif
#
# Libraries
#
lib_LIBRARIES = libGrid.a
libGrid_a_SOURCES = \
Grid_init.cc \
stencil/Grid_lebesgue.cc \
stencil/Grid_stencil_common.cc \
algorithms/approx/Zolotarev.cc \
algorithms/approx/Remez.cc \
qcd/action/fermion/FiveDimWilsonFermion.cc\
qcd/action/fermion/WilsonFermion.cc\
qcd/action/fermion/WilsonKernels.cc\
qcd/action/fermion/WilsonKernelsHand.cc\
qcd/Dirac.cc\
$(extra_sources)
include Make.inc
lib_LIBRARIES = libGrid.a
libGrid_a_SOURCES = $(CCFILES) $(extra_sources)
# qcd/action/fermion/PartialFractionFermion5D.cc\ \
#
# Include files
#
nobase_include_HEADERS=\
./algorithms/approx/bigfloat.h\
./algorithms/approx/bigfloat_double.h\
./algorithms/approx/Chebyshev.h\
./algorithms/approx/Remez.h\
./algorithms/approx/Zolotarev.h\
./algorithms/iterative/ConjugateGradient.h\
./algorithms/iterative/NormalEquations.h\
./algorithms/iterative/SchurRedBlack.h\
./algorithms/LinearOperator.h\
./algorithms/SparseMatrix.h\
./cartesian/Grid_cartesian_base.h\
./cartesian/Grid_cartesian_full.h\
./cartesian/Grid_cartesian_red_black.h\
./communicator/Grid_communicator_base.h\
./cshift/Grid_cshift_common.h\
./cshift/Grid_cshift_mpi.h\
./cshift/Grid_cshift_none.h\
./Grid.h\
./Grid_algorithms.h\
./Grid_aligned_allocator.h\
./Grid_cartesian.h\
./Grid_communicator.h\
./Grid_comparison.h\
./Grid_config.h\
./Grid_cshift.h\
./Grid_extract.h\
./Grid_lattice.h\
./Grid_math.h\
./Grid_simd.h\
./Grid_stencil.h\
./Grid_threads.h\
./lattice/Grid_lattice_arith.h\
./lattice/Grid_lattice_base.h\
./lattice/Grid_lattice_comparison.h\
./lattice/Grid_lattice_conformable.h\
./lattice/Grid_lattice_coordinate.h\
./lattice/Grid_lattice_ET.h\
./lattice/Grid_lattice_local.h\
./lattice/Grid_lattice_overload.h\
./lattice/Grid_lattice_peekpoke.h\
./lattice/Grid_lattice_reality.h\
./lattice/Grid_lattice_reduction.h\
./lattice/Grid_lattice_rng.h\
./lattice/Grid_lattice_trace.h\
./lattice/Grid_lattice_transfer.h\
./lattice/Grid_lattice_transpose.h\
./lattice/Grid_lattice_where.h\
./math/Grid_math_arith.h\
./math/Grid_math_arith_add.h\
./math/Grid_math_arith_mac.h\
./math/Grid_math_arith_mul.h\
./math/Grid_math_arith_scalar.h\
./math/Grid_math_arith_sub.h\
./math/Grid_math_inner.h\
./math/Grid_math_outer.h\
./math/Grid_math_peek.h\
./math/Grid_math_poke.h\
./math/Grid_math_reality.h\
./math/Grid_math_tensors.h\
./math/Grid_math_trace.h\
./math/Grid_math_traits.h\
./math/Grid_math_transpose.h\
./parallelIO/GridNerscIO.h\
./qcd/action/Actions.h\
./qcd/action/fermion/FermionAction.h\
./qcd/action/fermion/FiveDimWilsonFermion.h\
./qcd/action/fermion/WilsonCompressor.h\
./qcd/action/fermion/WilsonFermion.h\
./qcd/action/fermion/WilsonKernels.h\
./qcd/Dirac.h\
./qcd/QCD.h\
./qcd/TwoSpinor.h\
./simd/Grid_avx.h\
./simd/Grid_avx512.h\
./simd/Grid_qpx.h\
./simd/Grid_sse4.h\
./simd/Grid_vector_types.h\
./simd/Old/Grid_vComplexD.h\
./simd/Old/Grid_vComplexF.h\
./simd/Old/Grid_vInteger.h\
./simd/Old/Grid_vRealD.h\
./simd/Old/Grid_vRealF.h\
./stencil/Grid_lebesgue.h
nobase_include_HEADERS=$(HFILES)

0
lib/Grid_simd.h → lib/Simd.h
View File

2
lib/Grid_stencil.h → lib/Stencil.h
View File

@ -1,7 +1,7 @@
#ifndef GRID_STENCIL_H
#define GRID_STENCIL_H
#include <stencil/Grid_lebesgue.h> // subdir aggregate
#include <stencil/Lebesgue.h> // subdir aggregate
//////////////////////////////////////////////////////////////////////////////////////////
// Must not lose sight that goal is to be able to construct really efficient

17
lib/Tensors.h Normal file
View File

@ -0,0 +1,17 @@
#ifndef GRID_MATH_H
#define GRID_MATH_H
#include <tensors/Tensor_traits.h>
#include <tensors/Tensor_class.h>
#include <tensors/Tensor_arith.h>
#include <tensors/Tensor_inner.h>
#include <tensors/Tensor_outer.h>
#include <tensors/Tensor_transpose.h>
#include <tensors/Tensor_trace.h>
#include <tensors/Tensor_Ta.h>
#include <tensors/Tensor_peek.h>
#include <tensors/Tensor_poke.h>
#include <tensors/Tensor_reality.h>
#include <tensors/Tensor_extract_merge.h>
#endif

0
lib/Grid_threads.h → lib/Threads.h
View File

32
lib/algorithms/LinearOperator.h
View File

@ -125,39 +125,7 @@ namespace Grid {
};
*/
// Chroma interface defining GaugeAction
/*
template<typename P, typename Q> class GaugeAction
virtual const CreateGaugeState<P,Q>& getCreateState() const = 0;
virtual GaugeState<P,Q>* createState(const Q& q) const
virtual const GaugeBC<P,Q>& getGaugeBC() const
virtual const Set& getSet(void) const = 0;
virtual void deriv(P& result, const Handle< GaugeState<P,Q> >& state) const
virtual Double S(const Handle< GaugeState<P,Q> >& state) const = 0;
class LinearGaugeAction : public GaugeAction< multi1d<LatticeColorMatrix>, multi1d<LatticeColorMatrix> >
typedef multi1d<LatticeColorMatrix> P;
typedef multi1d<LatticeColorMatrix> Q;
virtual void staple(LatticeColorMatrix& result,
const Handle< GaugeState<P,Q> >& state,
int mu, int cb) const = 0;
*/
// Chroma interface defining FermionAction
/*
template<typename T, typename P, typename Q> class FermAct4D : public FermionAction<T,P,Q>
virtual LinearOperator<T>* linOp(Handle< FermState<T,P,Q> > state) const = 0;
virtual LinearOperator<T>* lMdagM(Handle< FermState<T,P,Q> > state) const = 0;
virtual LinOpSystemSolver<T>* invLinOp(Handle< FermState<T,P,Q> > state,
virtual MdagMSystemSolver<T>* invMdagM(Handle< FermState<T,P,Q> > state,
virtual LinOpMultiSystemSolver<T>* mInvLinOp(Handle< FermState<T,P,Q> > state,
virtual MdagMMultiSystemSolver<T>* mInvMdagM(Handle< FermState<T,P,Q> > state,
virtual MdagMMultiSystemSolverAccumulate<T>* mInvMdagMAcc(Handle< FermState<T,P,Q> > state,
virtual SystemSolver<T>* qprop(Handle< FermState<T,P,Q> > state,
class DiffFermAct4D : public FermAct4D<T,P,Q>
virtual DiffLinearOperator<T,Q,P>* linOp(Handle< FermState<T,P,Q> > state) const = 0;
virtual DiffLinearOperator<T,Q,P>* lMdagM(Handle< FermState<T,P,Q> > state) const = 0;
*/
}
#endif

21
lib/algorithms/approx/Zolotarev.cc
View File

@ -58,6 +58,8 @@
/* Compute the partial fraction expansion coefficients (alpha) from the
* factored form */
namespace Grid {
namespace Approx {
static void construct_partfrac(izd *z) {
int dn = z -> dn, dd = z -> dd, type = z -> type;
@ -291,7 +293,7 @@ static void sncndnFK(INTERNAL_PRECISION u, INTERNAL_PRECISION k,
* Set type = 0 for the Zolotarev approximation, which is zero at x = 0, and
* type = 1 for the approximation which is infinite at x = 0. */
zolotarev_data* bfm_zolotarev(PRECISION epsilon, int n, int type) {
zolotarev_data* zolotarev(PRECISION epsilon, int n, int type) {
INTERNAL_PRECISION A, c, cp, kp, ksq, sn, cn, dn, Kp, Kj, z, z0, t, M, F,
l, invlambda, xi, xisq, *tv, s, opl;
int m, czero, ts;
@ -412,7 +414,19 @@ zolotarev_data* bfm_zolotarev(PRECISION epsilon, int n, int type) {
return zd;
}
zolotarev_data* bfm_higham(PRECISION epsilon, int n) {
void zolotarev_free(zolotarev_data *zdata)
{
free(zdata -> a);
free(zdata -> ap);
free(zdata -> alpha);
free(zdata -> beta);
free(zdata -> gamma);
free(zdata);
}
zolotarev_data* higham(PRECISION epsilon, int n) {
INTERNAL_PRECISION A, M, c, cp, z, z0, t, epssq;
int m, czero;
zolotarev_data *zd;
@ -502,6 +516,7 @@ zolotarev_data* bfm_higham(PRECISION epsilon, int n) {
free(d);
return zd;
}
}}
#ifdef TEST
@ -707,4 +722,6 @@ int main(int argc, char** argv) {
return EXIT_SUCCESS;
}
#endif /* TEST */

10
lib/algorithms/approx/Zolotarev.h
View File

@ -1,7 +1,8 @@
/* -*- Mode: C; comment-column: 22; fill-column: 79; -*- */
#ifdef __cplusplus
extern "C" {
namespace Grid {
namespace Approx {
#endif
#define HVERSION Header Time-stamp: <14-OCT-2004 09:26:51.00 adk@MISSCONTRARY>
@ -76,10 +77,11 @@ typedef struct {
* zolotarev_data structure. The arguments must satisfy the constraints that
* epsilon > 0, n > 0, and type = 0 or 1. */
ZOLOTAREV_DATA* bfm_higham(PRECISION epsilon, int n) ;
ZOLOTAREV_DATA* bfm_zolotarev(PRECISION epsilon, int n, int type);
ZOLOTAREV_DATA* higham(PRECISION epsilon, int n) ;
ZOLOTAREV_DATA* zolotarev(PRECISION epsilon, int n, int type);
void zolotarev_free(zolotarev_data *zdata);
#endif
#ifdef __cplusplus
}
}}
#endif

5
lib/algorithms/iterative/NormalEquations.h
View File

@ -4,8 +4,7 @@
namespace Grid {
///////////////////////////////////////////////////////////////////////////////////////////////////////
// Take a matrix and form a Red Black solver calling a Herm solver
// Use of RB info prevents making SchurRedBlackSolve conform to standard interface
// Take a matrix and form an NE solver calling a Herm solver
///////////////////////////////////////////////////////////////////////////////////////////////////////
template<class Field> class NormalEquations : public OperatorFunction<Field>{
private:
@ -15,7 +14,7 @@ namespace Grid {
public:
/////////////////////////////////////////////////////
// Wrap the usual normal equations Schur trick
// Wrap the usual normal equations trick
/////////////////////////////////////////////////////
NormalEquations(SparseMatrixBase<Field> &Matrix, OperatorFunction<Field> &HermitianSolver)
: _Matrix(Matrix), _HermitianSolver(HermitianSolver) {};

3
lib/cartesian/Grid_cartesian_base.h → lib/cartesian/Cartesian_base.h
View File

@ -2,7 +2,6 @@
#define GRID_CARTESIAN_BASE_H
#include <Grid.h>
#include <Grid_communicator.h>
namespace Grid{
@ -21,7 +20,7 @@ public:
// Give Lattice access
template<class object> friend class Lattice;
GridBase(std::vector<int> & processor_grid) : CartesianCommunicator(processor_grid) {};
GridBase(const std::vector<int> & processor_grid) : CartesianCommunicator(processor_grid) {};
// Physics Grid information.

6
lib/cartesian/Grid_cartesian_full.h → lib/cartesian/Cartesian_full.h
View File

@ -27,9 +27,9 @@ public:
virtual int CheckerBoardShift(int source_cb,int dim,int shift, int osite){
return shift;
}
GridCartesian(std::vector<int> &dimensions,
std::vector<int> &simd_layout,
std::vector<int> &processor_grid
GridCartesian(const std::vector<int> &dimensions,
const std::vector<int> &simd_layout,
const std::vector<int> &processor_grid
) : GridBase(processor_grid)
{
///////////////////////

24
lib/cartesian/Grid_cartesian_red_black.h → lib/cartesian/Cartesian_red_black.h
View File

@ -81,28 +81,28 @@ public:
}
};
GridRedBlackCartesian(GridBase *base) : GridRedBlackCartesian(base->_fdimensions,base->_simd_layout,base->_processors) {};
GridRedBlackCartesian(const GridBase *base) : GridRedBlackCartesian(base->_fdimensions,base->_simd_layout,base->_processors) {};
GridRedBlackCartesian(std::vector<int> &dimensions,
std::vector<int> &simd_layout,
std::vector<int> &processor_grid,
std::vector<int> &checker_dim_mask,
GridRedBlackCartesian(const std::vector<int> &dimensions,
const std::vector<int> &simd_layout,
const std::vector<int> &processor_grid,
const std::vector<int> &checker_dim_mask,
int checker_dim
) : GridBase(processor_grid)
{
Init(dimensions,simd_layout,processor_grid,checker_dim_mask,checker_dim);
}
GridRedBlackCartesian(std::vector<int> &dimensions,
std::vector<int> &simd_layout,
std::vector<int> &processor_grid) : GridBase(processor_grid)
GridRedBlackCartesian(const std::vector<int> &dimensions,
const std::vector<int> &simd_layout,
const std::vector<int> &processor_grid) : GridBase(processor_grid)
{
std::vector<int> checker_dim_mask(dimensions.size(),1);
Init(dimensions,simd_layout,processor_grid,checker_dim_mask,0);
}
void Init(std::vector<int> &dimensions,
std::vector<int> &simd_layout,
std::vector<int> &processor_grid,
std::vector<int> &checker_dim_mask,
void Init(const std::vector<int> &dimensions,
const std::vector<int> &simd_layout,
const std::vector<int> &processor_grid,
const std::vector<int> &checker_dim_mask,
int checker_dim)
{
///////////////////////

2
lib/communicator/Grid_communicator_base.h → lib/communicator/Communicator_base.h
View File

@ -27,7 +27,7 @@ class CartesianCommunicator {
#endif
// Constructor
CartesianCommunicator(std::vector<int> &pdimensions_in);
CartesianCommunicator(const std::vector<int> &pdimensions_in);
// Wraps MPI_Cart routines
void ShiftedRanks(int dim,int shift,int & source, int & dest);

2
lib/communicator/Grid_communicator_mpi.cc → lib/communicator/Communicator_mpi.cc
View File

@ -5,7 +5,7 @@ namespace Grid {
// Should error check all MPI calls.
CartesianCommunicator::CartesianCommunicator(std::vector<int> &processors)
CartesianCommunicator::CartesianCommunicator(const std::vector<int> &processors)
{
_ndimension = processors.size();
std::vector<int> periodic(_ndimension,1);

2
lib/communicator/Grid_communicator_none.cc → lib/communicator/Communicator_none.cc
View File

@ -1,7 +1,7 @@
#include "Grid.h"
namespace Grid {
CartesianCommunicator::CartesianCommunicator(std::vector<int> &processors)
CartesianCommunicator::CartesianCommunicator(const std::vector<int> &processors)
{
_processors = processors;
_ndimension = processors.size();

0
lib/cshift/Grid_cshift_common.h → lib/cshift/Cshift_common.h
View File

0
lib/cshift/Grid_cshift_mpi.h → lib/cshift/Cshift_mpi.h
View File

0
lib/cshift/Grid_cshift_none.h → lib/cshift/Cshift_none.h
View File

0
lib/lattice/Grid_lattice_ET.h → lib/lattice/Lattice_ET.h
View File

0
lib/lattice/Grid_lattice_arith.h → lib/lattice/Lattice_arith.h
View File

27
lib/lattice/Grid_lattice_base.h → lib/lattice/Lattice_base.h
View File

@ -283,24 +283,23 @@ PARALLEL_FOR_LOOP
#include <lattice/Grid_lattice_conformable.h>
#include <lattice/Lattice_conformable.h>
#define GRID_LATTICE_EXPRESSION_TEMPLATES
#ifdef GRID_LATTICE_EXPRESSION_TEMPLATES
#include <lattice/Grid_lattice_ET.h>
#include <lattice/Lattice_ET.h>
#else
#include <lattice/Grid_lattice_overload.h>
#include <lattice/Lattice_overload.h>
#endif
#include <lattice/Grid_lattice_arith.h>
#include <lattice/Grid_lattice_trace.h>
#include <lattice/Grid_lattice_transpose.h>
#include <lattice/Grid_lattice_local.h>
#include <lattice/Grid_lattice_reduction.h>
#include <lattice/Grid_lattice_peekpoke.h>
#include <lattice/Grid_lattice_reality.h>
#include <Grid_extract.h>
#include <lattice/Grid_lattice_coordinate.h>
#include <lattice/Grid_lattice_rng.h>
#include <lattice/Grid_lattice_transfer.h>
#include <lattice/Lattice_arith.h>
#include <lattice/Lattice_trace.h>
#include <lattice/Lattice_transpose.h>
#include <lattice/Lattice_local.h>
#include <lattice/Lattice_reduction.h>
#include <lattice/Lattice_peekpoke.h>
#include <lattice/Lattice_reality.h>
#include <lattice/Lattice_coordinate.h>
#include <lattice/Lattice_rng.h>
#include <lattice/Lattice_transfer.h>

0
lib/lattice/Grid_lattice_comparison.h → lib/lattice/Lattice_comparison.h
View File

0
lib/lattice/Grid_lattice_conformable.h → lib/lattice/Lattice_conformable.h
View File

0
lib/lattice/Grid_lattice_coordinate.h → lib/lattice/Lattice_coordinate.h
View File

0
lib/lattice/Grid_lattice_local.h → lib/lattice/Lattice_local.h
View File

0
lib/lattice/Grid_lattice_overload.h → lib/lattice/Lattice_overload.h
View File

0
lib/lattice/Grid_lattice_peekpoke.h → lib/lattice/Lattice_peekpoke.h
View File

11
lib/lattice/Grid_lattice_reality.h → lib/lattice/Lattice_reality.h
View File

@ -48,5 +48,16 @@ PARALLEL_FOR_LOOP
}
template<class vobj> inline auto Ta(const Lattice<vobj> &z) -> Lattice<decltype(Ta(z._odata[0]))>
{
Lattice<decltype(Ta(z._odata[0]))> ret(z._grid);
PARALLEL_FOR_LOOP
for(int ss=0;ss<z._grid->oSites();ss++){
ret._odata[ss] = Ta(z._odata[ss]);
}
return ret;
}
}
#endif

0
lib/lattice/Grid_lattice_reduction.h → lib/lattice/Lattice_reduction.h
View File

0
lib/lattice/Grid_lattice_rng.h → lib/lattice/Lattice_rng.h
View File

0
lib/lattice/Grid_lattice_trace.h → lib/lattice/Lattice_trace.h
View File

0
lib/lattice/Grid_lattice_transfer.h → lib/lattice/Lattice_transfer.h
View File

0
lib/lattice/Grid_lattice_transpose.h → lib/lattice/Lattice_transpose.h
View File

0
lib/lattice/Grid_lattice_where.h → lib/lattice/Lattice_where.h
View File

11
lib/math/Grid_math_arith.h
View File

@ -1,11 +0,0 @@
#ifndef GRID_MATH_ARITH_H
#define GRID_MATH_ARITH_H
#include <math/Grid_math_arith_add.h>
#include <math/Grid_math_arith_sub.h>
#include <math/Grid_math_arith_mac.h>
#include <math/Grid_math_arith_mul.h>
#include <math/Grid_math_arith_scalar.h>
#endif

0
lib/parallelIO/GridNerscIO.h → lib/parallelIO/NerscIO.h
View File

113
lib/qcd/LinalgUtils.h Normal file
View File

@ -0,0 +1,113 @@
#ifndef GRID_QCD_LINALG_UTILS_H
#define GRID_QCD_LINALG_UTILS_H
namespace Grid{
namespace QCD{
////////////////////////////////////////////////////////////////////////
//This file brings additional linear combination assist that is helpful
//to QCD such as chiral projectors and spin matrices applied to one of the inputs.
//These routines support five-D chiral fermions and contain s-subslice indexing
//on the 5d (rb4d) checkerboarded lattices
////////////////////////////////////////////////////////////////////////
template<class vobj>
void axpby_ssp(Lattice<vobj> &z, RealD a,const Lattice<vobj> &x,RealD b,const Lattice<vobj> &y,int s,int sp)
{
z.checkerboard = x.checkerboard;
conformable(x,y);
conformable(x,z);
GridBase *grid=x._grid;
int Ls = grid->_rdimensions[0];
PARALLEL_FOR_LOOP
for(int ss=0;ss<grid->oSites();ss+=Ls){ // adds Ls
vobj tmp = a*x._odata[ss+s]+b*y._odata[ss+sp];
vstream(z._odata[ss+s],tmp);
}
}
template<class vobj>
void ag5xpby_ssp(Lattice<vobj> &z,RealD a,const Lattice<vobj> &x,RealD b,const Lattice<vobj> &y,int s,int sp)
{
z.checkerboard = x.checkerboard;
conformable(x,y);
conformable(x,z);
GridBase *grid=x._grid;
int Ls = grid->_rdimensions[0];
PARALLEL_FOR_LOOP
for(int ss=0;ss<grid->oSites();ss+=Ls){ // adds Ls
vobj tmp;
multGamma5(tmp(),a*x._odata[ss+s]());
tmp = tmp + b*y._odata[ss+sp];
vstream(z._odata[ss+s],tmp);
}
}
template<class vobj>
void axpbg5y_ssp(Lattice<vobj> &z,RealD a,const Lattice<vobj> &x,RealD b,const Lattice<vobj> &y,int s,int sp)
{
z.checkerboard = x.checkerboard;
conformable(x,y);
conformable(x,z);
GridBase *grid=x._grid;
int Ls = grid->_rdimensions[0];
PARALLEL_FOR_LOOP
for(int ss=0;ss<grid->oSites();ss+=Ls){ // adds Ls
vobj tmp;
multGamma5(tmp(),b*y._odata[ss+sp]());
tmp = tmp + a*x._odata[ss+s];
vstream(z._odata[ss+s],tmp);
}
}
template<class vobj>
void ag5xpbg5y_ssp(Lattice<vobj> &z,RealD a,const Lattice<vobj> &x,RealD b,const Lattice<vobj> &y,int s,int sp)
{
z.checkerboard = x.checkerboard;
conformable(x,y);
conformable(x,z);
GridBase *grid=x._grid;
int Ls = grid->_rdimensions[0];
PARALLEL_FOR_LOOP
for(int ss=0;ss<grid->oSites();ss+=Ls){ // adds Ls
vobj tmp1;
vobj tmp2;
tmp1 = a*x._odata[ss+s]+b*y._odata[ss+sp];
multGamma5(tmp2(),tmp1());
vstream(z._odata[ss+s],tmp2);
}
}
template<class vobj>
void axpby_ssp_pminus(Lattice<vobj> &z,RealD a,const Lattice<vobj> &x,RealD b,const Lattice<vobj> &y,int s,int sp)
{
z.checkerboard = x.checkerboard;
conformable(x,y);
conformable(x,z);
GridBase *grid=x._grid;
int Ls = grid->_rdimensions[0];
PARALLEL_FOR_LOOP
for(int ss=0;ss<grid->oSites();ss+=Ls){ // adds Ls
vobj tmp;
spProj5m(tmp,y._odata[ss+sp]);
tmp = a*x._odata[ss+s]+b*tmp;
vstream(z._odata[ss+s],tmp);
}
}
template<class vobj>
void axpby_ssp_pplus(Lattice<vobj> &z,RealD a,const Lattice<vobj> &x,RealD b,const Lattice<vobj> &y,int s,int sp)
{
z.checkerboard = x.checkerboard;
conformable(x,y);
conformable(x,z);
GridBase *grid=x._grid;
int Ls = grid->_rdimensions[0];
PARALLEL_FOR_LOOP
for(int ss=0;ss<grid->oSites();ss+=Ls){ // adds Ls
vobj tmp;
spProj5p(tmp,y._odata[ss+sp]);
tmp = a*x._odata[ss+s]+b*tmp;
vstream(z._odata[ss+s],tmp);
}
}
}}
#endif

2
lib/qcd/QCD.h
View File

@ -307,8 +307,10 @@ namespace QCD {
} //namespace QCD
} // Grid
#include <qcd/SpaceTimeGrid.h>
#include <qcd/Dirac.h>
#include <qcd/TwoSpinor.h>
#include <qcd/LinalgUtils.h>
#include <qcd/action/Actions.h>
#endif

52
lib/qcd/SpaceTimeGrid.cc Normal file
View File

@ -0,0 +1,52 @@
#include <Grid.h>
namespace Grid {
namespace QCD {
/////////////////////////////////////////////////////////////////
// Public interface
/////////////////////////////////////////////////////////////////
GridCartesian *SpaceTimeGrid::makeFourDimGrid(const std::vector<int> & latt,const std::vector<int> &simd,const std::vector<int> &mpi)
{
return new GridCartesian(latt,simd,mpi);
}
GridRedBlackCartesian *SpaceTimeGrid::makeFourDimRedBlackGrid(const GridCartesian *FourDimGrid)
{
return new GridRedBlackCartesian(FourDimGrid);
}
GridCartesian *SpaceTimeGrid::makeFiveDimGrid(int Ls,const GridCartesian *FourDimGrid)
{
int N4=FourDimGrid->_ndimension;
std::vector<int> latt5(1,Ls);
std::vector<int> simd5(1,1);
std::vector<int> mpi5(1,1);
for(int d=0;d<N4;d++){
latt5.push_back(FourDimGrid->_fdimensions[d]);
simd5.push_back(FourDimGrid->_simd_layout[d]);
mpi5.push_back(FourDimGrid->_processors[d]);
}
return new GridCartesian(latt5,simd5,mpi5);
}
GridRedBlackCartesian *SpaceTimeGrid::makeFiveDimRedBlackGrid(int Ls,const GridCartesian *FourDimGrid)
{
int N4=FourDimGrid->_ndimension;
int cbd=1;
std::vector<int> latt5(1,Ls);
std::vector<int> simd5(1,1);
std::vector<int> mpi5(1,1);
std::vector<int> cb5(1,0);
for(int d=0;d<N4;d++){
latt5.push_back(FourDimGrid->_fdimensions[d]);
simd5.push_back(FourDimGrid->_simd_layout[d]);
mpi5.push_back(FourDimGrid->_processors[d]);
cb5.push_back( 1);
}
return new GridRedBlackCartesian(latt5,simd5,mpi5,cb5,cbd);
}
}}

18
lib/qcd/SpaceTimeGrid.h Normal file
View File

@ -0,0 +1,18 @@
#ifndef GRID_QCD_SPACE_TIME_GRID_H
#define GRID_QCD_SPACE_TIME_GRID_H
namespace Grid {
namespace QCD {
class SpaceTimeGrid {
public:
static GridCartesian *makeFourDimGrid(const std::vector<int> & latt,const std::vector<int> &simd,const std::vector<int> &mpi);
static GridRedBlackCartesian *makeFourDimRedBlackGrid (const GridCartesian *FourDimGrid);
static GridCartesian *makeFiveDimGrid (int Ls,const GridCartesian *FourDimGrid);
static GridRedBlackCartesian *makeFiveDimRedBlackGrid(int Ls,const GridCartesian *FourDimGrid);
};
}}
#endif

98
lib/qcd/action/Actions.h
View File

@ -1,10 +1,100 @@
#ifndef GRID_QCD_ACTIONS_H
#define GRID_QCD_ACTIONS_H
#include <qcd/action/fermion/FermionAction.h>
#include <qcd/action/fermion/WilsonCompressor.h>
#include <qcd/action/fermion/WilsonKernels.h>
// Some reorganisation likely required as both Chroma and IroIro
// are separating the concept of the operator from that of action.
//
// The FermAction contains methods to create
//
// * Linear operators (Hermitian and non-hermitian) .. my LinearOperator
// * System solvers (Hermitian and non-hermitian) .. my OperatorFunction
// * MultiShift System solvers (Hermitian and non-hermitian) .. my OperatorFunction
////////////////////////////////////////////
// Abstract base interface
////////////////////////////////////////////
#include <qcd/action/fermion/FermionOperator.h>
////////////////////////////////////////////
// Utility functions
////////////////////////////////////////////
#include <qcd/action/fermion/WilsonCompressor.h> //used by all wilson type fermions
#include <qcd/action/fermion/WilsonKernels.h> //used by all wilson type fermions
////////////////////////////////////////////
// 4D formulations
////////////////////////////////////////////
#include <qcd/action/fermion/WilsonFermion.h>
#include <qcd/action/fermion/FiveDimWilsonFermion.h>
//#include <qcd/action/fermion/CloverFermion.h>
////////////////////////////////////////////
// 5D formulations...
////////////////////////////////////////////
#include <qcd/action/fermion/WilsonFermion5D.h> // used by all 5d overlap types
//////////
// Cayley
//////////
#include <qcd/action/fermion/CayleyFermion5D.h>
#include <qcd/action/fermion/DomainWallFermion.h>
#include <qcd/action/fermion/DomainWallFermion.h>
#include <qcd/action/fermion/MobiusFermion.h>
#include <qcd/action/fermion/ScaledShamirFermion.h>
#include <qcd/action/fermion/OverlapWilsonCayleyTanhFermion.h>
#include <qcd/action/fermion/MobiusZolotarevFermion.h>
#include <qcd/action/fermion/ShamirZolotarevFermion.h>
#include <qcd/action/fermion/OverlapWilsonCayleyZolotarevFermion.h>
//////////////////////
// Continued fraction
//////////////////////
#include <qcd/action/fermion/ContinuedFractionFermion5D.h>
#include <qcd/action/fermion/OverlapWilsonContfracTanhFermion.h>
#include <qcd/action/fermion/OverlapWilsonContfracZolotarevFermion.h>
//////////////////////
// Partial fraction
//////////////////////
#include <qcd/action/fermion/PartialFractionFermion5D.h>
#include <qcd/action/fermion/OverlapWilsonPartialFractionTanhFermion.h>
#include <qcd/action/fermion/OverlapWilsonPartialFractionZolotarevFermion.h>
// Chroma interface defining FermionAction
/*
template<typename T, typename P, typename Q> class FermAct4D : public FermionAction<T,P,Q>
virtual LinearOperator<T>* linOp(Handle< FermState<T,P,Q> > state) const = 0;
virtual LinearOperator<T>* lMdagM(Handle< FermState<T,P,Q> > state) const = 0;
virtual LinOpSystemSolver<T>* invLinOp(Handle< FermState<T,P,Q> > state,
virtual MdagMSystemSolver<T>* invMdagM(Handle< FermState<T,P,Q> > state,
virtual LinOpMultiSystemSolver<T>* mInvLinOp(Handle< FermState<T,P,Q> > state,
virtual MdagMMultiSystemSolver<T>* mInvMdagM(Handle< FermState<T,P,Q> > state,
virtual MdagMMultiSystemSolverAccumulate<T>* mInvMdagMAcc(Handle< FermState<T,P,Q> > state,
virtual SystemSolver<T>* qprop(Handle< FermState<T,P,Q> > state,
class DiffFermAct4D : public FermAct4D<T,P,Q>
virtual DiffLinearOperator<T,Q,P>* linOp(Handle< FermState<T,P,Q> > state) const = 0;
virtual DiffLinearOperator<T,Q,P>* lMdagM(Handle< FermState<T,P,Q> > state) const = 0;
*/
// Chroma interface defining GaugeAction
/*
template<typename P, typename Q> class GaugeAction
virtual const CreateGaugeState<P,Q>& getCreateState() const = 0;
virtual GaugeState<P,Q>* createState(const Q& q) const
virtual const GaugeBC<P,Q>& getGaugeBC() const
virtual const Set& getSet(void) const = 0;
virtual void deriv(P& result, const Handle< GaugeState<P,Q> >& state) const
virtual Double S(const Handle< GaugeState<P,Q> >& state) const = 0;
class LinearGaugeAction : public GaugeAction< multi1d<LatticeColorMatrix>, multi1d<LatticeColorMatrix> >
typedef multi1d<LatticeColorMatrix> P;
*/
#endif

0
lib/qcd/action/fermion/.dirstamp Normal file
View File

346
lib/qcd/action/fermion/CayleyFermion5D.cc Normal file
View File

@ -0,0 +1,346 @@
#include <Grid.h>
namespace Grid {
namespace QCD {
CayleyFermion5D::CayleyFermion5D(LatticeGaugeField &_Umu,
GridCartesian &FiveDimGrid,
GridRedBlackCartesian &FiveDimRedBlackGrid,
GridCartesian &FourDimGrid,
GridRedBlackCartesian &FourDimRedBlackGrid,
RealD _mass,RealD _M5) :
WilsonFermion5D(_Umu,
FiveDimGrid,
FiveDimRedBlackGrid,
FourDimGrid,
FourDimRedBlackGrid,_M5),
mass(_mass)
{
}
// override multiply
RealD CayleyFermion5D::M (const LatticeFermion &psi, LatticeFermion &chi)
{
LatticeFermion Din(psi._grid);
// Assemble Din
for(int s=0;s<Ls;s++){
if ( s==0 ) {
// Din = bs psi[s] + cs[s] psi[s+1}
axpby_ssp_pminus(Din,bs[s],psi,cs[s],psi,s,s+1);
// Din+= -mass*cs[s] psi[s+1}
axpby_ssp_pplus (Din,1.0,Din,-mass*cs[s],psi,s,Ls-1);
} else if ( s==(Ls-1)) {
axpby_ssp_pminus(Din,bs[s],psi,-mass*cs[s],psi,s,0);
axpby_ssp_pplus (Din,1.0,Din,cs[s],psi,s,s-1);
} else {
axpby_ssp_pminus(Din,bs[s],psi,cs[s],psi,s,s+1);
axpby_ssp_pplus(Din,1.0,Din,cs[s],psi,s,s-1);
}
}
DW(Din,chi,DaggerNo);
// ((b D_W + D_w hop terms +1) on s-diag
axpby(chi,1.0,1.0,chi,psi);
for(int s=0;s<Ls;s++){
if ( s==0 ){
axpby_ssp_pminus(chi,1.0,chi,-1.0,psi,s,s+1);
axpby_ssp_pplus (chi,1.0,chi,mass,psi,s,Ls-1);
} else if ( s==(Ls-1)) {
axpby_ssp_pminus(chi,1.0,chi,mass,psi,s,0);
axpby_ssp_pplus (chi,1.0,chi,-1.0,psi,s,s-1);
} else {
axpby_ssp_pminus(chi,1.0,chi,-1.0,psi,s,s+1);
axpby_ssp_pplus (chi,1.0,chi,-1.0,psi,s,s-1);
}
}
return norm2(chi);
}
RealD CayleyFermion5D::Mdag (const LatticeFermion &psi, LatticeFermion &chi)
{
// Under adjoint
//D1+ D1- P- -> D1+^dag P+ D2-^dag
//D2- P+ D2+ P-D1-^dag D2+dag
LatticeFermion Din(psi._grid);
// Apply Dw
DW(psi,Din,DaggerYes);
for(int s=0;s<Ls;s++){
// Collect the terms in DW
// Chi = bs Din[s] + cs[s] Din[s+1}
// Chi+= -mass*cs[s] psi[s+1}
if ( s==0 ) {
axpby_ssp_pplus (chi,bs[s],Din,cs[s+1],Din,s,s+1);
axpby_ssp_pminus(chi,1.0,chi,-mass*cs[Ls-1],Din,s,Ls-1);
} else if ( s==(Ls-1)) {
axpby_ssp_pplus (chi,bs[s],Din,-mass*cs[0],Din,s,0);
axpby_ssp_pminus(chi,1.0,chi,cs[s-1],Din,s,s-1);
} else {
axpby_ssp_pplus (chi,bs[s],Din,cs[s+1],Din,s,s+1);
axpby_ssp_pminus(chi,1.0,chi,cs[s-1],Din,s,s-1);
}
// Collect the terms indept of DW
if ( s==0 ){
axpby_ssp_pplus (chi,1.0,chi,-1.0,psi,s,s+1);
axpby_ssp_pminus(chi,1.0,chi,mass,psi,s,Ls-1);
} else if ( s==(Ls-1)) {
axpby_ssp_pplus (chi,1.0,chi,mass,psi,s,0);
axpby_ssp_pminus(chi,1.0,chi,-1.0,psi,s,s-1);
} else {
axpby_ssp_pplus(chi,1.0,chi,-1.0,psi,s,s+1);
axpby_ssp_pminus(chi,1.0,chi,-1.0,psi,s,s-1);
}
}
// ((b D_W + D_w hop terms +1) on s-diag
axpby (chi,1.0,1.0,chi,psi);
return norm2(chi);
}
// half checkerboard operations
void CayleyFermion5D::Meooe (const LatticeFermion &psi, LatticeFermion &chi)
{
LatticeFermion tmp(psi._grid);
// Assemble the 5d matrix
for(int s=0;s<Ls;s++){
if ( s==0 ) {
// tmp = bs psi[s] + cs[s] psi[s+1}
// tmp+= -mass*cs[s] psi[s+1}
axpby_ssp_pminus(tmp,beo[s],psi,-ceo[s],psi ,s, s+1);
axpby_ssp_pplus(tmp,1.0,tmp,mass*ceo[s],psi,s,Ls-1);
} else if ( s==(Ls-1)) {
axpby_ssp_pminus(tmp,beo[s],psi,mass*ceo[s],psi,s,0);
axpby_ssp_pplus(tmp,1.0,tmp,-ceo[s],psi,s,s-1);
} else {
axpby_ssp_pminus(tmp,beo[s],psi,-ceo[s],psi,s,s+1);
axpby_ssp_pplus (tmp,1.0,tmp,-ceo[s],psi,s,s-1);
}
}
// Apply 4d dslash
if ( psi.checkerboard == Odd ) {
DhopEO(tmp,chi,DaggerNo);
} else {
DhopOE(tmp,chi,DaggerNo);
}
}
void CayleyFermion5D::MeooeDag (const LatticeFermion &psi, LatticeFermion &chi)
{
LatticeFermion tmp(psi._grid);
// Apply 4d dslash
if ( psi.checkerboard == Odd ) {
DhopEO(psi,tmp,DaggerYes);
} else {
DhopOE(psi,tmp,DaggerYes);
}
// Assemble the 5d matrix
for(int s=0;s<Ls;s++){
if ( s==0 ) {
axpby_ssp_pplus(chi,beo[s],tmp, -ceo[s+1] ,tmp,s,s+1);
axpby_ssp_pminus(chi, 1.0,chi,mass*ceo[Ls-1],tmp,s,Ls-1);
} else if ( s==(Ls-1)) {
axpby_ssp_pplus(chi,beo[s],tmp,mass*ceo[0],tmp,s,0);
axpby_ssp_pminus(chi,1.0,chi,-ceo[s-1],tmp,s,s-1);
} else {
axpby_ssp_pplus(chi,beo[s],tmp,-ceo[s+1],tmp,s,s+1);
axpby_ssp_pminus(chi,1.0 ,chi,-ceo[s-1],tmp,s,s-1);
}
}
}
void CayleyFermion5D::Mooee (const LatticeFermion &psi, LatticeFermion &chi)
{
for (int s=0;s<Ls;s++){
if ( s==0 ) {
axpby_ssp_pminus(chi,bee[s],psi ,-cee[s],psi,s,s+1);
axpby_ssp_pplus (chi,1.0,chi,mass*cee[s],psi,s,Ls-1);
} else if ( s==(Ls-1)) {
axpby_ssp_pminus(chi,bee[s],psi,mass*cee[s],psi,s,0);
axpby_ssp_pplus (chi,1.0,chi,-cee[s],psi,s,s-1);
} else {
axpby_ssp_pminus(chi,bee[s],psi,-cee[s],psi,s,s+1);
axpby_ssp_pplus (chi,1.0,chi,-cee[s],psi,s,s-1);
}
}
}
void CayleyFermion5D::MooeeDag (const LatticeFermion &psi, LatticeFermion &chi)
{
for (int s=0;s<Ls;s++){
// Assemble the 5d matrix
if ( s==0 ) {
axpby_ssp_pplus(chi,bee[s],psi,-cee[s+1] ,psi,s,s+1);
axpby_ssp_pminus(chi,1.0,chi,mass*cee[Ls-1],psi,s,Ls-1);
} else if ( s==(Ls-1)) {
axpby_ssp_pplus(chi,bee[s],psi,mass*cee[0],psi,s,0);
axpby_ssp_pminus(chi,1.0,chi,-cee[s-1],psi,s,s-1);
} else {
axpby_ssp_pplus(chi,bee[s],psi,-cee[s+1],psi,s,s+1);
axpby_ssp_pminus(chi,1.0 ,chi,-cee[s-1],psi,s,s-1);
}
}
}
void CayleyFermion5D::MooeeInv (const LatticeFermion &psi, LatticeFermion &chi)
{
// Apply (L^{\prime})^{-1}
axpby_ssp (chi,1.0,psi, 0.0,psi,0,0); // chi[0]=psi[0]
for (int s=1;s<Ls;s++){
axpby_ssp_pplus(chi,1.0,psi,-lee[s-1],chi,s,s-1);// recursion Psi[s] -lee P_+ chi[s-1]
}
// L_m^{-1}
for (int s=0;s<Ls-1;s++){ // Chi[ee] = 1 - sum[s<Ls-1] -leem[s]P_- chi
axpby_ssp_pminus(chi,1.0,chi,-leem[s],chi,Ls-1,s);
}
// U_m^{-1} D^{-1}
for (int s=0;s<Ls-1;s++){
// Chi[s] + 1/d chi[s]
axpby_ssp_pplus(chi,1.0/dee[s],chi,-ueem[s]/dee[Ls-1],chi,s,Ls-1);
}
axpby_ssp(chi,1.0/dee[Ls-1],chi,0.0,chi,Ls-1,Ls-1); // Modest avoidable
// Apply U^{-1}
for (int s=Ls-2;s>=0;s--){
axpby_ssp_pminus (chi,1.0,chi,-uee[s],chi,s,s+1); // chi[Ls]
}
}
void CayleyFermion5D::MooeeInvDag (const LatticeFermion &psi, LatticeFermion &chi)
{
// Apply (U^{\prime})^{-dagger}
axpby_ssp (chi,1.0,psi, 0.0,psi,0,0); // chi[0]=psi[0]
for (int s=1;s<Ls;s++){
axpby_ssp_pminus(chi,1.0,psi,-uee[s-1],chi,s,s-1);
}
// U_m^{-\dagger}
for (int s=0;s<Ls-1;s++){
axpby_ssp_pplus(chi,1.0,chi,-ueem[s],chi,Ls-1,s);
}
// L_m^{-\dagger} D^{-dagger}
for (int s=0;s<Ls-1;s++){
axpby_ssp_pminus(chi,1.0/dee[s],chi,-leem[s]/dee[Ls-1],chi,s,Ls-1);
}
axpby_ssp(chi,1.0/dee[Ls-1],chi,0.0,chi,Ls-1,Ls-1); // Modest avoidable
// Apply L^{-dagger}
for (int s=Ls-2;s>=0;s--){
axpby_ssp_pplus (chi,1.0,chi,-lee[s],chi,s,s+1); // chi[Ls]
}
}
// Tanh
void CayleyFermion5D::SetCoefficientsTanh(Approx::zolotarev_data *zdata,RealD b,RealD c)
{
SetCoefficientsZolotarev(1.0,zdata,b,c);
}
//Zolo
void CayleyFermion5D::SetCoefficientsZolotarev(RealD zolo_hi,Approx::zolotarev_data *zdata,RealD b,RealD c)
{
///////////////////////////////////////////////////////////
// The Cayley coeffs (unprec)
///////////////////////////////////////////////////////////
omega.resize(Ls);
bs.resize(Ls);
cs.resize(Ls);
as.resize(Ls);
//
// Ts = ( [bs+cs]Dw )^-1 ( (bs+cs) Dw )
// -(g5 ------- -1 ) ( g5 --------- + 1 )
// ( {2+(bs-cs)Dw} ) ( 2+(bs-cs) Dw )
//
// bs = 1/2( (1/omega_s + 1)*b + (1/omega - 1)*c ) = 1/2( 1/omega(b+c) + (b-c) )
// cs = 1/2( (1/omega_s - 1)*b + (1/omega + 1)*c ) = 1/2( 1/omega(b+c) - (b-c) )
//
// bs+cs = 0.5*( 1/omega(b+c) + (b-c) + 1/omega(b+c) - (b-c) ) = 1/omega(b+c)
// bs-cs = 0.5*( 1/omega(b+c) + (b-c) - 1/omega(b+c) + (b-c) ) = b-c
//
// So
//
// Ts = ( [b+c]Dw/omega_s )^-1 ( (b+c) Dw /omega_s )
// -(g5 ------- -1 ) ( g5 --------- + 1 )
// ( {2+(b-c)Dw} ) ( 2+(b-c) Dw )
//
// Ts = ( [b+c]Dw )^-1 ( (b+c) Dw )
// -(g5 ------- -omega_s) ( g5 --------- + omega_s )
// ( {2+(b-c)Dw} ) ( 2+(b-c) Dw )
//
double bpc = b+c;
double bmc = b-c;
for(int i=0; i < Ls; i++){
as[i] = 1.0;
omega[i] = ((double)zdata->gamma[i])*zolo_hi; //NB reciprocal relative to Chroma NEF code
bs[i] = 0.5*(bpc/omega[i] + bmc);
cs[i] = 0.5*(bpc/omega[i] - bmc);
}
////////////////////////////////////////////////////////
// Constants for the preconditioned matrix Cayley form
////////////////////////////////////////////////////////
bee.resize(Ls);
cee.resize(Ls);
beo.resize(Ls);
ceo.resize(Ls);
for(int i=0;i<Ls;i++){
bee[i]=as[i]*(bs[i]*(4.0-M5) +1.0);
cee[i]=as[i]*(1.0-cs[i]*(4.0-M5));
beo[i]=as[i]*bs[i];
ceo[i]=-as[i]*cs[i];
}
aee.resize(Ls);
aeo.resize(Ls);
for(int i=0;i<Ls;i++){
aee[i]=cee[i];
aeo[i]=ceo[i];
}
//////////////////////////////////////////
// LDU decomposition of eeoo
//////////////////////////////////////////
dee.resize(Ls);
lee.resize(Ls);
leem.resize(Ls);
uee.resize(Ls);
ueem.resize(Ls);
for(int i=0;i<Ls;i++){
dee[i] = bee[i];
if ( i < Ls-1 ) {
lee[i] =-cee[i+1]/bee[i]; // sub-diag entry on the ith column
leem[i]=mass*cee[Ls-1]/bee[0];
for(int j=0;j<i;j++) leem[i]*= aee[j]/bee[j+1];
uee[i] =-aee[i]/bee[i]; // up-diag entry on the ith row
ueem[i]=mass;
for(int j=1;j<=i;j++) ueem[i]*= cee[j]/bee[j];
ueem[i]*= aee[0]/bee[0];
} else {
lee[i] =0.0;
leem[i]=0.0;
uee[i] =0.0;
ueem[i]=0.0;
}
}
{
double delta_d=mass*cee[Ls-1];
for(int j=0;j<Ls-1;j++) delta_d *= cee[j]/bee[j];
dee[Ls-1] += delta_d;
}
}
}}

62
lib/qcd/action/fermion/CayleyFermion5D.h Normal file
View File

@ -0,0 +1,62 @@
#ifndef GRID_QCD_CAYLEY_FERMION_H
#define GRID_QCD_CAYLEY_FERMION_H
namespace Grid {
namespace QCD {
class CayleyFermion5D : public WilsonFermion5D
{
public:
// override multiply
virtual RealD M (const LatticeFermion &in, LatticeFermion &out);
virtual RealD Mdag (const LatticeFermion &in, LatticeFermion &out);
// half checkerboard operations
virtual void Meooe (const LatticeFermion &in, LatticeFermion &out);
virtual void MeooeDag (const LatticeFermion &in, LatticeFermion &out);
virtual void Mooee (const LatticeFermion &in, LatticeFermion &out);
virtual void MooeeDag (const LatticeFermion &in, LatticeFermion &out);
virtual void MooeeInv (const LatticeFermion &in, LatticeFermion &out);
virtual void MooeeInvDag (const LatticeFermion &in, LatticeFermion &out);
virtual void Instantiatable(void)=0;
// protected:
RealD mass;
// Cayley form Moebius (tanh and zolotarev)
std::vector<RealD> omega;
std::vector<RealD> bs; // S dependent coeffs
std::vector<RealD> cs;
std::vector<RealD> as;
// For preconditioning Cayley form
std::vector<RealD> bee;
std::vector<RealD> cee;
std::vector<RealD> aee;
std::vector<RealD> beo;
std::vector<RealD> ceo;
std::vector<RealD> aeo;
// LDU factorisation of the eeoo matrix
std::vector<RealD> lee;
std::vector<RealD> leem;
std::vector<RealD> uee;
std::vector<RealD> ueem;
std::vector<RealD> dee;
// Constructors
CayleyFermion5D(LatticeGaugeField &_Umu,
GridCartesian &FiveDimGrid,
GridRedBlackCartesian &FiveDimRedBlackGrid,
GridCartesian &FourDimGrid,
GridRedBlackCartesian &FourDimRedBlackGrid,
RealD _mass,RealD _M5);
protected:
void SetCoefficientsZolotarev(RealD zolohi,Approx::zolotarev_data *zdata,RealD b,RealD c);
void SetCoefficientsTanh(Approx::zolotarev_data *zdata,RealD b,RealD c);
};
}
}
#endif

179
lib/qcd/action/fermion/ContinuedFractionFermion5D.cc Normal file
View File

@ -0,0 +1,179 @@
#include <Grid.h>
namespace Grid {
namespace QCD {
void ContinuedFractionFermion5D::SetCoefficientsTanh(Approx::zolotarev_data *zdata,RealD scale)
{
SetCoefficientsZolotarev(1.0/scale,zdata);
}
void ContinuedFractionFermion5D::SetCoefficientsZolotarev(RealD zolo_hi,Approx::zolotarev_data *zdata)
{
// How to check Ls matches??
// std::cout << Ls << " Ls"<<std::endl;
// std::cout << zdata->n << " - n"<<std::endl;
// std::cout << zdata->da << " -da "<<std::endl;
// std::cout << zdata->db << " -db"<<std::endl;
// std::cout << zdata->dn << " -dn"<<std::endl;
// std::cout << zdata->dd << " -dd"<<std::endl;
assert(zdata->db==Ls);// Beta has Ls coeffs
R=(1+this->mass)/(1-this->mass);
Beta.resize(Ls);
cc.resize(Ls);
cc_d.resize(Ls);
sqrt_cc.resize(Ls);
for(int i=0; i < Ls ; i++){
Beta[i] = zdata -> beta[i];
cc[i] = 1.0/Beta[i];
cc_d[i]=sqrt(cc[i]);
}
cc_d[Ls-1]=1.0;
for(int i=0; i < Ls-1 ; i++){
sqrt_cc[i]= sqrt(cc[i]*cc[i+1]);
}
sqrt_cc[Ls-2]=sqrt(cc[Ls-2]);
ZoloHiInv =1.0/zolo_hi;
dw_diag = (4.0-M5)*ZoloHiInv;
See.resize(Ls);
Aee.resize(Ls);
int sign=1;
for(int s=0;s<Ls;s++){
Aee[s] = sign * Beta[s] * dw_diag;
sign = - sign;
}
Aee[Ls-1] += R;
See[0] = Aee[0];
for(int s=1;s<Ls;s++){
See[s] = Aee[s] - 1.0/See[s-1];
}
for(int s=0;s<Ls;s++){
std::cout <<"s = "<<s<<" Beta "<<Beta[s]<<" Aee "<<Aee[s] <<" See "<<See[s] <<std::endl;
}
}
RealD ContinuedFractionFermion5D::M (const LatticeFermion &psi, LatticeFermion &chi)
{
LatticeFermion D(psi._grid);
DW(psi,D,DaggerNo);
int sign=1;
for(int s=0;s<Ls;s++){
if ( s==0 ) {
ag5xpby_ssp(chi,cc[0]*Beta[0]*sign*ZoloHiInv,D,sqrt_cc[0],psi,s,s+1); // Multiplies Dw by G5 so Hw
} else if ( s==(Ls-1) ){
RealD R=(1.0+mass)/(1.0-mass);
ag5xpby_ssp(chi,Beta[s]*ZoloHiInv,D,sqrt_cc[s-1],psi,s,s-1);
ag5xpby_ssp(chi,R,psi,1.0,chi,s,s);
} else {
ag5xpby_ssp(chi,cc[s]*Beta[s]*sign*ZoloHiInv,D,sqrt_cc[s],psi,s,s+1);
axpby_ssp(chi,1.0,chi,sqrt_cc[s-1],psi,s,s-1);
}
sign=-sign;
}
return norm2(chi);
}
RealD ContinuedFractionFermion5D::Mdag (const LatticeFermion &psi, LatticeFermion &chi)
{
// This matrix is already hermitian. (g5 Dw) = Dw dag g5 = (g5 Dw)dag
// The rest of matrix is symmetric.
// Can ignore "dag"
return M(psi,chi);
}
void ContinuedFractionFermion5D::Meooe (const LatticeFermion &psi, LatticeFermion &chi)
{
// Apply 4d dslash
if ( psi.checkerboard == Odd ) {
DhopEO(psi,chi,DaggerNo); // Dslash on diagonal. g5 Dslash is hermitian
} else {
DhopOE(psi,chi,DaggerNo); // Dslash on diagonal. g5 Dslash is hermitian
}
int sign=1;
for(int s=0;s<Ls;s++){
if ( s==(Ls-1) ){
ag5xpby_ssp(chi,Beta[s]*ZoloHiInv,chi,0.0,chi,s,s);
} else {
ag5xpby_ssp(chi,cc[s]*Beta[s]*sign*ZoloHiInv,chi,0.0,chi,s,s);
}
sign=-sign;
}
}
void ContinuedFractionFermion5D::MeooeDag (const LatticeFermion &psi, LatticeFermion &chi)
{
Meooe(psi,chi);
}
void ContinuedFractionFermion5D::Mooee (const LatticeFermion &psi, LatticeFermion &chi)
{
int sign=1;
for(int s=0;s<Ls;s++){
if ( s==0 ) {
ag5xpby_ssp(chi,cc[0]*Beta[0]*sign*dw_diag,psi,sqrt_cc[0],psi,s,s+1); // Multiplies Dw by G5 so Hw
} else if ( s==(Ls-1) ){
// Drop the CC here.
double R=(1+mass)/(1-mass);
ag5xpby_ssp(chi,Beta[s]*dw_diag,psi,sqrt_cc[s-1],psi,s,s-1);
ag5xpby_ssp(chi,R,psi,1.0,chi,s,s);
} else {
ag5xpby_ssp(chi,cc[s]*Beta[s]*sign*dw_diag,psi,sqrt_cc[s],psi,s,s+1);
axpby_ssp(chi,1.0,chi,sqrt_cc[s-1],psi,s,s-1);
}
sign=-sign;
}
}
void ContinuedFractionFermion5D::MooeeDag (const LatticeFermion &psi, LatticeFermion &chi)
{
Mooee(psi,chi);
}
void ContinuedFractionFermion5D::MooeeInv (const LatticeFermion &psi, LatticeFermion &chi)
{
// Apply Linv
axpby_ssp(chi,1.0/cc_d[0],psi,0.0,psi,0,0);
for(int s=1;s<Ls;s++){
axpbg5y_ssp(chi,1.0/cc_d[s],psi,-1.0/See[s-1],chi,s,s-1);
}
// Apply Dinv
for(int s=0;s<Ls;s++){
ag5xpby_ssp(chi,1.0/See[s],chi,0.0,chi,s,s); //only appearance of See[0]
}
// Apply Uinv = (Linv)^T
axpby_ssp(chi,1.0/cc_d[Ls-1],chi,0.0,chi,Ls-1,Ls-1);
for(int s=Ls-2;s>=0;s--){
axpbg5y_ssp(chi,1.0/cc_d[s],chi,-1.0*cc_d[s+1]/See[s]/cc_d[s],chi,s,s+1);
}
}
void ContinuedFractionFermion5D::MooeeInvDag (const LatticeFermion &psi, LatticeFermion &chi)
{
MooeeInv(psi,chi);
}
// Constructors
ContinuedFractionFermion5D::ContinuedFractionFermion5D(
LatticeGaugeField &_Umu,
GridCartesian &FiveDimGrid,
GridRedBlackCartesian &FiveDimRedBlackGrid,
GridCartesian &FourDimGrid,
GridRedBlackCartesian &FourDimRedBlackGrid,
RealD _mass,RealD M5) :
WilsonFermion5D(_Umu,
FiveDimGrid, FiveDimRedBlackGrid,
FourDimGrid, FourDimRedBlackGrid,M5),
mass(_mass)
{
assert((Ls&0x1)==1); // Odd Ls required
}
}
}

58
lib/qcd/action/fermion/ContinuedFractionFermion5D.h Normal file
View File

@ -0,0 +1,58 @@
#ifndef GRID_QCD_CONTINUED_FRACTION_H
#define GRID_QCD_CONTINUED_FRACTION_H
namespace Grid {
namespace QCD {
class ContinuedFractionFermion5D : public WilsonFermion5D
{
public:
// override multiply
virtual RealD M (const LatticeFermion &in, LatticeFermion &out);
virtual RealD Mdag (const LatticeFermion &in, LatticeFermion &out);
// half checkerboard operaions
virtual void Meooe (const LatticeFermion &in, LatticeFermion &out);
virtual void MeooeDag (const LatticeFermion &in, LatticeFermion &out);
virtual void Mooee (const LatticeFermion &in, LatticeFermion &out);
virtual void MooeeDag (const LatticeFermion &in, LatticeFermion &out);
virtual void MooeeInv (const LatticeFermion &in, LatticeFermion &out);
virtual void MooeeInvDag (const LatticeFermion &in, LatticeFermion &out);
// virtual void Instantiatable(void)=0;
virtual void Instantiatable(void) =0;
// Constructors
ContinuedFractionFermion5D(LatticeGaugeField &_Umu,
GridCartesian &FiveDimGrid,
GridRedBlackCartesian &FiveDimRedBlackGrid,
GridCartesian &FourDimGrid,
GridRedBlackCartesian &FourDimRedBlackGrid,
RealD _mass,RealD M5);
protected:
void SetCoefficientsTanh(Approx::zolotarev_data *zdata,RealD scale);
void SetCoefficientsZolotarev(RealD zolo_hi,Approx::zolotarev_data *zdata);;
// Cont frac
RealD dw_diag;
RealD mass;
RealD R;
RealD ZoloHiInv;
std::vector<double> Beta;
std::vector<double> cc;;
std::vector<double> cc_d;;
std::vector<double> sqrt_cc;
std::vector<double> See;
std::vector<double> Aee;
};
}
}
#endif

47
lib/qcd/action/fermion/DomainWallFermion.h Normal file
View File

@ -0,0 +1,47 @@
#ifndef GRID_QCD_DOMAIN_WALL_FERMION_H
#define GRID_QCD_DOMAIN_WALL_FERMION_H
#include <Grid.h>
namespace Grid {
namespace QCD {
class DomainWallFermion : public CayleyFermion5D
{
public:
virtual void Instantiatable(void) {};
// Constructors
DomainWallFermion(LatticeGaugeField &_Umu,
GridCartesian &FiveDimGrid,
GridRedBlackCartesian &FiveDimRedBlackGrid,
GridCartesian &FourDimGrid,
GridRedBlackCartesian &FourDimRedBlackGrid,
RealD _mass,RealD _M5) :
CayleyFermion5D(_Umu,
FiveDimGrid,
FiveDimRedBlackGrid,
FourDimGrid,
FourDimRedBlackGrid,_mass,_M5)
{
RealD eps = 1.0;
Approx::zolotarev_data *zdata = Approx::higham(eps,this->Ls);// eps is ignored for higham
assert(zdata->n==this->Ls);
std::cout << "DomainWallFermion with Ls="<<Ls<<std::endl;
// Call base setter
this->CayleyFermion5D::SetCoefficientsTanh(zdata,1.0,0.0);
Approx::zolotarev_free(zdata);
}
};
}
}
#endif

7
lib/qcd/action/fermion/FermionAction.h → lib/qcd/action/fermion/FermionOperator.h
View File

@ -1,5 +1,5 @@
#ifndef GRID_QCD_WILSON_DOP_H
#define GRID_QCD_WILSON_DOP_H
#ifndef GRID_QCD_FERMION_OPERATOR_H
#define GRID_QCD_FERMION_OPERATOR_H
namespace Grid {
@ -11,7 +11,7 @@ namespace Grid {
// Think about multiple representations
//////////////////////////////////////////////////////////////////////////////
template<class FermionField,class GaugeField>
class FermionAction : public CheckerBoardedSparseMatrixBase<FermionField>
class FermionOperator : public CheckerBoardedSparseMatrixBase<FermionField>
{
public:
@ -40,6 +40,7 @@ namespace Grid {
virtual void DhopOE(const FermionField &in, FermionField &out,int dag)=0;
virtual void DhopEO(const FermionField &in, FermionField &out,int dag)=0;
};
}

49
lib/qcd/action/fermion/MobiusFermion.h Normal file
View File

@ -0,0 +1,49 @@
#ifndef GRID_QCD_MOBIUS_FERMION_H
#define GRID_QCD_MOBIUS_FERMION_H
#include <Grid.h>
namespace Grid {
namespace QCD {
class MobiusFermion : public CayleyFermion5D
{
public:
virtual void Instantiatable(void) {};
// Constructors
MobiusFermion(LatticeGaugeField &_Umu,
GridCartesian &FiveDimGrid,
GridRedBlackCartesian &FiveDimRedBlackGrid,
GridCartesian &FourDimGrid,
GridRedBlackCartesian &FourDimRedBlackGrid,
RealD _mass,RealD _M5,
RealD b, RealD c) :
CayleyFermion5D(_Umu,
FiveDimGrid,
FiveDimRedBlackGrid,
FourDimGrid,
FourDimRedBlackGrid,_mass,_M5)
{
RealD eps = 1.0;
std::cout << "MobiusFermion (b="<<b<<",c="<<c<<") with Ls= "<<Ls<<" Tanh approx"<<std::endl;
Approx::zolotarev_data *zdata = Approx::higham(eps,this->Ls);// eps is ignored for higham
assert(zdata->n==this->Ls);
// Call base setter
this->CayleyFermion5D::SetCoefficientsTanh(zdata,b,c);
Approx::zolotarev_free(zdata);
}
};
}
}
#endif

50
lib/qcd/action/fermion/MobiusZolotarevFermion.h Normal file
View File

@ -0,0 +1,50 @@
#ifndef GRID_QCD_MOBIUS_ZOLOTAREV_FERMION_H
#define GRID_QCD_MOBIUS_ZOLOTAREV_FERMION_H
#include <Grid.h>
namespace Grid {
namespace QCD {
class MobiusZolotarevFermion : public CayleyFermion5D
{
public:
virtual void Instantiatable(void) {};
// Constructors
MobiusZolotarevFermion(LatticeGaugeField &_Umu,
GridCartesian &FiveDimGrid,
GridRedBlackCartesian &FiveDimRedBlackGrid,
GridCartesian &FourDimGrid,
GridRedBlackCartesian &FourDimRedBlackGrid,
RealD _mass,RealD _M5,
RealD b, RealD c,
RealD lo, RealD hi) :
CayleyFermion5D(_Umu,
FiveDimGrid,
FiveDimRedBlackGrid,
FourDimGrid,
FourDimRedBlackGrid,_mass,_M5)
{
RealD eps = lo/hi;
Approx::zolotarev_data *zdata = Approx::zolotarev(eps,this->Ls,0);
assert(zdata->n==this->Ls);
std::cout << "MobiusZolotarevFermion (b="<<b<<",c="<<c<<") with Ls= "<<Ls<<" Zolotarev range ["<<lo<<","<<hi<<"]"<<std::endl;
// Call base setter
this->CayleyFermion5D::SetCoefficientsZolotarev(hi,zdata,b,c);
Approx::zolotarev_free(zdata);
}
};
}
}
#endif

34
lib/qcd/action/fermion/OverlapWilsonCayleyTanhFermion.h Normal file
View File

@ -0,0 +1,34 @@
#ifndef OVERLAP_WILSON_CAYLEY_TANH_FERMION_H
#define OVERLAP_WILSON_CAYLEY_TANH_FERMION_H
#include <Grid.h>
namespace Grid {
namespace QCD {
class OverlapWilsonCayleyTanhFermion : public MobiusFermion
{
public:
// Constructors
OverlapWilsonCayleyTanhFermion(LatticeGaugeField &_Umu,
GridCartesian &FiveDimGrid,
GridRedBlackCartesian &FiveDimRedBlackGrid,
GridCartesian &FourDimGrid,
GridRedBlackCartesian &FourDimRedBlackGrid,
RealD _mass,RealD _M5,
RealD scale) :
// b+c=scale, b-c = 0 <=> b =c = scale/2
MobiusFermion(_Umu,
FiveDimGrid,
FiveDimRedBlackGrid,
FourDimGrid,
FourDimRedBlackGrid,_mass,_M5,0.5*scale,0.5*scale)
{
}
};
}
}
#endif

37
lib/qcd/action/fermion/OverlapWilsonCayleyZolotarevFermion.h Normal file
View File

@ -0,0 +1,37 @@
#ifndef OVERLAP_WILSON_CAYLEY_ZOLOTAREV_FERMION_H
#define OVERLAP_WILSON_CAYLEY_ZOLOTAREV_FERMION_H
#include <Grid.h>
namespace Grid {
namespace QCD {
class OverlapWilsonCayleyZolotarevFermion : public MobiusZolotarevFermion
{
public:
// Constructors
OverlapWilsonCayleyZolotarevFermion(LatticeGaugeField &_Umu,
GridCartesian &FiveDimGrid,
GridRedBlackCartesian &FiveDimRedBlackGrid,
GridCartesian &FourDimGrid,
GridRedBlackCartesian &FourDimRedBlackGrid,
RealD _mass,RealD _M5,
RealD lo, RealD hi) :
// b+c=1.0, b-c = 0 <=> b =c = 1/2
MobiusZolotarevFermion(_Umu,
FiveDimGrid,
FiveDimRedBlackGrid,
FourDimGrid,
FourDimRedBlackGrid,_mass,_M5,0.5,0.5,lo,hi)
{}
};
}
}
#endif

40
lib/qcd/action/fermion/OverlapWilsonContfracTanhFermion.h Normal file
View File

@ -0,0 +1,40 @@
#ifndef OVERLAP_WILSON_CONTFRAC_TANH_FERMION_H
#define OVERLAP_WILSON_CONTFRAC_TANH_FERMION_H
#include <Grid.h>
namespace Grid {
namespace QCD {
class OverlapWilsonContFracTanhFermion : public ContinuedFractionFermion5D
{
public:
virtual void Instantiatable(void){};
// Constructors
OverlapWilsonContFracTanhFermion(LatticeGaugeField &_Umu,
GridCartesian &FiveDimGrid,
GridRedBlackCartesian &FiveDimRedBlackGrid,
GridCartesian &FourDimGrid,
GridRedBlackCartesian &FourDimRedBlackGrid,
RealD _mass,RealD _M5,
RealD scale) :
// b+c=scale, b-c = 0 <=> b =c = scale/2
ContinuedFractionFermion5D(_Umu,
FiveDimGrid,
FiveDimRedBlackGrid,
FourDimGrid,
FourDimRedBlackGrid,_mass,_M5)
{
assert((Ls&0x1)==1); // Odd Ls required
int nrational=Ls-1;// Even rational order
Approx::zolotarev_data *zdata = Approx::higham(1.0,nrational);// eps is ignored for higham
SetCoefficientsTanh(zdata,scale);
Approx::zolotarev_free(zdata);
}
};
}
}
#endif

44
lib/qcd/action/fermion/OverlapWilsonContfracZolotarevFermion.h Normal file
View File

@ -0,0 +1,44 @@
#ifndef OVERLAP_WILSON_CONTFRAC_ZOLOTAREV_FERMION_H
#define OVERLAP_WILSON_CONTFRAC_ZOLOTAREV_FERMION_H
#include <Grid.h>
namespace Grid {
namespace QCD {
class OverlapWilsonContFracZolotarevFermion : public ContinuedFractionFermion5D
{
public:
virtual void Instantiatable(void){};
// Constructors
OverlapWilsonContFracZolotarevFermion(LatticeGaugeField &_Umu,
GridCartesian &FiveDimGrid,
GridRedBlackCartesian &FiveDimRedBlackGrid,
GridCartesian &FourDimGrid,
GridRedBlackCartesian &FourDimRedBlackGrid,
RealD _mass,RealD _M5,
RealD lo,RealD hi):
// b+c=scale, b-c = 0 <=> b =c = scale/2
ContinuedFractionFermion5D(_Umu,
FiveDimGrid,
FiveDimRedBlackGrid,
FourDimGrid,
FourDimRedBlackGrid,_mass,_M5)
{
assert((Ls&0x1)==1); // Odd Ls required
int nrational=Ls;// Odd rational order
RealD eps = lo/hi;
Approx::zolotarev_data *zdata = Approx::zolotarev(eps,nrational,0);
SetCoefficientsZolotarev(hi,zdata);
Approx::zolotarev_free(zdata);
}
};
}
}
#endif

40
lib/qcd/action/fermion/OverlapWilsonPartialFractionTanhFermion.h Normal file
View File

@ -0,0 +1,40 @@
#ifndef OVERLAP_WILSON_PARTFRAC_TANH_FERMION_H
#define OVERLAP_WILSON_PARTFRAC_TANH_FERMION_H
#include <Grid.h>
namespace Grid {
namespace QCD {
class OverlapWilsonPartialFractionTanhFermion : public PartialFractionFermion5D
{
public:
virtual void Instantiatable(void){};
// Constructors
OverlapWilsonPartialFractionTanhFermion(LatticeGaugeField &_Umu,
GridCartesian &FiveDimGrid,
GridRedBlackCartesian &FiveDimRedBlackGrid,
GridCartesian &FourDimGrid,
GridRedBlackCartesian &FourDimRedBlackGrid,
RealD _mass,RealD _M5,
RealD scale) :
// b+c=scale, b-c = 0 <=> b =c = scale/2
PartialFractionFermion5D(_Umu,
FiveDimGrid,
FiveDimRedBlackGrid,
FourDimGrid,
FourDimRedBlackGrid,_mass,_M5)
{
assert((Ls&0x1)==1); // Odd Ls required
int nrational=Ls-1;// Even rational order
Approx::zolotarev_data *zdata = Approx::higham(1.0,nrational);// eps is ignored for higham
SetCoefficientsTanh(zdata,scale);
Approx::zolotarev_free(zdata);
}
};
}
}
#endif

44
lib/qcd/action/fermion/OverlapWilsonPartialFractionZolotarevFermion.h Normal file
View File

@ -0,0 +1,44 @@
#ifndef OVERLAP_WILSON_PARTFRAC_ZOLOTAREV_FERMION_H
#define OVERLAP_WILSON_PARTFRAC_ZOLOTAREV_FERMION_H
#include <Grid.h>
namespace Grid {
namespace QCD {
class OverlapWilsonPartialFractionZolotarevFermion : public PartialFractionFermion5D
{
public:
virtual void Instantiatable(void){};
// Constructors
OverlapWilsonPartialFractionZolotarevFermion(LatticeGaugeField &_Umu,
GridCartesian &FiveDimGrid,
GridRedBlackCartesian &FiveDimRedBlackGrid,
GridCartesian &FourDimGrid,
GridRedBlackCartesian &FourDimRedBlackGrid,
RealD _mass,RealD _M5,
RealD lo,RealD hi):
// b+c=scale, b-c = 0 <=> b =c = scale/2
PartialFractionFermion5D(_Umu,
FiveDimGrid,
FiveDimRedBlackGrid,
FourDimGrid,
FourDimRedBlackGrid,_mass,_M5)
{
assert((Ls&0x1)==1); // Odd Ls required
int nrational=Ls;// Odd rational order
RealD eps = lo/hi;
Approx::zolotarev_data *zdata = Approx::zolotarev(eps,nrational,0);
SetCoefficientsZolotarev(hi,zdata);
Approx::zolotarev_free(zdata);
}
};
}
}
#endif

310
lib/qcd/action/fermion/PartialFractionFermion5D.cc Normal file
View File

@ -0,0 +1,310 @@
#include <Grid.h>
namespace Grid {
namespace QCD {
void PartialFractionFermion5D::Meooe_internal(const LatticeFermion &psi, LatticeFermion &chi,int dag)
{
// this does both dag and undag but is trivial; make a common helper routing
int sign = dag ? (-1) : 1;
if ( psi.checkerboard == Odd ) {
DhopEO(psi,chi,DaggerNo);
} else {
DhopOE(psi,chi,DaggerNo);
}
int nblock=(Ls-1)/2;
for(int b=0;b<nblock;b++){
int s = 2*b;
ag5xpby_ssp(chi,-scale,chi,0.0,chi,s,s);
ag5xpby_ssp(chi, scale,chi,0.0,chi,s+1,s+1);
}
ag5xpby_ssp(chi,p[nblock]*scale/amax,chi,0.0,chi,Ls-1,Ls-1);
}
void PartialFractionFermion5D::Mooee_internal(const LatticeFermion &psi, LatticeFermion &chi,int dag)
{
// again dag and undag are trivially related
int sign = dag ? (-1) : 1;
int nblock=(Ls-1)/2;
for(int b=0;b<nblock;b++){
int s = 2*b;
RealD pp = p[nblock-1-b];
RealD qq = q[nblock-1-b];
// Do each 2x2 block aligned at s and multiplies Dw site diagonal by G5 so Hw
ag5xpby_ssp(chi,-dw_diag*scale,psi,amax*sqrt(qq)*scale,psi, s ,s+1);
ag5xpby_ssp(chi, dw_diag*scale,psi,amax*sqrt(qq)*scale,psi, s+1,s);
axpby_ssp (chi, 1.0, chi,sqrt(amax*pp)*scale*sign,psi,s+1,Ls-1);
}
{
RealD R=(1+mass)/(1-mass);
//R g5 psi[Ls-1] + p[0] H
ag5xpbg5y_ssp(chi,R*scale,psi,p[nblock]*scale*dw_diag/amax,psi,Ls-1,Ls-1);
for(int b=0;b<nblock;b++){
int s = 2*b+1;
RealD pp = p[nblock-1-b];
axpby_ssp(chi,1.0,chi,-sqrt(amax*pp)*scale*sign,psi,Ls-1,s);
}
}
}
void PartialFractionFermion5D::MooeeInv_internal(const LatticeFermion &psi, LatticeFermion &chi,int dag)
{
int sign = dag ? (-1) : 1;
LatticeFermion tmp(psi._grid);
///////////////////////////////////////////////////////////////////////////////////////
//Linv
///////////////////////////////////////////////////////////////////////////////////////
int nblock=(Ls-1)/2;
axpy(chi,0.0,psi,psi); // Identity piece
for(int b=0;b<nblock;b++){
int s = 2*b;
RealD pp = p[nblock-1-b];
RealD qq = q[nblock-1-b];
RealD coeff1=sign*sqrt(amax*amax*amax*pp*qq) / ( dw_diag*dw_diag + amax*amax* qq);
RealD coeff2=sign*sqrt(amax*pp)*dw_diag / ( dw_diag*dw_diag + amax*amax* qq); // Implicit g5 here
axpby_ssp (chi,1.0,chi,coeff1,psi,Ls-1,s);
axpbg5y_ssp(chi,1.0,chi,coeff2,psi,Ls-1,s+1);
}
///////////////////////////////////////////////////////////////////////////////////////
//Dinv (note D isn't really diagonal -- just diagonal enough that we can still invert)
// Compute Seeinv (coeff of gamma5)
///////////////////////////////////////////////////////////////////////////////////////
RealD R=(1+mass)/(1-mass);
RealD Seeinv = R + p[nblock]*dw_diag/amax;
for(int b=0;b<nblock;b++){
Seeinv += p[nblock-1-b]*dw_diag/amax / ( dw_diag*dw_diag/amax/amax + q[nblock-1-b]);
}
Seeinv = 1.0/Seeinv;
for(int b=0;b<nblock;b++){
int s = 2*b;
RealD pp = p[nblock-1-b];
RealD qq = q[nblock-1-b];
RealD coeff1=dw_diag / ( dw_diag*dw_diag + amax*amax* qq); // Implicit g5 here
RealD coeff2=amax*sqrt(qq) / ( dw_diag*dw_diag + amax*amax* qq);
ag5xpby_ssp (tmp,-coeff1,chi,coeff2,chi,s,s+1);
ag5xpby_ssp (tmp, coeff1,chi,coeff2,chi,s+1,s);
}
ag5xpby_ssp (tmp, Seeinv,chi,0.0,chi,Ls-1,Ls-1);
///////////////////////////////////////////////////////////////////////////////////////
// Uinv
///////////////////////////////////////////////////////////////////////////////////////
for(int b=0;b<nblock;b++){
int s = 2*b;
RealD pp = p[nblock-1-b];
RealD qq = q[nblock-1-b];
RealD coeff1=-sign*sqrt(amax*amax*amax*pp*qq) / ( dw_diag*dw_diag + amax*amax* qq);
RealD coeff2=-sign*sqrt(amax*pp)*dw_diag / ( dw_diag*dw_diag + amax*amax* qq); // Implicit g5 here
axpby_ssp (chi,1.0/scale,tmp,coeff1/scale,tmp,s,Ls-1);
axpbg5y_ssp(chi,1.0/scale,tmp,coeff2/scale,tmp,s+1,Ls-1);
}
axpby_ssp (chi, 1.0/scale,tmp,0.0,tmp,Ls-1,Ls-1);
}
void PartialFractionFermion5D::M_internal(const LatticeFermion &psi, LatticeFermion &chi,int dag)
{
LatticeFermion D(psi._grid);
int sign = dag ? (-1) : 1;
// For partial frac Hw case (b5=c5=1) chroma quirkily computes
//
// Conventions for partfrac appear to be a mess.
// Tony's Nara lectures have
//
// BlockDiag( H/p_i 1 | 1 )
// ( 1 p_i H / q_i^2 | 0 )
// ---------------------------------
// ( -1 0 | R +p0 H )
//
//Chroma ( -2H 2sqrt(q_i) | 0 )
// (2 sqrt(q_i) 2H | 2 sqrt(p_i) )
// ---------------------------------
// ( 0 -2 sqrt(p_i) | 2 R gamma_5 + p0 2H
//
// Edwards/Joo/Kennedy/Wenger
//
// Here, the "beta's" selected by chroma to scale the unphysical bulk constraint fields
// incorporate the approx scale factor. This is obtained by propagating the
// scale on "H" out to the off diagonal elements as follows:
//
// BlockDiag( H/p_i 1 | 1 )
// ( 1 p_i H / q_i^2 | 0 )
// ---------------------------------
// ( -1 0 | R + p_0 H )
//
// becomes:
// BlockDiag( H/ sp_i 1 | 1 )
// ( 1 sp_i H / s^2q_i^2 | 0 )
// ---------------------------------
// ( -1 0 | R + p_0/s H )
//
//
// This is implemented in Chroma by
// p0' = p0/approxMax
// p_i' = p_i*approxMax
// q_i' = q_i*approxMax*approxMax
//
// After the equivalence transform is applied the matrix becomes
//
//Chroma ( -2H sqrt(q'_i) | 0 )
// (sqrt(q'_i) 2H | sqrt(p'_i) )
// ---------------------------------
// ( 0 -sqrt(p'_i) | 2 R gamma_5 + p'0 2H
//
// = ( -2H sqrt(q_i)amax | 0 )
// (sqrt(q_i)amax 2H | sqrt(p_i*amax) )
// ---------------------------------
// ( 0 -sqrt(p_i)*amax | 2 R gamma_5 + p0/amax 2H
//
DW(psi,D,DaggerNo);
int nblock=(Ls-1)/2;
for(int b=0;b<nblock;b++){
int s = 2*b;
double pp = p[nblock-1-b];
double qq = q[nblock-1-b];
// Do each 2x2 block aligned at s and
ag5xpby_ssp(chi,-1.0*scale,D,amax*sqrt(qq)*scale,psi, s ,s+1); // Multiplies Dw by G5 so Hw
ag5xpby_ssp(chi, 1.0*scale,D,amax*sqrt(qq)*scale,psi, s+1,s);
// Pick up last column
axpby_ssp (chi, 1.0, chi,sqrt(amax*pp)*scale*sign,psi,s+1,Ls-1);
}
{
double R=(1+this->mass)/(1-this->mass);
//R g5 psi[Ls] + p[0] H
ag5xpbg5y_ssp(chi,R*scale,psi,p[nblock]*scale/amax,D,Ls-1,Ls-1);
for(int b=0;b<nblock;b++){
int s = 2*b+1;
double pp = p[nblock-1-b];
axpby_ssp(chi,1.0,chi,-sqrt(amax*pp)*scale*sign,psi,Ls-1,s);
}
}
}
RealD PartialFractionFermion5D::M (const LatticeFermion &in, LatticeFermion &out)
{
M_internal(in,out,DaggerNo);
return norm2(out);
}
RealD PartialFractionFermion5D::Mdag (const LatticeFermion &in, LatticeFermion &out)
{
M_internal(in,out,DaggerYes);
return norm2(out);
}
void PartialFractionFermion5D::Meooe (const LatticeFermion &in, LatticeFermion &out)
{
Meooe_internal(in,out,DaggerNo);
}
void PartialFractionFermion5D::MeooeDag (const LatticeFermion &in, LatticeFermion &out)
{
Meooe_internal(in,out,DaggerYes);
}
void PartialFractionFermion5D::Mooee (const LatticeFermion &in, LatticeFermion &out)
{
Mooee_internal(in,out,DaggerNo);
}
void PartialFractionFermion5D::MooeeDag (const LatticeFermion &in, LatticeFermion &out)
{
Mooee_internal(in,out,DaggerYes);
}
void PartialFractionFermion5D::MooeeInv (const LatticeFermion &in, LatticeFermion &out)
{
MooeeInv_internal(in,out,DaggerNo);
}
void PartialFractionFermion5D::MooeeInvDag (const LatticeFermion &in, LatticeFermion &out)
{
MooeeInv_internal(in,out,DaggerYes);
}
void PartialFractionFermion5D::SetCoefficientsTanh(Approx::zolotarev_data *zdata,RealD scale){
SetCoefficientsZolotarev(1.0/scale,zdata);
}
void PartialFractionFermion5D::SetCoefficientsZolotarev(RealD zolo_hi,Approx::zolotarev_data *zdata){
// check on degree matching
// std::cout << Ls << " Ls"<<std::endl;
// std::cout << zdata->n << " - n"<<std::endl;
// std::cout << zdata->da << " -da "<<std::endl;
// std::cout << zdata->db << " -db"<<std::endl;
// std::cout << zdata->dn << " -dn"<<std::endl;
// std::cout << zdata->dd << " -dd"<<std::endl;
assert(Ls == (2*zdata->da -1) );
// Part frac
// RealD R;
R=(1+mass)/(1-mass);
dw_diag = (4.0-M5);
// std::vector<RealD> p;
// std::vector<RealD> q;
p.resize(zdata->da);
q.resize(zdata->dd);
for(int n=0;n<zdata->da;n++){
p[n] = zdata -> alpha[n];
}
for(int n=0;n<zdata->dd;n++){
q[n] = -zdata -> ap[n];
}
scale= part_frac_chroma_convention ? 2.0 : 1.0; // Chroma conventions annoy me
amax=zolo_hi;
}
// Constructors
PartialFractionFermion5D::PartialFractionFermion5D(LatticeGaugeField &_Umu,
GridCartesian &FiveDimGrid,
GridRedBlackCartesian &FiveDimRedBlackGrid,
GridCartesian &FourDimGrid,
GridRedBlackCartesian &FourDimRedBlackGrid,
RealD _mass,RealD M5) :
WilsonFermion5D(_Umu,
FiveDimGrid, FiveDimRedBlackGrid,
FourDimGrid, FourDimRedBlackGrid,M5),
mass(_mass)
{
assert((Ls&0x1)==1); // Odd Ls required
int nrational=Ls-1;
Approx::zolotarev_data *zdata = Approx::higham(1.0,nrational);
// NB: chroma uses a cast to "float" for the zolotarev range(!?).
// this creates a real difference in the operator which I do not like but we can replicate here
// to demonstrate compatibility
// RealD eps = (zolo_lo / zolo_hi);
// zdata = bfm_zolotarev(eps,nrational,0);
SetCoefficientsTanh(zdata,1.0);
Approx::zolotarev_free(zdata);
}
}
}

61
lib/qcd/action/fermion/PartialFractionFermion5D.h Normal file
View File

@ -0,0 +1,61 @@
#ifndef GRID_QCD_PARTIAL_FRACTION_H
#define GRID_QCD_PARTIAL_FRACTION_H
namespace Grid {
namespace QCD {
class PartialFractionFermion5D : public WilsonFermion5D
{
public:
const int part_frac_chroma_convention=1;
void Meooe_internal(const LatticeFermion &in, LatticeFermion &out,int dag);
void Mooee_internal(const LatticeFermion &in, LatticeFermion &out,int dag);
void MooeeInv_internal(const LatticeFermion &in, LatticeFermion &out,int dag);
void M_internal(const LatticeFermion &in, LatticeFermion &out,int dag);
// override multiply
virtual RealD M (const LatticeFermion &in, LatticeFermion &out);
virtual RealD Mdag (const LatticeFermion &in, LatticeFermion &out);
// half checkerboard operaions
virtual void Meooe (const LatticeFermion &in, LatticeFermion &out);
virtual void MeooeDag (const LatticeFermion &in, LatticeFermion &out);
virtual void Mooee (const LatticeFermion &in, LatticeFermion &out);
virtual void MooeeDag (const LatticeFermion &in, LatticeFermion &out);
virtual void MooeeInv (const LatticeFermion &in, LatticeFermion &out);
virtual void MooeeInvDag (const LatticeFermion &in, LatticeFermion &out);
virtual void Instantiatable(void) =0; // ensure no make-eee
// Constructors
PartialFractionFermion5D(LatticeGaugeField &_Umu,
GridCartesian &FiveDimGrid,
GridRedBlackCartesian &FiveDimRedBlackGrid,
GridCartesian &FourDimGrid,
GridRedBlackCartesian &FourDimRedBlackGrid,
RealD _mass,RealD M5);
protected:
virtual void SetCoefficientsTanh(Approx::zolotarev_data *zdata,RealD scale);
virtual void SetCoefficientsZolotarev(RealD zolo_hi,Approx::zolotarev_data *zdata);
// Part frac
RealD mass;
RealD dw_diag;
RealD R;
RealD amax;
RealD scale;
std::vector<double> p;
std::vector<double> q;
};
}
}
#endif

37
lib/qcd/action/fermion/ScaledShamirFermion.h Normal file
View File

@ -0,0 +1,37 @@
#ifndef GRID_QCD_SCALED_SHAMIR_FERMION_H
#define GRID_QCD_SCALED_SHAMIR_FERMION_H
#include <Grid.h>
namespace Grid {
namespace QCD {
class ScaledShamirFermion : public MobiusFermion
{
public:
// Constructors
ScaledShamirFermion(LatticeGaugeField &_Umu,
GridCartesian &FiveDimGrid,
GridRedBlackCartesian &FiveDimRedBlackGrid,
GridCartesian &FourDimGrid,
GridRedBlackCartesian &FourDimRedBlackGrid,
RealD _mass,RealD _M5,
RealD scale) :
// b+c=scale, b-c = 1 <=> 2b = scale+1; 2c = scale-1
MobiusFermion(_Umu,
FiveDimGrid,
FiveDimRedBlackGrid,
FourDimGrid,
FourDimRedBlackGrid,_mass,_M5,0.5*(scale+1.0),0.5*(scale-1.0))
{
}
};
}
}
#endif

39
lib/qcd/action/fermion/ShamirZolotarevFermion.h Normal file
View File

@ -0,0 +1,39 @@
#ifndef GRID_QCD_SHAMIR_ZOLOTAREV_FERMION_H
#define GRID_QCD_SHAMIR_ZOLOTAREV_FERMION_H
#include <Grid.h>
namespace Grid {
namespace QCD {
class ShamirZolotarevFermion : public MobiusZolotarevFermion
{
public:
// Constructors
ShamirZolotarevFermion(LatticeGaugeField &_Umu,
GridCartesian &FiveDimGrid,
GridRedBlackCartesian &FiveDimRedBlackGrid,
GridCartesian &FourDimGrid,
GridRedBlackCartesian &FourDimRedBlackGrid,
RealD _mass,RealD _M5,
RealD lo, RealD hi) :
// b+c = 1; b-c = 1 => b=1, c=0
MobiusZolotarevFermion(_Umu,
FiveDimGrid,
FiveDimRedBlackGrid,
FourDimGrid,
FourDimRedBlackGrid,_mass,_M5,1.0,0.0,lo,hi)
{}
};
}
}
#endif

6
lib/qcd/action/fermion/WilsonFermion.cc
View File

@ -9,9 +9,9 @@ const std::vector<int> WilsonFermion::displacements({1,1,1,1,-1,-1,-1,-1});
int WilsonFermion::HandOptDslash;
WilsonFermion::WilsonFermion(LatticeGaugeField &_Umu,
GridCartesian &Fgrid,
GridRedBlackCartesian &Hgrid,
double _mass) :
GridCartesian &Fgrid,
GridRedBlackCartesian &Hgrid,
RealD _mass) :
_grid(&Fgrid),
_cbgrid(&Hgrid),
Stencil (&Fgrid,npoint,Even,directions,displacements),

6
lib/qcd/action/fermion/WilsonFermion.h
View File

@ -5,7 +5,7 @@ namespace Grid {
namespace QCD {
class WilsonFermion : public FermionAction<LatticeFermion,LatticeGaugeField>
class WilsonFermion : public FermionOperator<LatticeFermion,LatticeGaugeField>
{
public:
@ -44,7 +44,7 @@ namespace Grid {
int dag);
// Constructor
WilsonFermion(LatticeGaugeField &_Umu,GridCartesian &Fgrid,GridRedBlackCartesian &Hgrid,double _mass);
WilsonFermion(LatticeGaugeField &_Umu,GridCartesian &Fgrid,GridRedBlackCartesian &Hgrid,RealD _mass);
// DoubleStore
void DoubleStore(LatticeDoubledGaugeField &Uds,const LatticeGaugeField &Umu);
@ -57,7 +57,7 @@ namespace Grid {
protected:
double mass;
RealD mass;
GridBase * _grid;
GridBase * _cbgrid;

106
lib/qcd/action/fermion/FiveDimWilsonFermion.cc → lib/qcd/action/fermion/WilsonFermion5D.cc
View File

@ -4,18 +4,18 @@ namespace Grid {
namespace QCD {
// S-direction is INNERMOST and takes no part in the parity.
const std::vector<int> FiveDimWilsonFermion::directions ({1,2,3,4, 1, 2, 3, 4});
const std::vector<int> FiveDimWilsonFermion::displacements({1,1,1,1,-1,-1,-1,-1});
const std::vector<int> WilsonFermion5D::directions ({1,2,3,4, 1, 2, 3, 4});
const std::vector<int> WilsonFermion5D::displacements({1,1,1,1,-1,-1,-1,-1});
int FiveDimWilsonFermion::HandOptDslash;
int WilsonFermion5D::HandOptDslash;
// 5d lattice for DWF.
FiveDimWilsonFermion::FiveDimWilsonFermion(LatticeGaugeField &_Umu,
WilsonFermion5D::WilsonFermion5D(LatticeGaugeField &_Umu,
GridCartesian &FiveDimGrid,
GridRedBlackCartesian &FiveDimRedBlackGrid,
GridCartesian &FourDimGrid,
GridRedBlackCartesian &FourDimRedBlackGrid,
double _mass) :
RealD _M5) :
_FiveDimGrid(&FiveDimGrid),
_FiveDimRedBlackGrid(&FiveDimRedBlackGrid),
_FourDimGrid(&FourDimGrid),
@ -23,7 +23,7 @@ namespace QCD {
Stencil (_FiveDimGrid,npoint,Even,directions,displacements),
StencilEven(_FiveDimRedBlackGrid,npoint,Even,directions,displacements), // source is Even
StencilOdd (_FiveDimRedBlackGrid,npoint,Odd ,directions,displacements), // source is Odd
mass(_mass),
M5(_M5),
Umu(_FourDimGrid),
UmuEven(_FourDimRedBlackGrid),
UmuOdd (_FourDimRedBlackGrid),
@ -70,7 +70,7 @@ namespace QCD {
pickCheckerboard(Even,UmuEven,Umu);
pickCheckerboard(Odd ,UmuOdd,Umu);
}
void FiveDimWilsonFermion::DoubleStore(LatticeDoubledGaugeField &Uds,const LatticeGaugeField &Umu)
void WilsonFermion5D::DoubleStore(LatticeDoubledGaugeField &Uds,const LatticeGaugeField &Umu)
{
conformable(Uds._grid,GaugeGrid());
conformable(Umu._grid,GaugeGrid());
@ -82,60 +82,9 @@ void FiveDimWilsonFermion::DoubleStore(LatticeDoubledGaugeField &Uds,const Latti
pokeIndex<LorentzIndex>(Uds,U,mu+4);
}
}
RealD FiveDimWilsonFermion::M(const LatticeFermion &in, LatticeFermion &out)
{
out.checkerboard=in.checkerboard;
Dhop(in,out,DaggerNo);
return axpy_norm(out,5.0-M5,in,out);
}
RealD FiveDimWilsonFermion::Mdag(const LatticeFermion &in, LatticeFermion &out)
{
out.checkerboard=in.checkerboard;
Dhop(in,out,DaggerYes);
return axpy_norm(out,5.0-M5,in,out);
}
void FiveDimWilsonFermion::Meooe(const LatticeFermion &in, LatticeFermion &out)
{
if ( in.checkerboard == Odd ) {
DhopEO(in,out,DaggerNo);
} else {
DhopOE(in,out,DaggerNo);
}
}
void FiveDimWilsonFermion::MeooeDag(const LatticeFermion &in, LatticeFermion &out)
{
if ( in.checkerboard == Odd ) {
DhopEO(in,out,DaggerYes);
} else {
DhopOE(in,out,DaggerYes);
}
}
void FiveDimWilsonFermion::Mooee(const LatticeFermion &in, LatticeFermion &out)
{
out.checkerboard = in.checkerboard;
out = (5.0-M5)*in;
return ;
}
void FiveDimWilsonFermion::MooeeDag(const LatticeFermion &in, LatticeFermion &out)
{
out.checkerboard = in.checkerboard;
Mooee(in,out);
}
void FiveDimWilsonFermion::MooeeInv(const LatticeFermion &in, LatticeFermion &out)
{
out.checkerboard = in.checkerboard;
out = (1.0/(5.0-M5))*in;
return ;
}
void FiveDimWilsonFermion::MooeeInvDag(const LatticeFermion &in, LatticeFermion &out)
{
out.checkerboard = in.checkerboard;
MooeeInv(in,out);
}
void FiveDimWilsonFermion::DhopInternal(CartesianStencil & st, LebesgueOrder &lo,
LatticeDoubledGaugeField & U,
const LatticeFermion &in, LatticeFermion &out,int dag)
void WilsonFermion5D::DhopInternal(CartesianStencil & st, LebesgueOrder &lo,
LatticeDoubledGaugeField & U,
const LatticeFermion &in, LatticeFermion &out,int dag)
{
assert((dag==DaggerNo) ||(dag==DaggerYes));
@ -150,19 +99,21 @@ void FiveDimWilsonFermion::DhopInternal(CartesianStencil & st, LebesgueOrder &lo
// - 8 linear access unit stride streams per thread for Fermion for hw prefetchable.
if ( dag == DaggerYes ) {
if( HandOptDslash ) {
for(int ss=0;ss<U._grid->oSites();ss++){
int sU=lo.Reorder(ss);
PARALLEL_FOR_LOOP
for(int ss=0;ss<U._grid->oSites();ss++){
for(int s=0;s<Ls;s++){
//int sU=lo.Reorder(ss);
int sU=ss;
int sF = s+Ls*sU;
DiracOptHand::DhopSiteDag(st,U,comm_buf,sF,sU,in,out);
}
}
} else {
for(int ss=0;ss<U._grid->oSites();ss++){
int sU=lo.Reorder(ss);
PARALLEL_FOR_LOOP
for(int ss=0;ss<U._grid->oSites();ss++){
for(int s=0;s<Ls;s++){
// int sU=lo.Reorder(ss);
int sU=ss;
int sF = s+Ls*sU;
DiracOpt::DhopSiteDag(st,U,comm_buf,sF,sU,in,out);
}
@ -170,21 +121,22 @@ PARALLEL_FOR_LOOP
}
} else {
if( HandOptDslash ) {
PARALLEL_FOR_LOOP
for(int ss=0;ss<U._grid->oSites();ss++){
int sU=lo.Reorder(ss);
for(int s=0;s<Ls;s++){
// int sU=lo.Reorder(ss);
int sU=ss;
int sF = s+Ls*sU;
DiracOptHand::DhopSite(st,U,comm_buf,sF,sU,in,out);
}
}
} else {
for(int ss=0;ss<U._grid->oSites();ss++){
int sU=lo.Reorder(ss);
PARALLEL_FOR_LOOP
for(int ss=0;ss<U._grid->oSites();ss++){
for(int s=0;s<Ls;s++){
// int sU=lo.Reorder(ss);
int sU=ss;
int sF = s+Ls*sU;
DiracOpt::DhopSite(st,U,comm_buf,sF,sU,in,out);
}
@ -192,7 +144,7 @@ PARALLEL_FOR_LOOP
}
}
}
void FiveDimWilsonFermion::DhopOE(const LatticeFermion &in, LatticeFermion &out,int dag)
void WilsonFermion5D::DhopOE(const LatticeFermion &in, LatticeFermion &out,int dag)
{
conformable(in._grid,FermionRedBlackGrid()); // verifies half grid
conformable(in._grid,out._grid); // drops the cb check
@ -202,7 +154,7 @@ void FiveDimWilsonFermion::DhopOE(const LatticeFermion &in, LatticeFermion &out,
DhopInternal(StencilEven,LebesgueEvenOdd,UmuOdd,in,out,dag);
}
void FiveDimWilsonFermion::DhopEO(const LatticeFermion &in, LatticeFermion &out,int dag)
void WilsonFermion5D::DhopEO(const LatticeFermion &in, LatticeFermion &out,int dag)
{
conformable(in._grid,FermionRedBlackGrid()); // verifies half grid
conformable(in._grid,out._grid); // drops the cb check
@ -212,7 +164,7 @@ void FiveDimWilsonFermion::DhopEO(const LatticeFermion &in, LatticeFermion &out,
DhopInternal(StencilOdd,LebesgueEvenOdd,UmuEven,in,out,dag);
}
void FiveDimWilsonFermion::Dhop(const LatticeFermion &in, LatticeFermion &out,int dag)
void WilsonFermion5D::Dhop(const LatticeFermion &in, LatticeFermion &out,int dag)
{
conformable(in._grid,FermionGrid()); // verifies full grid
conformable(in._grid,out._grid);
@ -221,8 +173,14 @@ void FiveDimWilsonFermion::Dhop(const LatticeFermion &in, LatticeFermion &out,in
DhopInternal(Stencil,Lebesgue,Umu,in,out,dag);
}
}}
void WilsonFermion5D::DW(const LatticeFermion &in, LatticeFermion &out,int dag)
{
out.checkerboard=in.checkerboard;
Dhop(in,out,dag); // -0.5 is included
axpy(out,4.0-M5,in,out);
}
}
}

47
lib/qcd/action/fermion/FiveDimWilsonFermion.h → lib/qcd/action/fermion/WilsonFermion5D.h
View File

@ -1,5 +1,5 @@
#ifndef GRID_QCD_DWF_H
#define GRID_QCD_DWF_H
#ifndef GRID_QCD_WILSON_FERMION_5D_H
#define GRID_QCD_WILSON_FERMION_5D_H
namespace Grid {
@ -14,8 +14,20 @@ namespace Grid {
// i.e. even even contains fifth dim hopping term.
//
// [DIFFERS from original CPS red black implementation parity = (x+y+z+t+s)|2 ]
////////////////////////////
//ContFrac:
// Ls always odd. Rational poly deg is either Ls or Ls-1
//PartFrac
// Ls always odd. Rational poly deg is either Ls or Ls-1
//
//Cayley: Ls always even, Rational poly deg is Ls
//
// Just set nrational as Ls. Forget about Ls-1 cases.
//
// Require odd Ls for cont and part frac
////////////////////////////
////////////////////////////////////////////////////////////////////////////////
class FiveDimWilsonFermion : public FermionAction<LatticeFermion,LatticeGaugeField>
class WilsonFermion5D : public FermionOperator<LatticeFermion,LatticeGaugeField>
{
public:
///////////////////////////////////////////////////////////////
@ -26,19 +38,21 @@ namespace Grid {
GridBase *FermionGrid(void) { return _FiveDimGrid;}
GridBase *FermionRedBlackGrid(void) { return _FiveDimRedBlackGrid;}
// override multiply
virtual RealD M (const LatticeFermion &in, LatticeFermion &out);
virtual RealD Mdag (const LatticeFermion &in, LatticeFermion &out);
// full checkerboard operations; leave unimplemented as abstract for now
//virtual RealD M (const LatticeFermion &in, LatticeFermion &out)=0;
//virtual RealD Mdag (const LatticeFermion &in, LatticeFermion &out)=0;
// half checkerboard operaions
virtual void Meooe (const LatticeFermion &in, LatticeFermion &out);
virtual void MeooeDag (const LatticeFermion &in, LatticeFermion &out);
virtual void Mooee (const LatticeFermion &in, LatticeFermion &out);
virtual void MooeeDag (const LatticeFermion &in, LatticeFermion &out);
virtual void MooeeInv (const LatticeFermion &in, LatticeFermion &out);
virtual void MooeeInvDag (const LatticeFermion &in, LatticeFermion &out);
// half checkerboard operations; leave unimplemented as abstract for now
// virtual void Meooe (const LatticeFermion &in, LatticeFermion &out)=0;
// virtual void MeooeDag (const LatticeFermion &in, LatticeFermion &out)=0;
// virtual void Mooee (const LatticeFermion &in, LatticeFermion &out)=0;
// virtual void MooeeDag (const LatticeFermion &in, LatticeFermion &out)=0;
// virtual void MooeeInv (const LatticeFermion &in, LatticeFermion &out)=0;
// virtual void MooeeInvDag (const LatticeFermion &in, LatticeFermion &out)=0;
// non-hermitian hopping term; half cb or both
// Implement hopping term non-hermitian hopping term; half cb or both
// Implement s-diagonal DW
void DW (const LatticeFermion &in, LatticeFermion &out,int dag);
void Dhop (const LatticeFermion &in, LatticeFermion &out,int dag);
void DhopOE(const LatticeFermion &in, LatticeFermion &out,int dag);
void DhopEO(const LatticeFermion &in, LatticeFermion &out,int dag);
@ -54,12 +68,12 @@ namespace Grid {
int dag);
// Constructors
FiveDimWilsonFermion(LatticeGaugeField &_Umu,
WilsonFermion5D(LatticeGaugeField &_Umu,
GridCartesian &FiveDimGrid,
GridRedBlackCartesian &FiveDimRedBlackGrid,
GridCartesian &FourDimGrid,
GridRedBlackCartesian &FourDimRedBlackGrid,
double _mass);
double _M5);
// DoubleStore
void DoubleStore(LatticeDoubledGaugeField &Uds,const LatticeGaugeField &Umu);
@ -82,7 +96,6 @@ namespace Grid {
static const std::vector<int> displacements;
double M5;
double mass;
int Ls;
//Defines the stencils for even and odd

10
lib/simd/Grid_avx.h
View File

@ -4,7 +4,7 @@
Using intrinsics
*/
// Time-stamp: <2015-05-27 12:07:15 neo>
// Time-stamp: <2015-05-29 14:13:30 neo>
//----------------------------------------------------------------------
#include <immintrin.h>
@ -261,13 +261,7 @@ namespace Optimization {
}
// Complex double
inline __m256d operator()(__m256d in){
return _mm256_xor_pd(_mm256_addsub_pd(_mm256_setzero_pd(),in), _mm256_set1_pd(-0.f));//untested
/*
// original
// addsubps 0, inv=>0+in.v[3] 0-in.v[2], 0+in.v[1], 0-in.v[0], ...
__m256d tmp = _mm256_addsub_pd(_mm256_setzero_pd(),_mm256_shuffle_pd(in,in,0x5));
return _mm256_shuffle_pd(tmp,tmp,0x5);
*/
return _mm256_xor_pd(_mm256_addsub_pd(_mm256_setzero_pd(),in), _mm256_set1_pd(-0.f));
}
// do not define for integer input
};

3
lib/simd/Grid_vector_types.h
View File

@ -2,7 +2,7 @@
/*! @file Grid_vector_types.h
@brief Defines templated class Grid_simd to deal with inner vector types
*/
// Time-stamp: <2015-05-27 12:04:06 neo>
// Time-stamp: <2015-05-29 14:19:48 neo>
//---------------------------------------------------------------------------
#ifndef GRID_VECTOR_TYPES
#define GRID_VECTOR_TYPES
@ -55,7 +55,6 @@ namespace Grid {
// general forms to allow for vsplat syntax
// need explicit declaration of types when used since
// clang cannot automatically determine the output type sometimes
// use decltype?
template < class Out, class Input1, class Input2, class Operation >
Out binary(Input1 src_1, Input2 src_2, Operation op){
return op(src_1, src_2);

2
lib/stamp-h1
View File

@ -1 +1 @@
timestamp for lib/Grid_config.h
timestamp for lib/GridConfig.h

0
lib/stencil/Grid_lebesgue.cc → lib/stencil/Lebesgue.cc
View File

0
lib/stencil/Grid_lebesgue.h → lib/stencil/Lebesgue.h
View File

12
lib/stencil/Grid_stencil_common.cc → lib/stencil/Stencil_common.cc
View File

@ -52,8 +52,8 @@ namespace Grid {
// up a table containing the npoint "neighbours" and whether they
// live in lattice or a comms buffer.
if ( !comm_dim ) {
sshift[0] = _grid->CheckerBoardShift(_checkerboard,dimension,shift,Even);
sshift[1] = _grid->CheckerBoardShift(_checkerboard,dimension,shift,Odd);
sshift[0] = _grid->CheckerBoardShiftForCB(_checkerboard,dimension,shift,Even);
sshift[1] = _grid->CheckerBoardShiftForCB(_checkerboard,dimension,shift,Odd);
if ( sshift[0] == sshift[1] ) {
Local(point,dimension,shift,0x3);
@ -63,8 +63,8 @@ namespace Grid {
}
} else { // All permute extract done in comms phase prior to Stencil application
// So tables are the same whether comm_dim or splice_dim
sshift[0] = _grid->CheckerBoardShift(_checkerboard,dimension,shift,Even);
sshift[1] = _grid->CheckerBoardShift(_checkerboard,dimension,shift,Odd);
sshift[0] = _grid->CheckerBoardShiftForCB(_checkerboard,dimension,shift,Even);
sshift[1] = _grid->CheckerBoardShiftForCB(_checkerboard,dimension,shift,Odd);
if ( sshift[0] == sshift[1] ) {
Comms(point,dimension,shift,0x3);
} else {
@ -96,7 +96,7 @@ namespace Grid {
int cb= (cbmask==0x2)? Odd : Even;
int sshift = _grid->CheckerBoardShift(_checkerboard,dimension,shift,cb);
int sshift = _grid->CheckerBoardShiftForCB(_checkerboard,dimension,shift,cb);
int sx = (x+sshift)%rd;
int permute_slice=0;
@ -134,7 +134,7 @@ namespace Grid {
// send to one or more remote nodes.
int cb= (cbmask==0x2)? Odd : Even;
int sshift= _grid->CheckerBoardShift(_checkerboard,dimension,shift,cb);
int sshift= _grid->CheckerBoardShiftForCB(_checkerboard,dimension,shift,cb);
for(int x=0;x<rd;x++){

43
lib/tensors/Tensor_Ta.h Normal file
View File

@ -0,0 +1,43 @@
#ifndef GRID_MATH_TA_H
#define GRID_MATH_TA_H
namespace Grid {
///////////////////////////////////////////////
// Ta function for scalar, vector, matrix
///////////////////////////////////////////////
inline ComplexF Ta( const ComplexF &arg){ return arg;}
inline ComplexD Ta( const ComplexD &arg){ return arg;}
inline RealF Ta( const RealF &arg){ return arg;}
inline RealD Ta( const RealD &arg){ return arg;}
template<class vtype> inline iScalar<vtype> Ta(const iScalar<vtype>&r)
{
iScalar<vtype> ret;
ret._internal = Ta(r._internal);
return ret;
}
template<class vtype,int N> inline iVector<vtype,N> Ta(const iVector<vtype,N>&r)
{
iVector<vtype,N> ret;
for(int i=0;i<N;i++){
ret._internal[i] = Ta(r._internal[i]);
}
return ret;
}
template<class vtype,int N> inline iMatrix<vtype,N> Ta(const iMatrix<vtype,N> &arg)
{
iMatrix<vtype,N> ret(arg);
double factor = (1/(double)N);
for(int c1=0;c1<N;c1++){
for(int c2=0;c2<N;c2++){
ret._internal[c1][c2]= (ret._internal[c1][c2] - adj(arg._internal[c2][c1]));
ret._internal[c1][c2] *= 0.5;
}}
//ret = (ret - adj(arg))*0.5;
ret -= trace(ret)*factor;
return ret;
}
}
#endif

11
lib/tensors/Tensor_arith.h Normal file
View File

@ -0,0 +1,11 @@
#ifndef GRID_MATH_ARITH_H
#define GRID_MATH_ARITH_H
#include <tensors/Tensor_arith_add.h>
#include <tensors/Tensor_arith_sub.h>
#include <tensors/Tensor_arith_mac.h>
#include <tensors/Tensor_arith_mul.h>
#include <tensors/Tensor_arith_scalar.h>
#endif

5
lib/math/Grid_math_arith_add.h → lib/tensors/Tensor_arith_add.h
View File

@ -45,7 +45,10 @@ namespace Grid {
{
for(int c2=0;c2<N;c2++){
for(int c1=0;c1<N;c1++){
add(&ret->_internal[c1][c2],&lhs->_internal,&rhs->_internal[c1][c2]);
if ( c1==c2)
add(&ret->_internal[c1][c2],&lhs->_internal,&rhs->_internal[c1][c2]);
else
ret->_internal[c1][c2]=lhs->_internal[c1][c2];
}}
return;
}

0
lib/math/Grid_math_arith_mac.h → lib/tensors/Tensor_arith_mac.h
View File

0
lib/math/Grid_math_arith_mul.h → lib/tensors/Tensor_arith_mul.h
View File

0
lib/math/Grid_math_arith_scalar.h → lib/tensors/Tensor_arith_scalar.h
View File

4
lib/math/Grid_math_arith_sub.h → lib/tensors/Tensor_arith_sub.h
View File

@ -44,7 +44,7 @@ template<class vtype,class ltype,class rtype, int N> strong_inline void sub(iMat
const iMatrix<rtype,N> * __restrict__ rhs){
for(int c2=0;c2<N;c2++){
for(int c1=0;c1<N;c1++){
if ( c1!=c2) {
if ( c1==c2) {
sub(&ret->_internal[c1][c2],&lhs->_internal,&rhs->_internal[c1][c2]);
} else {
// Fails -- need unary minus. Catalogue other unops?
@ -60,7 +60,7 @@ template<class vtype,class ltype,class rtype, int N> strong_inline void sub(iMat
const iScalar<rtype> * __restrict__ rhs){
for(int c2=0;c2<N;c2++){
for(int c1=0;c1<N;c1++){
if ( c1!=c2)
if ( c1==c2)
sub(&ret->_internal[c1][c2],&lhs->_internal[c1][c2],&rhs->_internal);
else
ret->_internal[c1][c2]=lhs->_internal[c1][c2];

0
lib/math/Grid_math_tensors.h → lib/tensors/Tensor_class.h
View File

0
lib/Grid_extract.h → lib/tensors/Tensor_extract_merge.h
View File

0
lib/math/Grid_math_inner.h → lib/tensors/Tensor_inner.h
View File

0
lib/math/Grid_math_outer.h → lib/tensors/Tensor_outer.h
View File

0
lib/math/Grid_math_peek.h → lib/tensors/Tensor_peek.h
View File

0
lib/math/Grid_math_poke.h → lib/tensors/Tensor_poke.h
View File

7
lib/math/Grid_math_reality.h → lib/tensors/Tensor_reality.h
View File

@ -2,10 +2,6 @@
#define GRID_MATH_REALITY_H
namespace Grid {
///////////////////////////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////// CONJ ///////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////
// multiply by I; make recursive.
///////////////////////////////////////////////
@ -151,6 +147,9 @@ template<class vtype,int N> inline iMatrix<vtype,N> adj(const iMatrix<vtype,N> &
/////////////////////////////////////////////////////////////////
// Can only take the real/imag part of scalar objects, since
// lattice objects of different complex nature are non-conformable.

13
lib/math/Grid_math_trace.h → lib/tensors/Tensor_trace.h
View File

@ -72,5 +72,18 @@ auto traceIndex(const iMatrix<vtype,N> &arg) -> iMatrix<decltype(traceIndex<Lev
return ret;
}
// Allow to recurse if vector, but never terminate on a vector
// trace of a different index can distribute across the vector index in a replicated way
// but we do not trace a vector index.
template<int Level,class vtype,int N,typename std::enable_if< iVector<vtype, N>::TensorLevel != Level >::type * =nullptr> inline
auto traceIndex(const iVector<vtype,N> &arg) -> iVector<decltype(traceIndex<Level>(arg._internal[0])),N>
{
iVector<decltype(traceIndex<Level>(arg._internal[0])),N> ret;
for(int i=0;i<N;i++){
ret._internal[i] = traceIndex<Level>(arg._internal[i]);
}
return ret;
}
}
#endif

Some files were not shown because too many files have changed in this diff Show More