mirror of
https://github.com/paboyle/Grid.git
synced 2024-09-20 01:05:38 +01:00
Merge pull request #210 from grid-test-organisation/feature/gpu-port-develop
Cayley fermion functions for GPUs
This commit is contained in:
Peter Boyle
GitHub
commit
ee6f96d85c
@ -220,7 +220,10 @@ public:
|
||||
#endif
|
||||
#endif
|
||||
}
|
||||
void construct(pointer __p, const _Tp& __val) { };
|
||||
|
||||
// FIXME: hack for the copy constructor, eventually it must be avoided
|
||||
void construct(pointer __p, const _Tp& __val) { new((void *)__p) _Tp(__val); };
|
||||
//void construct(pointer __p, const _Tp& __val) { };
|
||||
void construct(pointer __p) { };
|
||||
void destroy(pointer __p) { };
|
||||
};
|
||||
|
||||
@ -23,12 +23,7 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
|
||||
|
||||
#include <Grid/Grid_Eigen_Dense.h>
|
||||
#ifdef GRID_NVCC
|
||||
#include <thrust/host_vector.h>
|
||||
#include <thrust/device_vector.h>
|
||||
#include <thrust/generate.h>
|
||||
#include <thrust/reduce.h>
|
||||
#include <thrust/functional.h>
|
||||
#include <thrust/reduce.h>
|
||||
#include <thrust/inner_product.h>
|
||||
#endif
|
||||
|
||||
NAMESPACE_BEGIN(Grid);
|
||||
@ -41,23 +36,12 @@ template<class vobj> inline RealD norm2(const Lattice<vobj> &arg){
|
||||
return real(nrm);
|
||||
}
|
||||
|
||||
#if 0
|
||||
//#warning "ThrustReduce compiled"
|
||||
//#include <thrust/execution_policy.h>
|
||||
template<class vobj>
|
||||
vobj ThrustNorm(const Lattice<vobj> &lat)
|
||||
#ifdef GRID_NVCC
|
||||
template<class T, class R>
|
||||
struct innerProductFunctor : public thrust::binary_function<T,T,R>
|
||||
{
|
||||
typedef typename vobj::scalar_type scalar_type;
|
||||
auto lat_v=lat.View();
|
||||
Integer s0=0;
|
||||
Integer sN=lat_v.end();
|
||||
scalar_type sum = 0;
|
||||
scalar_type * begin = (scalar_type *)&lat_v[s0];
|
||||
scalar_type * end = (scalar_type *)&lat_v[sN];
|
||||
thrust::reduce(begin,end,sum);
|
||||
std::cout <<" thrust::reduce sum "<< sum << std::endl;
|
||||
return sum;
|
||||
}
|
||||
accelerator R operator()(T x, T y) { return innerProduct(x,y); }
|
||||
};
|
||||
#endif
|
||||
|
||||
// Double inner product
|
||||
@ -75,24 +59,17 @@ inline ComplexD innerProduct(const Lattice<vobj> &left,const Lattice<vobj> &righ
|
||||
auto left_v = left.View();
|
||||
auto right_v=right.View();
|
||||
|
||||
#if 0
|
||||
#ifdef GRID_NVCC
|
||||
|
||||
typedef decltype(TensorRemove(innerProduct(left_v[0],right_v[0]))) inner_t;
|
||||
typedef decltype(innerProduct(left_v[0],right_v[0])) inner_t;
|
||||
thrust::plus<inner_t> binary_sum;
|
||||
innerProductFunctor<vobj,inner_t> binary_inner_p;
|
||||
Integer sN = left_v.end();
|
||||
inner_t zero = Zero();
|
||||
// is there a way of using the efficient thrust reduction while maintaining memory coalescing?
|
||||
inner_t vnrm = thrust::inner_product(thrust::device, &left_v[0], &left_v[sN], &right_v[0], zero, binary_sum, binary_inner_p);
|
||||
nrm = Reduce(TensorRemove(vnrm));// sum across simd
|
||||
|
||||
Lattice<inner_t> inner_tmp(grid);
|
||||
|
||||
/////////////////////////
|
||||
// localInnerProduct
|
||||
/////////////////////////
|
||||
auto inner_tmp_v = inner_tmp.View();
|
||||
accelerator_loop(ss,left_v,{
|
||||
inner_tmp_v[ss] = TensorRemove(innerProduct(left_v[ss],right_v[ss]));
|
||||
});
|
||||
/////////////////////////
|
||||
// and site sum the scalars
|
||||
/////////////////////////
|
||||
inner_t vnrm = ThrustNorm(inner_tmp);
|
||||
auto vvnrm = vnrm;
|
||||
#else
|
||||
thread_loop( (int thr=0;thr<grid->SumArraySize();thr++),{
|
||||
int mywork, myoff;
|
||||
@ -108,9 +85,9 @@ inline ComplexD innerProduct(const Lattice<vobj> &left,const Lattice<vobj> &righ
|
||||
vector_type vvnrm; vvnrm=Zero(); // sum across threads
|
||||
for(int i=0;i<grid->SumArraySize();i++){
|
||||
vvnrm = vvnrm+sumarray[i];
|
||||
}
|
||||
#endif
|
||||
}
|
||||
nrm = Reduce(vvnrm);// sum across simd
|
||||
#endif
|
||||
right.Grid()->GlobalSum(nrm);
|
||||
return nrm;
|
||||
}
|
||||
|
||||
@ -163,10 +163,16 @@ template<class Impl> void CayleyFermion5D<Impl>::CayleyReport(void)
|
||||
std::cout << GridLogMessage << "CayleyFermion5D Number of M5D Calls : " << M5Dcalls << std::endl;
|
||||
std::cout << GridLogMessage << "CayleyFermion5D ComputeTime/Calls : " << M5Dtime / M5Dcalls << " us" << std::endl;
|
||||
|
||||
// Flops = 6.0*(Nc*Ns) *Ls*vol
|
||||
RealD mflops = 6.0*12*volume*M5Dcalls/M5Dtime/2; // 2 for red black counting
|
||||
// Flops = 10.0*(Nc*Ns) *Ls*vol
|
||||
RealD mflops = 10.0*(Nc*Ns)*volume*M5Dcalls/M5Dtime/2; // 2 for red black counting
|
||||
std::cout << GridLogMessage << "Average mflops/s per call : " << mflops << std::endl;
|
||||
std::cout << GridLogMessage << "Average mflops/s per call per rank : " << mflops/NP << std::endl;
|
||||
|
||||
// Bytes = sizeof(Real) * (Nc*Ns*Nreim) * Ls * vol * (read+write) (/2 for red black counting)
|
||||
// read = 2 ( psi[ss+s+1] and psi[ss+s-1] count as 1 )
|
||||
// write = 1
|
||||
RealD Gbytes = sizeof(Real) * (Nc*Ns*2) * volume * 3 /2. * 1.e-9;
|
||||
std::cout << GridLogMessage << "Average bandwidth (GB/s) : " << Gbytes/M5Dtime*M5Dcalls*1.e6 << std::endl;
|
||||
}
|
||||
|
||||
if ( MooeeInvCalls > 0 ) {
|
||||
@ -174,11 +180,16 @@ template<class Impl> void CayleyFermion5D<Impl>::CayleyReport(void)
|
||||
std::cout << GridLogMessage << "#### MooeeInv calls report " << std::endl;
|
||||
std::cout << GridLogMessage << "CayleyFermion5D Number of MooeeInv Calls : " << MooeeInvCalls << std::endl;
|
||||
std::cout << GridLogMessage << "CayleyFermion5D ComputeTime/Calls : " << MooeeInvTime / MooeeInvCalls << " us" << std::endl;
|
||||
|
||||
#ifdef GRID_NVCC
|
||||
RealD mflops = ( -16.*Nc*Ns+this->Ls*(1.+18.*Nc*Ns) )*volume*MooeeInvCalls/MooeeInvTime/2; // 2 for red black counting
|
||||
std::cout << GridLogMessage << "Average mflops/s per call : " << mflops << std::endl;
|
||||
std::cout << GridLogMessage << "Average mflops/s per call per rank : " << mflops/NP << std::endl;
|
||||
#else
|
||||
// Flops = MADD * Ls *Ls *4dvol * spin/colour/complex
|
||||
RealD mflops = 2.0*24*this->Ls*volume*MooeeInvCalls/MooeeInvTime/2; // 2 for red black counting
|
||||
std::cout << GridLogMessage << "Average mflops/s per call : " << mflops << std::endl;
|
||||
std::cout << GridLogMessage << "Average mflops/s per call per rank : " << mflops/NP << std::endl;
|
||||
#endif
|
||||
}
|
||||
|
||||
}
|
||||
@ -197,18 +208,18 @@ template<class Impl>
|
||||
void CayleyFermion5D<Impl>::M5D (const FermionField &psi, FermionField &chi)
|
||||
{
|
||||
int Ls=this->Ls;
|
||||
std::vector<Coeff_t> diag (Ls,1.0);
|
||||
std::vector<Coeff_t> upper(Ls,-1.0); upper[Ls-1]=mass;
|
||||
std::vector<Coeff_t> lower(Ls,-1.0); lower[0] =mass;
|
||||
Vector<Coeff_t> diag (Ls,1.0);
|
||||
Vector<Coeff_t> upper(Ls,-1.0); upper[Ls-1]=mass;
|
||||
Vector<Coeff_t> lower(Ls,-1.0); lower[0] =mass;
|
||||
M5D(psi,chi,chi,lower,diag,upper);
|
||||
}
|
||||
template<class Impl>
|
||||
void CayleyFermion5D<Impl>::Meooe5D (const FermionField &psi, FermionField &Din)
|
||||
{
|
||||
int Ls=this->Ls;
|
||||
std::vector<Coeff_t> diag = bs;
|
||||
std::vector<Coeff_t> upper= cs;
|
||||
std::vector<Coeff_t> lower= cs;
|
||||
Vector<Coeff_t> diag = bs;
|
||||
Vector<Coeff_t> upper= cs;
|
||||
Vector<Coeff_t> lower= cs;
|
||||
upper[Ls-1]=-mass*upper[Ls-1];
|
||||
lower[0] =-mass*lower[0];
|
||||
M5D(psi,psi,Din,lower,diag,upper);
|
||||
@ -217,9 +228,9 @@ void CayleyFermion5D<Impl>::Meooe5D (const FermionField &psi, FermionField &D
|
||||
template<class Impl> void CayleyFermion5D<Impl>::Meo5D (const FermionField &psi, FermionField &chi)
|
||||
{
|
||||
int Ls=this->Ls;
|
||||
std::vector<Coeff_t> diag = beo;
|
||||
std::vector<Coeff_t> upper(Ls);
|
||||
std::vector<Coeff_t> lower(Ls);
|
||||
Vector<Coeff_t> diag = beo;
|
||||
Vector<Coeff_t> upper(Ls);
|
||||
Vector<Coeff_t> lower(Ls);
|
||||
for(int i=0;i<Ls;i++) {
|
||||
upper[i]=-ceo[i];
|
||||
lower[i]=-ceo[i];
|
||||
@ -232,9 +243,9 @@ template<class Impl>
|
||||
void CayleyFermion5D<Impl>::Mooee (const FermionField &psi, FermionField &chi)
|
||||
{
|
||||
int Ls=this->Ls;
|
||||
std::vector<Coeff_t> diag = bee;
|
||||
std::vector<Coeff_t> upper(Ls);
|
||||
std::vector<Coeff_t> lower(Ls);
|
||||
Vector<Coeff_t> diag = bee;
|
||||
Vector<Coeff_t> upper(Ls);
|
||||
Vector<Coeff_t> lower(Ls);
|
||||
for(int i=0;i<Ls;i++) {
|
||||
upper[i]=-cee[i];
|
||||
lower[i]=-cee[i];
|
||||
@ -247,9 +258,9 @@ template<class Impl>
|
||||
void CayleyFermion5D<Impl>::MooeeDag (const FermionField &psi, FermionField &chi)
|
||||
{
|
||||
int Ls=this->Ls;
|
||||
std::vector<Coeff_t> diag = bee;
|
||||
std::vector<Coeff_t> upper(Ls);
|
||||
std::vector<Coeff_t> lower(Ls);
|
||||
Vector<Coeff_t> diag = bee;
|
||||
Vector<Coeff_t> upper(Ls);
|
||||
Vector<Coeff_t> lower(Ls);
|
||||
|
||||
for (int s=0;s<Ls;s++){
|
||||
// Assemble the 5d matrix
|
||||
@ -277,9 +288,9 @@ template<class Impl>
|
||||
void CayleyFermion5D<Impl>::M5Ddag (const FermionField &psi, FermionField &chi)
|
||||
{
|
||||
int Ls=this->Ls;
|
||||
std::vector<Coeff_t> diag(Ls,1.0);
|
||||
std::vector<Coeff_t> upper(Ls,-1.0);
|
||||
std::vector<Coeff_t> lower(Ls,-1.0);
|
||||
Vector<Coeff_t> diag(Ls,1.0);
|
||||
Vector<Coeff_t> upper(Ls,-1.0);
|
||||
Vector<Coeff_t> lower(Ls,-1.0);
|
||||
upper[Ls-1]=-mass*upper[Ls-1];
|
||||
lower[0] =-mass*lower[0];
|
||||
M5Ddag(psi,chi,chi,lower,diag,upper);
|
||||
@ -289,9 +300,9 @@ template<class Impl>
|
||||
void CayleyFermion5D<Impl>::MeooeDag5D (const FermionField &psi, FermionField &Din)
|
||||
{
|
||||
int Ls=this->Ls;
|
||||
std::vector<Coeff_t> diag =bs;
|
||||
std::vector<Coeff_t> upper=cs;
|
||||
std::vector<Coeff_t> lower=cs;
|
||||
Vector<Coeff_t> diag =bs;
|
||||
Vector<Coeff_t> upper=cs;
|
||||
Vector<Coeff_t> lower=cs;
|
||||
|
||||
for (int s=0;s<Ls;s++){
|
||||
if ( s== 0 ) {
|
||||
@ -428,7 +439,7 @@ void CayleyFermion5D<Impl>::MeoDeriv(GaugeField &mat,const FermionField &U,const
|
||||
template<class Impl>
|
||||
void CayleyFermion5D<Impl>::SetCoefficientsTanh(Approx::zolotarev_data *zdata,RealD b,RealD c)
|
||||
{
|
||||
std::vector<Coeff_t> gamma(this->Ls);
|
||||
Vector<Coeff_t> gamma(this->Ls);
|
||||
for(int s=0;s<this->Ls;s++) gamma[s] = zdata->gamma[s];
|
||||
SetCoefficientsInternal(1.0,gamma,b,c);
|
||||
}
|
||||
@ -436,13 +447,13 @@ void CayleyFermion5D<Impl>::SetCoefficientsTanh(Approx::zolotarev_data *zdata,Re
|
||||
template<class Impl>
|
||||
void CayleyFermion5D<Impl>::SetCoefficientsZolotarev(RealD zolo_hi,Approx::zolotarev_data *zdata,RealD b,RealD c)
|
||||
{
|
||||
std::vector<Coeff_t> gamma(this->Ls);
|
||||
Vector<Coeff_t> gamma(this->Ls);
|
||||
for(int s=0;s<this->Ls;s++) gamma[s] = zdata->gamma[s];
|
||||
SetCoefficientsInternal(zolo_hi,gamma,b,c);
|
||||
}
|
||||
//Zolo
|
||||
template<class Impl>
|
||||
void CayleyFermion5D<Impl>::SetCoefficientsInternal(RealD zolo_hi,std::vector<Coeff_t> & gamma,RealD b,RealD c)
|
||||
void CayleyFermion5D<Impl>::SetCoefficientsInternal(RealD zolo_hi,Vector<Coeff_t> & gamma,RealD b,RealD c)
|
||||
{
|
||||
int Ls=this->Ls;
|
||||
|
||||
|
||||
@ -108,16 +108,16 @@ public:
|
||||
void M5D(const FermionField &psi,
|
||||
const FermionField &phi,
|
||||
FermionField &chi,
|
||||
std::vector<Coeff_t> &lower,
|
||||
std::vector<Coeff_t> &diag,
|
||||
std::vector<Coeff_t> &upper);
|
||||
Vector<Coeff_t> &lower,
|
||||
Vector<Coeff_t> &diag,
|
||||
Vector<Coeff_t> &upper);
|
||||
|
||||
void M5Ddag(const FermionField &psi,
|
||||
const FermionField &phi,
|
||||
FermionField &chi,
|
||||
std::vector<Coeff_t> &lower,
|
||||
std::vector<Coeff_t> &diag,
|
||||
std::vector<Coeff_t> &upper);
|
||||
Vector<Coeff_t> &lower,
|
||||
Vector<Coeff_t> &diag,
|
||||
Vector<Coeff_t> &upper);
|
||||
|
||||
void MooeeInternal(const FermionField &in, FermionField &out,int dag,int inv);
|
||||
void MooeeInternalCompute(int dag, int inv, Vector<iSinglet<Simd> > & Matp, Vector<iSinglet<Simd> > & Matm);
|
||||
@ -149,29 +149,29 @@ public:
|
||||
RealD mass;
|
||||
|
||||
// Save arguments to SetCoefficientsInternal
|
||||
std::vector<Coeff_t> _gamma;
|
||||
Vector<Coeff_t> _gamma;
|
||||
RealD _zolo_hi;
|
||||
RealD _b;
|
||||
RealD _c;
|
||||
|
||||
// Cayley form Moebius (tanh and zolotarev)
|
||||
std::vector<Coeff_t> omega;
|
||||
std::vector<Coeff_t> bs; // S dependent coeffs
|
||||
std::vector<Coeff_t> cs;
|
||||
std::vector<Coeff_t> as;
|
||||
Vector<Coeff_t> omega;
|
||||
Vector<Coeff_t> bs; // S dependent coeffs
|
||||
Vector<Coeff_t> cs;
|
||||
Vector<Coeff_t> as;
|
||||
// For preconditioning Cayley form
|
||||
std::vector<Coeff_t> bee;
|
||||
std::vector<Coeff_t> cee;
|
||||
std::vector<Coeff_t> aee;
|
||||
std::vector<Coeff_t> beo;
|
||||
std::vector<Coeff_t> ceo;
|
||||
std::vector<Coeff_t> aeo;
|
||||
Vector<Coeff_t> bee;
|
||||
Vector<Coeff_t> cee;
|
||||
Vector<Coeff_t> aee;
|
||||
Vector<Coeff_t> beo;
|
||||
Vector<Coeff_t> ceo;
|
||||
Vector<Coeff_t> aeo;
|
||||
// LDU factorisation of the eeoo matrix
|
||||
std::vector<Coeff_t> lee;
|
||||
std::vector<Coeff_t> leem;
|
||||
std::vector<Coeff_t> uee;
|
||||
std::vector<Coeff_t> ueem;
|
||||
std::vector<Coeff_t> dee;
|
||||
Vector<Coeff_t> lee;
|
||||
Vector<Coeff_t> leem;
|
||||
Vector<Coeff_t> uee;
|
||||
Vector<Coeff_t> ueem;
|
||||
Vector<Coeff_t> dee;
|
||||
|
||||
// Matrices of 5d ee inverse params
|
||||
Vector<iSinglet<Simd> > MatpInv;
|
||||
@ -203,22 +203,26 @@ public:
|
||||
protected:
|
||||
virtual void SetCoefficientsZolotarev(RealD zolohi,Approx::zolotarev_data *zdata,RealD b,RealD c);
|
||||
virtual void SetCoefficientsTanh(Approx::zolotarev_data *zdata,RealD b,RealD c);
|
||||
virtual void SetCoefficientsInternal(RealD zolo_hi,std::vector<Coeff_t> & gamma,RealD b,RealD c);
|
||||
virtual void SetCoefficientsInternal(RealD zolo_hi,Vector<Coeff_t> & gamma,RealD b,RealD c);
|
||||
};
|
||||
|
||||
NAMESPACE_END(Grid);
|
||||
|
||||
#define INSTANTIATE_DPERP(A) \
|
||||
template void CayleyFermion5D< A >::M5D(const FermionField &psi,const FermionField &phi,FermionField &chi, \
|
||||
std::vector<Coeff_t> &lower,std::vector<Coeff_t> &diag,std::vector<Coeff_t> &upper); \
|
||||
Vector<Coeff_t> &lower,Vector<Coeff_t> &diag,Vector<Coeff_t> &upper); \
|
||||
template void CayleyFermion5D< A >::M5Ddag(const FermionField &psi,const FermionField &phi,FermionField &chi, \
|
||||
std::vector<Coeff_t> &lower,std::vector<Coeff_t> &diag,std::vector<Coeff_t> &upper); \
|
||||
Vector<Coeff_t> &lower,Vector<Coeff_t> &diag,Vector<Coeff_t> &upper); \
|
||||
template void CayleyFermion5D< A >::MooeeInv (const FermionField &psi, FermionField &chi); \
|
||||
template void CayleyFermion5D< A >::MooeeInvDag (const FermionField &psi, FermionField &chi);
|
||||
|
||||
#ifdef GRID_NVCC
|
||||
#define CAYLEY_DPERP_GPU
|
||||
#else
|
||||
#undef CAYLEY_DPERP_DENSE
|
||||
#define CAYLEY_DPERP_CACHE
|
||||
#undef CAYLEY_DPERP_LINALG
|
||||
#endif
|
||||
#define CAYLEY_DPERP_VEC
|
||||
|
||||
#endif
|
||||
|
||||
@ -41,9 +41,9 @@ template<class Impl>
|
||||
void CayleyFermion5D<Impl>::M5D(const FermionField &psi_i,
|
||||
const FermionField &phi_i,
|
||||
FermionField &chi_i,
|
||||
std::vector<Coeff_t> &lower,
|
||||
std::vector<Coeff_t> &diag,
|
||||
std::vector<Coeff_t> &upper)
|
||||
Vector<Coeff_t> &lower,
|
||||
Vector<Coeff_t> &diag,
|
||||
Vector<Coeff_t> &upper)
|
||||
{
|
||||
chi_i.Checkerboard()=psi_i.Checkerboard();
|
||||
GridBase *grid=psi_i.Grid();
|
||||
@ -52,7 +52,8 @@ void CayleyFermion5D<Impl>::M5D(const FermionField &psi_i,
|
||||
auto chi = chi_i.View();
|
||||
int Ls =this->Ls;
|
||||
assert(phi.Checkerboard() == psi.Checkerboard());
|
||||
// Flops = 6.0*(Nc*Ns) *Ls*vol
|
||||
// 10 = 3 complex mult + 2 complex add
|
||||
// Flops = 10.0*(Nc*Ns) *Ls*vol (/2 for red black counting)
|
||||
M5Dcalls++;
|
||||
M5Dtime-=usecond();
|
||||
|
||||
@ -87,9 +88,9 @@ template<class Impl>
|
||||
void CayleyFermion5D<Impl>::M5Ddag(const FermionField &psi_i,
|
||||
const FermionField &phi_i,
|
||||
FermionField &chi_i,
|
||||
std::vector<Coeff_t> &lower,
|
||||
std::vector<Coeff_t> &diag,
|
||||
std::vector<Coeff_t> &upper)
|
||||
Vector<Coeff_t> &lower,
|
||||
Vector<Coeff_t> &diag,
|
||||
Vector<Coeff_t> &upper)
|
||||
{
|
||||
chi_i.Checkerboard()=psi_i.Checkerboard();
|
||||
GridBase *grid=psi_i.Grid();
|
||||
|
||||
284
Grid/qcd/action/fermion/CayleyFermion5Dgpu.cc
Normal file
284
Grid/qcd/action/fermion/CayleyFermion5Dgpu.cc
Normal file
@ -0,0 +1,284 @@
|
||||
/*************************************************************************************
|
||||
|
||||
Grid physics library, www.github.com/paboyle/Grid
|
||||
|
||||
Source file: ./lib/qcd/action/fermion/CayleyFermion5D.cc
|
||||
|
||||
Copyright (C) 2015
|
||||
|
||||
Author: Peter Boyle <pabobyle@ph.ed.ac.uk>
|
||||
Author: Peter Boyle <paboyle@ph.ed.ac.uk>
|
||||
Author: Peter Boyle <peterboyle@Peters-MacBook-Pro-2.local>
|
||||
Author: paboyle <paboyle@ph.ed.ac.uk>
|
||||
|
||||
This program is free software; you can redistribute it and/or modify
|
||||
it under the terms of the GNU General Public License as published by
|
||||
the Free Software Foundation; either version 2 of the License, or
|
||||
(at your option) any later version.
|
||||
|
||||
This program is distributed in the hope that it will be useful,
|
||||
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||
GNU General Public License for more details.
|
||||
|
||||
You should have received a copy of the GNU General Public License along
|
||||
with this program; if not, write to the Free Software Foundation, Inc.,
|
||||
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
|
||||
|
||||
See the full license in the file "LICENSE" in the top level distribution directory
|
||||
*************************************************************************************/
|
||||
/* END LEGAL */
|
||||
|
||||
#include <Grid/qcd/action/fermion/FermionCore.h>
|
||||
#include <Grid/qcd/action/fermion/CayleyFermion5D.h>
|
||||
|
||||
|
||||
NAMESPACE_BEGIN(Grid);
|
||||
|
||||
// Pminus fowards
|
||||
// Pplus backwards..
|
||||
template<class Impl>
|
||||
void CayleyFermion5D<Impl>::M5D(const FermionField &psi_i,
|
||||
const FermionField &phi_i,
|
||||
FermionField &chi_i,
|
||||
Vector<Coeff_t> &lower,
|
||||
Vector<Coeff_t> &diag,
|
||||
Vector<Coeff_t> &upper)
|
||||
{
|
||||
chi_i.Checkerboard()=psi_i.Checkerboard();
|
||||
GridBase *grid=psi_i.Grid();
|
||||
auto psi = psi_i.View();
|
||||
auto phi = phi_i.View();
|
||||
auto chi = chi_i.View();
|
||||
Coeff_t *lower_v = &lower[0];
|
||||
Coeff_t *diag_v = &diag[0];
|
||||
Coeff_t *upper_v = &upper[0];
|
||||
int Ls =this->Ls;
|
||||
assert(phi.Checkerboard() == psi.Checkerboard());
|
||||
|
||||
const uint64_t nsimd = grid->Nsimd();
|
||||
const uint64_t sites4d = nsimd * grid->oSites() / Ls;
|
||||
|
||||
// 10 = 3 complex mult + 2 complex add
|
||||
// Flops = 10.0*(Nc*Ns) *Ls*vol (/2 for red black counting)
|
||||
M5Dcalls++;
|
||||
M5Dtime-=usecond();
|
||||
|
||||
accelerator_loopN( sss, sites4d ,{
|
||||
uint64_t lane = sss % nsimd;
|
||||
uint64_t ss = Ls * (sss / nsimd);
|
||||
|
||||
for(int s=0;s<Ls;s++){
|
||||
auto res = extractLane(lane,phi[ss+s]);
|
||||
res = diag_v[s]*res;
|
||||
|
||||
auto tmp = extractLane(lane,psi[ss+(s+1)%Ls]);
|
||||
spProj5m(tmp,tmp);
|
||||
res += upper_v[s]*tmp;
|
||||
|
||||
tmp = extractLane(lane,psi[ss+(s+Ls-1)%Ls]);
|
||||
spProj5p(tmp,tmp);
|
||||
res += lower_v[s]*tmp;
|
||||
|
||||
insertLane(lane,chi[ss+s],res);
|
||||
}
|
||||
});
|
||||
M5Dtime+=usecond();
|
||||
}
|
||||
|
||||
template<class Impl>
|
||||
void CayleyFermion5D<Impl>::M5Ddag(const FermionField &psi_i,
|
||||
const FermionField &phi_i,
|
||||
FermionField &chi_i,
|
||||
Vector<Coeff_t> &lower,
|
||||
Vector<Coeff_t> &diag,
|
||||
Vector<Coeff_t> &upper)
|
||||
{
|
||||
chi_i.Checkerboard()=psi_i.Checkerboard();
|
||||
GridBase *grid=psi_i.Grid();
|
||||
auto psi = psi_i.View();
|
||||
auto phi = phi_i.View();
|
||||
auto chi = chi_i.View();
|
||||
Coeff_t *lower_v = &lower[0];
|
||||
Coeff_t *diag_v = &diag[0];
|
||||
Coeff_t *upper_v = &upper[0];
|
||||
int Ls =this->Ls;
|
||||
assert(phi.Checkerboard() == psi.Checkerboard());
|
||||
|
||||
const uint64_t nsimd = grid->Nsimd();
|
||||
const uint64_t sites4d = nsimd * grid->oSites() / Ls;
|
||||
|
||||
// 10 = 3 complex mult + 2 complex add
|
||||
// Flops = 10.0*(Nc*Ns) *Ls*vol (/2 for red black counting)
|
||||
M5Dcalls++;
|
||||
M5Dtime-=usecond();
|
||||
|
||||
accelerator_loopN( sss, sites4d ,{
|
||||
uint64_t lane = sss % nsimd;
|
||||
uint64_t ss = Ls * (sss / nsimd);
|
||||
|
||||
for(int s=0;s<Ls;s++){
|
||||
auto res = extractLane(lane,phi[ss+s]);
|
||||
res = diag_v[s]*res;
|
||||
|
||||
auto tmp = extractLane(lane,psi[ss+(s+1)%Ls]);
|
||||
spProj5p(tmp,tmp);
|
||||
res += upper_v[s]*tmp;
|
||||
|
||||
tmp = extractLane(lane,psi[ss+(s+Ls-1)%Ls]);
|
||||
spProj5m(tmp,tmp);
|
||||
res += lower_v[s]*tmp;
|
||||
|
||||
insertLane(lane,chi[ss+s],res);
|
||||
}
|
||||
});
|
||||
M5Dtime+=usecond();
|
||||
}
|
||||
|
||||
template<class Impl>
|
||||
void CayleyFermion5D<Impl>::MooeeInv (const FermionField &psi_i, FermionField &chi_i)
|
||||
{
|
||||
chi_i.Checkerboard()=psi_i.Checkerboard();
|
||||
GridBase *grid=psi_i.Grid();
|
||||
|
||||
auto psi = psi_i.View();
|
||||
auto chi = chi_i.View();
|
||||
Coeff_t *lee_v = &lee[0];
|
||||
Coeff_t *leem_v = &leem[0];
|
||||
Coeff_t *uee_v = &uee[0];
|
||||
Coeff_t *ueem_v = &ueem[0];
|
||||
Coeff_t *dee_v = &dee[0];
|
||||
|
||||
int Ls=this->Ls;
|
||||
const uint64_t nsimd = grid->Nsimd();
|
||||
const uint64_t sites4d = nsimd * grid->oSites() / Ls;
|
||||
|
||||
typedef typename SiteSpinor::scalar_object ScalarSiteSpinor;
|
||||
|
||||
MooeeInvCalls++;
|
||||
MooeeInvTime-=usecond();
|
||||
|
||||
accelerator_loopN( sss, sites4d ,{
|
||||
uint64_t lane = sss % nsimd;
|
||||
uint64_t ss = Ls * (sss / nsimd);
|
||||
ScalarSiteSpinor res, tmp, acc;
|
||||
|
||||
// X = Nc*Ns
|
||||
// flops = 2X + (Ls-2)(4X + 4X) + 6X + 1 + 2X + (Ls-1)(10X + 1) = -16X + Ls(1+18X) = -192 + 217*Ls flops
|
||||
// Apply (L^{\prime})^{-1} L_m^{-1}
|
||||
res = extractLane(lane,psi[ss]);
|
||||
spProj5m(tmp,res);
|
||||
acc = leem_v[0]*tmp;
|
||||
spProj5p(tmp,res);
|
||||
insertLane(lane,chi[ss],res);
|
||||
|
||||
for(int s=1;s<Ls-1;s++){
|
||||
res = extractLane(lane,psi[ss+s]);
|
||||
res -= lee_v[s-1]*tmp;
|
||||
spProj5m(tmp,res);
|
||||
acc += leem_v[s]*tmp;
|
||||
spProj5p(tmp,res);
|
||||
insertLane(lane,chi[ss+s],res);
|
||||
}
|
||||
res = extractLane(lane,psi[ss+Ls-1]);
|
||||
res = res - lee_v[Ls-2]*tmp - acc;
|
||||
|
||||
// Apply U_m^{-1} D^{-1} U^{-1}
|
||||
res = (1.0/dee_v[Ls-1])*res;
|
||||
insertLane(lane,chi[ss+Ls-1],res);
|
||||
spProj5p(acc,res);
|
||||
spProj5m(tmp,res);
|
||||
for (int s=Ls-2;s>=0;s--){
|
||||
res = extractLane(lane,chi[ss+s]);
|
||||
res = (1.0/dee_v[s])*res - uee_v[s]*tmp - ueem_v[s]*acc;
|
||||
spProj5m(tmp,res);
|
||||
insertLane(lane,chi[ss+s],res);
|
||||
}
|
||||
});
|
||||
|
||||
MooeeInvTime+=usecond();
|
||||
|
||||
}
|
||||
|
||||
template<class Impl>
|
||||
void CayleyFermion5D<Impl>::MooeeInvDag (const FermionField &psi_i, FermionField &chi_i)
|
||||
{
|
||||
chi_i.Checkerboard()=psi_i.Checkerboard();
|
||||
GridBase *grid=psi_i.Grid();
|
||||
|
||||
auto psi = psi_i.View();
|
||||
auto chi = chi_i.View();
|
||||
Coeff_t *lee_v = &lee[0];
|
||||
Coeff_t *leem_v = &leem[0];
|
||||
Coeff_t *uee_v = &uee[0];
|
||||
Coeff_t *ueem_v = &ueem[0];
|
||||
Coeff_t *dee_v = &dee[0];
|
||||
|
||||
int Ls=this->Ls;
|
||||
const uint64_t nsimd = grid->Nsimd();
|
||||
const uint64_t sites4d = nsimd * grid->oSites() / Ls;
|
||||
|
||||
typedef typename SiteSpinor::scalar_object ScalarSiteSpinor;
|
||||
|
||||
MooeeInvCalls++;
|
||||
MooeeInvTime-=usecond();
|
||||
|
||||
accelerator_loopN( sss, sites4d ,{
|
||||
uint64_t lane = sss % nsimd;
|
||||
uint64_t ss = Ls * (sss / nsimd);
|
||||
ScalarSiteSpinor res, tmp, acc;
|
||||
|
||||
// X = Nc*Ns
|
||||
// flops = 2X + (Ls-2)(4X + 4X) + 6X + 1 + 2X + (Ls-1)(10X + 1) = -16X + Ls(1+18X) = -192 + 217*Ls flops
|
||||
// Apply (U^{\prime})^{-dagger} U_m^{-\dagger}
|
||||
res = extractLane(lane,psi[ss]);
|
||||
spProj5p(tmp,res);
|
||||
acc = conjugate(ueem_v[0])*tmp;
|
||||
spProj5m(tmp,res);
|
||||
insertLane(lane,chi[ss],res);
|
||||
|
||||
for(int s=1;s<Ls-1;s++){
|
||||
res = extractLane(lane,psi[ss+s]);
|
||||
res -= conjugate(uee_v[s-1])*tmp;
|
||||
spProj5p(tmp,res);
|
||||
acc += conjugate(ueem_v[s])*tmp;
|
||||
spProj5m(tmp,res);
|
||||
insertLane(lane,chi[ss+s],res);
|
||||
}
|
||||
res = extractLane(lane,psi[ss+Ls-1]);
|
||||
res = res - conjugate(uee_v[Ls-2])*tmp - acc;
|
||||
|
||||
// Apply L_m^{-\dagger} D^{-dagger} L^{-dagger}
|
||||
res = conjugate(1.0/dee_v[Ls-1])*res;
|
||||
insertLane(lane,chi[ss+Ls-1],res);
|
||||
spProj5m(acc,res);
|
||||
spProj5p(tmp,res);
|
||||
for (int s=Ls-2;s>=0;s--){
|
||||
res = extractLane(lane,chi[ss+s]);
|
||||
res = conjugate(1.0/dee_v[s])*res - conjugate(lee_v[s])*tmp - conjugate(leem_v[s])*acc;
|
||||
spProj5p(tmp,res);
|
||||
insertLane(lane,chi[ss+s],res);
|
||||
}
|
||||
});
|
||||
|
||||
MooeeInvTime+=usecond();
|
||||
|
||||
}
|
||||
|
||||
#ifdef CAYLEY_DPERP_GPU
|
||||
INSTANTIATE_DPERP(WilsonImplF);
|
||||
INSTANTIATE_DPERP(WilsonImplD);
|
||||
INSTANTIATE_DPERP(GparityWilsonImplF);
|
||||
INSTANTIATE_DPERP(GparityWilsonImplD);
|
||||
INSTANTIATE_DPERP(ZWilsonImplF);
|
||||
INSTANTIATE_DPERP(ZWilsonImplD);
|
||||
|
||||
INSTANTIATE_DPERP(WilsonImplFH);
|
||||
INSTANTIATE_DPERP(WilsonImplDF);
|
||||
INSTANTIATE_DPERP(GparityWilsonImplFH);
|
||||
INSTANTIATE_DPERP(GparityWilsonImplDF);
|
||||
INSTANTIATE_DPERP(ZWilsonImplFH);
|
||||
INSTANTIATE_DPERP(ZWilsonImplDF);
|
||||
#endif
|
||||
|
||||
NAMESPACE_END(Grid);
|
||||
@ -41,9 +41,9 @@ template<class Impl>
|
||||
void CayleyFermion5D<Impl>::M5D(const FermionField &psi,
|
||||
const FermionField &phi,
|
||||
FermionField &chi,
|
||||
std::vector<Coeff_t> &lower,
|
||||
std::vector<Coeff_t> &diag,
|
||||
std::vector<Coeff_t> &upper)
|
||||
Vector<Coeff_t> &lower,
|
||||
Vector<Coeff_t> &diag,
|
||||
Vector<Coeff_t> &upper)
|
||||
{
|
||||
Coeff_t one(1.0);
|
||||
int Ls=this->Ls;
|
||||
@ -64,9 +64,9 @@ template<class Impl>
|
||||
void CayleyFermion5D<Impl>::M5Ddag(const FermionField &psi,
|
||||
const FermionField &phi,
|
||||
FermionField &chi,
|
||||
std::vector<Coeff_t> &lower,
|
||||
std::vector<Coeff_t> &diag,
|
||||
std::vector<Coeff_t> &upper)
|
||||
Vector<Coeff_t> &lower,
|
||||
Vector<Coeff_t> &diag,
|
||||
Vector<Coeff_t> &upper)
|
||||
{
|
||||
Coeff_t one(1.0);
|
||||
int Ls=this->Ls;
|
||||
|
||||
@ -54,9 +54,9 @@ template<class Impl>
|
||||
void CayleyFermion5D<Impl>::M5D(const FermionField &psi_i,
|
||||
const FermionField &phi_i,
|
||||
FermionField &chi_i,
|
||||
std::vector<Coeff_t> &lower,
|
||||
std::vector<Coeff_t> &diag,
|
||||
std::vector<Coeff_t> &upper)
|
||||
Vector<Coeff_t> &lower,
|
||||
Vector<Coeff_t> &diag,
|
||||
Vector<Coeff_t> &upper)
|
||||
{
|
||||
chi_i.Checkerboard()=psi_i.Checkerboard();
|
||||
GridBase *grid=psi_i.Grid();
|
||||
@ -200,9 +200,9 @@ template<class Impl>
|
||||
void CayleyFermion5D<Impl>::M5Ddag(const FermionField &psi_i,
|
||||
const FermionField &phi_i,
|
||||
FermionField &chi_i,
|
||||
std::vector<Coeff_t> &lower,
|
||||
std::vector<Coeff_t> &diag,
|
||||
std::vector<Coeff_t> &upper)
|
||||
Vector<Coeff_t> &lower,
|
||||
Vector<Coeff_t> &diag,
|
||||
Vector<Coeff_t> &upper)
|
||||
{
|
||||
chi_i.Checkerboard()=psi_i.Checkerboard();
|
||||
GridBase *grid=psi_i.Grid();
|
||||
|
||||
@ -89,12 +89,12 @@ protected:
|
||||
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;
|
||||
Vector<double> Beta;
|
||||
Vector<double> cc;;
|
||||
Vector<double> cc_d;;
|
||||
Vector<double> sqrt_cc;
|
||||
Vector<double> See;
|
||||
Vector<double> Aee;
|
||||
|
||||
};
|
||||
|
||||
|
||||
@ -131,9 +131,9 @@ void DomainWallEOFAFermion<Impl>::M5D(const FermionField& psi, FermionField& chi
|
||||
else{ shiftm = -shift*(mq3-mq2); }
|
||||
}
|
||||
|
||||
std::vector<Coeff_t> diag(Ls,1.0);
|
||||
std::vector<Coeff_t> upper(Ls,-1.0); upper[Ls-1] = mq1 + shiftm;
|
||||
std::vector<Coeff_t> lower(Ls,-1.0); lower[0] = mq1 + shiftp;
|
||||
Vector<Coeff_t> diag(Ls,1.0);
|
||||
Vector<Coeff_t> upper(Ls,-1.0); upper[Ls-1] = mq1 + shiftm;
|
||||
Vector<Coeff_t> lower(Ls,-1.0); lower[0] = mq1 + shiftp;
|
||||
|
||||
#if(0)
|
||||
std::cout << GridLogMessage << "DomainWallEOFAFermion::M5D(FF&,FF&):" << std::endl;
|
||||
@ -168,9 +168,9 @@ void DomainWallEOFAFermion<Impl>::M5Ddag(const FermionField& psi, FermionField&
|
||||
else{ shiftm = -shift*(mq3-mq2); }
|
||||
}
|
||||
|
||||
std::vector<Coeff_t> diag(Ls,1.0);
|
||||
std::vector<Coeff_t> upper(Ls,-1.0); upper[Ls-1] = mq1 + shiftp;
|
||||
std::vector<Coeff_t> lower(Ls,-1.0); lower[0] = mq1 + shiftm;
|
||||
Vector<Coeff_t> diag(Ls,1.0);
|
||||
Vector<Coeff_t> upper(Ls,-1.0); upper[Ls-1] = mq1 + shiftp;
|
||||
Vector<Coeff_t> lower(Ls,-1.0); lower[0] = mq1 + shiftm;
|
||||
|
||||
#if(0)
|
||||
std::cout << GridLogMessage << "DomainWallEOFAFermion::M5Ddag(FF&,FF&):" << std::endl;
|
||||
@ -194,9 +194,9 @@ void DomainWallEOFAFermion<Impl>::Mooee(const FermionField& psi, FermionField& c
|
||||
{
|
||||
int Ls = this->Ls;
|
||||
|
||||
std::vector<Coeff_t> diag = this->bee;
|
||||
std::vector<Coeff_t> upper(Ls);
|
||||
std::vector<Coeff_t> lower(Ls);
|
||||
Vector<Coeff_t> diag = this->bee;
|
||||
Vector<Coeff_t> upper(Ls);
|
||||
Vector<Coeff_t> lower(Ls);
|
||||
|
||||
for(int s=0; s<Ls; s++){
|
||||
upper[s] = -this->cee[s];
|
||||
@ -213,9 +213,9 @@ void DomainWallEOFAFermion<Impl>::MooeeDag(const FermionField& psi, FermionField
|
||||
{
|
||||
int Ls = this->Ls;
|
||||
|
||||
std::vector<Coeff_t> diag = this->bee;
|
||||
std::vector<Coeff_t> upper(Ls);
|
||||
std::vector<Coeff_t> lower(Ls);
|
||||
Vector<Coeff_t> diag = this->bee;
|
||||
Vector<Coeff_t> upper(Ls);
|
||||
Vector<Coeff_t> lower(Ls);
|
||||
|
||||
for(int s=0; s<Ls; s++){
|
||||
upper[s] = -this->cee[s];
|
||||
@ -231,7 +231,7 @@ void DomainWallEOFAFermion<Impl>::MooeeDag(const FermionField& psi, FermionField
|
||||
|
||||
//Zolo
|
||||
template<class Impl>
|
||||
void DomainWallEOFAFermion<Impl>::SetCoefficientsInternal(RealD zolo_hi, std::vector<Coeff_t>& gamma, RealD b, RealD c)
|
||||
void DomainWallEOFAFermion<Impl>::SetCoefficientsInternal(RealD zolo_hi, Vector<Coeff_t>& gamma, RealD b, RealD c)
|
||||
{
|
||||
int Ls = this->Ls;
|
||||
int pm = this->pm;
|
||||
|
||||
@ -70,10 +70,10 @@ public:
|
||||
// Instantiate different versions depending on Impl
|
||||
/////////////////////////////////////////////////////
|
||||
void M5D(const FermionField& psi, const FermionField& phi, FermionField& chi,
|
||||
std::vector<Coeff_t>& lower, std::vector<Coeff_t>& diag, std::vector<Coeff_t>& upper);
|
||||
Vector<Coeff_t>& lower, Vector<Coeff_t>& diag, Vector<Coeff_t>& upper);
|
||||
|
||||
void M5Ddag(const FermionField& psi, const FermionField& phi, FermionField& chi,
|
||||
std::vector<Coeff_t>& lower, std::vector<Coeff_t>& diag, std::vector<Coeff_t>& upper);
|
||||
Vector<Coeff_t>& lower, Vector<Coeff_t>& diag, Vector<Coeff_t>& upper);
|
||||
|
||||
void MooeeInternal(const FermionField& in, FermionField& out, int dag, int inv);
|
||||
|
||||
@ -94,16 +94,16 @@ public:
|
||||
RealD _M5, const ImplParams& p=ImplParams());
|
||||
|
||||
protected:
|
||||
void SetCoefficientsInternal(RealD zolo_hi, std::vector<Coeff_t>& gamma, RealD b, RealD c);
|
||||
void SetCoefficientsInternal(RealD zolo_hi, Vector<Coeff_t>& gamma, RealD b, RealD c);
|
||||
};
|
||||
|
||||
NAMESPACE_END(Grid);
|
||||
|
||||
#define INSTANTIATE_DPERP_DWF_EOFA(A) \
|
||||
template void DomainWallEOFAFermion<A>::M5D(const FermionField& psi, const FermionField& phi, FermionField& chi, \
|
||||
std::vector<Coeff_t>& lower, std::vector<Coeff_t>& diag, std::vector<Coeff_t>& upper); \
|
||||
Vector<Coeff_t>& lower, Vector<Coeff_t>& diag, Vector<Coeff_t>& upper); \
|
||||
template void DomainWallEOFAFermion<A>::M5Ddag(const FermionField& psi, const FermionField& phi, FermionField& chi, \
|
||||
std::vector<Coeff_t>& lower, std::vector<Coeff_t>& diag, std::vector<Coeff_t>& upper); \
|
||||
Vector<Coeff_t>& lower, Vector<Coeff_t>& diag, Vector<Coeff_t>& upper); \
|
||||
template void DomainWallEOFAFermion<A>::MooeeInv(const FermionField& psi, FermionField& chi); \
|
||||
template void DomainWallEOFAFermion<A>::MooeeInvDag(const FermionField& psi, FermionField& chi);
|
||||
|
||||
|
||||
@ -41,7 +41,7 @@ NAMESPACE_BEGIN(Grid);
|
||||
// Pplus backwards..
|
||||
template<class Impl>
|
||||
void DomainWallEOFAFermion<Impl>::M5D(const FermionField& psi_i, const FermionField& phi_i,FermionField& chi_i,
|
||||
std::vector<Coeff_t>& lower, std::vector<Coeff_t>& diag, std::vector<Coeff_t>& upper)
|
||||
Vector<Coeff_t>& lower, Vector<Coeff_t>& diag, Vector<Coeff_t>& upper)
|
||||
{
|
||||
chi_i.Checkerboard() = psi_i.Checkerboard();
|
||||
int Ls = this->Ls;
|
||||
@ -81,7 +81,7 @@ void DomainWallEOFAFermion<Impl>::M5D(const FermionField& psi_i, const FermionFi
|
||||
|
||||
template<class Impl>
|
||||
void DomainWallEOFAFermion<Impl>::M5Ddag(const FermionField& psi_i, const FermionField& phi_i, FermionField& chi_i,
|
||||
std::vector<Coeff_t>& lower, std::vector<Coeff_t>& diag, std::vector<Coeff_t>& upper)
|
||||
Vector<Coeff_t>& lower, Vector<Coeff_t>& diag, Vector<Coeff_t>& upper)
|
||||
{
|
||||
chi_i.Checkerboard() = psi_i.Checkerboard();
|
||||
GridBase* grid = psi_i.Grid();
|
||||
@ -180,11 +180,11 @@ void DomainWallEOFAFermion<Impl>::MooeeInvDag(const FermionField& psi_i, Fermion
|
||||
|
||||
assert(psi.Checkerboard() == psi.Checkerboard());
|
||||
|
||||
std::vector<Coeff_t> ueec(Ls);
|
||||
std::vector<Coeff_t> deec(Ls+1);
|
||||
std::vector<Coeff_t> leec(Ls);
|
||||
std::vector<Coeff_t> ueemc(Ls);
|
||||
std::vector<Coeff_t> leemc(Ls);
|
||||
Vector<Coeff_t> ueec(Ls);
|
||||
Vector<Coeff_t> deec(Ls+1);
|
||||
Vector<Coeff_t> leec(Ls);
|
||||
Vector<Coeff_t> ueemc(Ls);
|
||||
Vector<Coeff_t> leemc(Ls);
|
||||
|
||||
for(int s=0; s<ueec.size(); s++){
|
||||
ueec[s] = conjugate(this->uee[s]);
|
||||
|
||||
@ -40,7 +40,7 @@ NAMESPACE_BEGIN(Grid);
|
||||
// Pplus backwards
|
||||
template<class Impl>
|
||||
void DomainWallEOFAFermion<Impl>::M5D(const FermionField& psi, const FermionField& phi,
|
||||
FermionField& chi, std::vector<Coeff_t>& lower, std::vector<Coeff_t>& diag, std::vector<Coeff_t>& upper)
|
||||
FermionField& chi, Vector<Coeff_t>& lower, Vector<Coeff_t>& diag, Vector<Coeff_t>& upper)
|
||||
{
|
||||
Coeff_t one(1.0);
|
||||
int Ls = this->Ls;
|
||||
@ -60,7 +60,7 @@ void DomainWallEOFAFermion<Impl>::M5D(const FermionField& psi, const FermionFiel
|
||||
|
||||
template<class Impl>
|
||||
void DomainWallEOFAFermion<Impl>::M5Ddag(const FermionField& psi, const FermionField& phi,
|
||||
FermionField& chi, std::vector<Coeff_t>& lower, std::vector<Coeff_t>& diag, std::vector<Coeff_t>& upper)
|
||||
FermionField& chi, Vector<Coeff_t>& lower, Vector<Coeff_t>& diag, Vector<Coeff_t>& upper)
|
||||
{
|
||||
Coeff_t one(1.0);
|
||||
int Ls = this->Ls;
|
||||
|
||||
@ -53,7 +53,7 @@ void DomainWallEOFAFermion<Impl>::MooeeInv(const FermionField& psi, FermionField
|
||||
|
||||
template<class Impl>
|
||||
void DomainWallEOFAFermion<Impl>::M5D(const FermionField& psi_i, const FermionField& phi_i, FermionField& chi_i,
|
||||
std::vector<Coeff_t>& lower, std::vector<Coeff_t>& diag, std::vector<Coeff_t>& upper)
|
||||
Vector<Coeff_t>& lower, Vector<Coeff_t>& diag, Vector<Coeff_t>& upper)
|
||||
{
|
||||
chi_i.Checkerboard() = psi_i.Checkerboard();
|
||||
GridBase* grid = psi_i.Grid();
|
||||
@ -201,7 +201,7 @@ void DomainWallEOFAFermion<Impl>::M5D(const FermionField& psi_i, const FermionFi
|
||||
|
||||
template<class Impl>
|
||||
void DomainWallEOFAFermion<Impl>::M5Ddag(const FermionField& psi_i, const FermionField& phi_i,FermionField& chi_i,
|
||||
std::vector<Coeff_t>& lower, std::vector<Coeff_t>& diag, std::vector<Coeff_t>& upper)
|
||||
Vector<Coeff_t>& lower, Vector<Coeff_t>& diag, Vector<Coeff_t>& upper)
|
||||
{
|
||||
chi_i.Checkerboard() = psi_i.Checkerboard();
|
||||
GridBase* grid = psi_i.Grid();
|
||||
|
||||
@ -197,9 +197,9 @@ void MobiusEOFAFermion<Impl>::M5D(const FermionField& psi, FermionField& chi)
|
||||
{
|
||||
int Ls = this->Ls;
|
||||
|
||||
std::vector<Coeff_t> diag(Ls,1.0);
|
||||
std::vector<Coeff_t> upper(Ls,-1.0); upper[Ls-1] = this->mq1;
|
||||
std::vector<Coeff_t> lower(Ls,-1.0); lower[0] = this->mq1;
|
||||
Vector<Coeff_t> diag(Ls,1.0);
|
||||
Vector<Coeff_t> upper(Ls,-1.0); upper[Ls-1] = this->mq1;
|
||||
Vector<Coeff_t> lower(Ls,-1.0); lower[0] = this->mq1;
|
||||
|
||||
// no shift term
|
||||
if(this->shift == 0.0){ this->M5D(psi, chi, chi, lower, diag, upper); }
|
||||
@ -213,9 +213,9 @@ void MobiusEOFAFermion<Impl>::M5Ddag(const FermionField& psi, FermionField& chi)
|
||||
{
|
||||
int Ls = this->Ls;
|
||||
|
||||
std::vector<Coeff_t> diag(Ls,1.0);
|
||||
std::vector<Coeff_t> upper(Ls,-1.0); upper[Ls-1] = this->mq1;
|
||||
std::vector<Coeff_t> lower(Ls,-1.0); lower[0] = this->mq1;
|
||||
Vector<Coeff_t> diag(Ls,1.0);
|
||||
Vector<Coeff_t> upper(Ls,-1.0); upper[Ls-1] = this->mq1;
|
||||
Vector<Coeff_t> lower(Ls,-1.0); lower[0] = this->mq1;
|
||||
|
||||
// no shift term
|
||||
if(this->shift == 0.0){ this->M5Ddag(psi, chi, chi, lower, diag, upper); }
|
||||
@ -231,9 +231,9 @@ void MobiusEOFAFermion<Impl>::Mooee(const FermionField& psi, FermionField& chi)
|
||||
int Ls = this->Ls;
|
||||
|
||||
// coefficients of Mooee
|
||||
std::vector<Coeff_t> diag = this->bee;
|
||||
std::vector<Coeff_t> upper(Ls);
|
||||
std::vector<Coeff_t> lower(Ls);
|
||||
Vector<Coeff_t> diag = this->bee;
|
||||
Vector<Coeff_t> upper(Ls);
|
||||
Vector<Coeff_t> lower(Ls);
|
||||
for(int s=0; s<Ls; s++){
|
||||
upper[s] = -this->cee[s];
|
||||
lower[s] = -this->cee[s];
|
||||
@ -254,9 +254,9 @@ void MobiusEOFAFermion<Impl>::MooeeDag(const FermionField& psi, FermionField& ch
|
||||
int Ls = this->Ls;
|
||||
|
||||
// coefficients of MooeeDag
|
||||
std::vector<Coeff_t> diag = this->bee;
|
||||
std::vector<Coeff_t> upper(Ls);
|
||||
std::vector<Coeff_t> lower(Ls);
|
||||
Vector<Coeff_t> diag = this->bee;
|
||||
Vector<Coeff_t> upper(Ls);
|
||||
Vector<Coeff_t> lower(Ls);
|
||||
for(int s=0; s<Ls; s++){
|
||||
if(s==0) {
|
||||
upper[s] = -this->cee[s+1];
|
||||
@ -315,10 +315,10 @@ void MobiusEOFAFermion<Impl>::SetCoefficientsPrecondShiftOps()
|
||||
// Tridiagonal solve for MooeeInvDag_shift_lc
|
||||
{
|
||||
Coeff_t m(0.0);
|
||||
std::vector<Coeff_t> d = Mooee_shift;
|
||||
std::vector<Coeff_t> u(Ls,0.0);
|
||||
std::vector<Coeff_t> y(Ls,0.0);
|
||||
std::vector<Coeff_t> q(Ls,0.0);
|
||||
Vector<Coeff_t> d = Mooee_shift;
|
||||
Vector<Coeff_t> u(Ls,0.0);
|
||||
Vector<Coeff_t> y(Ls,0.0);
|
||||
Vector<Coeff_t> q(Ls,0.0);
|
||||
if(pm == 1){ u[0] = 1.0; }
|
||||
else{ u[Ls-1] = 1.0; }
|
||||
|
||||
|
||||
@ -42,11 +42,11 @@ public:
|
||||
|
||||
public:
|
||||
// Shift operator coefficients for red-black preconditioned Mobius EOFA
|
||||
std::vector<Coeff_t> Mooee_shift;
|
||||
std::vector<Coeff_t> MooeeInv_shift_lc;
|
||||
std::vector<Coeff_t> MooeeInv_shift_norm;
|
||||
std::vector<Coeff_t> MooeeInvDag_shift_lc;
|
||||
std::vector<Coeff_t> MooeeInvDag_shift_norm;
|
||||
Vector<Coeff_t> Mooee_shift;
|
||||
Vector<Coeff_t> MooeeInv_shift_lc;
|
||||
Vector<Coeff_t> MooeeInv_shift_norm;
|
||||
Vector<Coeff_t> MooeeInvDag_shift_lc;
|
||||
Vector<Coeff_t> MooeeInvDag_shift_norm;
|
||||
|
||||
virtual void Instantiatable(void) {};
|
||||
|
||||
@ -74,18 +74,18 @@ public:
|
||||
// Instantiate different versions depending on Impl
|
||||
/////////////////////////////////////////////////////
|
||||
void M5D(const FermionField& psi, const FermionField& phi, FermionField& chi,
|
||||
std::vector<Coeff_t>& lower, std::vector<Coeff_t>& diag, std::vector<Coeff_t>& upper);
|
||||
Vector<Coeff_t>& lower, Vector<Coeff_t>& diag, Vector<Coeff_t>& upper);
|
||||
|
||||
void M5D_shift(const FermionField& psi, const FermionField& phi, FermionField& chi,
|
||||
std::vector<Coeff_t>& lower, std::vector<Coeff_t>& diag, std::vector<Coeff_t>& upper,
|
||||
std::vector<Coeff_t>& shift_coeffs);
|
||||
Vector<Coeff_t>& lower, Vector<Coeff_t>& diag, Vector<Coeff_t>& upper,
|
||||
Vector<Coeff_t>& shift_coeffs);
|
||||
|
||||
void M5Ddag(const FermionField& psi, const FermionField& phi, FermionField& chi,
|
||||
std::vector<Coeff_t>& lower, std::vector<Coeff_t>& diag, std::vector<Coeff_t>& upper);
|
||||
Vector<Coeff_t>& lower, Vector<Coeff_t>& diag, Vector<Coeff_t>& upper);
|
||||
|
||||
void M5Ddag_shift(const FermionField& psi, const FermionField& phi, FermionField& chi,
|
||||
std::vector<Coeff_t>& lower, std::vector<Coeff_t>& diag, std::vector<Coeff_t>& upper,
|
||||
std::vector<Coeff_t>& shift_coeffs);
|
||||
Vector<Coeff_t>& lower, Vector<Coeff_t>& diag, Vector<Coeff_t>& upper,
|
||||
Vector<Coeff_t>& shift_coeffs);
|
||||
|
||||
void MooeeInternal(const FermionField& in, FermionField& out, int dag, int inv);
|
||||
|
||||
@ -113,13 +113,13 @@ NAMESPACE_END(Grid);
|
||||
|
||||
#define INSTANTIATE_DPERP_MOBIUS_EOFA(A) \
|
||||
template void MobiusEOFAFermion<A>::M5D(const FermionField& psi, const FermionField& phi, FermionField& chi, \
|
||||
std::vector<Coeff_t>& lower, std::vector<Coeff_t>& diag, std::vector<Coeff_t>& upper); \
|
||||
Vector<Coeff_t>& lower, Vector<Coeff_t>& diag, Vector<Coeff_t>& upper); \
|
||||
template void MobiusEOFAFermion<A>::M5D_shift(const FermionField& psi, const FermionField& phi, FermionField& chi, \
|
||||
std::vector<Coeff_t>& lower, std::vector<Coeff_t>& diag, std::vector<Coeff_t>& upper, std::vector<Coeff_t>& shift_coeffs); \
|
||||
Vector<Coeff_t>& lower, Vector<Coeff_t>& diag, Vector<Coeff_t>& upper, Vector<Coeff_t>& shift_coeffs); \
|
||||
template void MobiusEOFAFermion<A>::M5Ddag(const FermionField& psi, const FermionField& phi, FermionField& chi, \
|
||||
std::vector<Coeff_t>& lower, std::vector<Coeff_t>& diag, std::vector<Coeff_t>& upper); \
|
||||
Vector<Coeff_t>& lower, Vector<Coeff_t>& diag, Vector<Coeff_t>& upper); \
|
||||
template void MobiusEOFAFermion<A>::M5Ddag_shift(const FermionField& psi, const FermionField& phi, FermionField& chi, \
|
||||
std::vector<Coeff_t>& lower, std::vector<Coeff_t>& diag, std::vector<Coeff_t>& upper, std::vector<Coeff_t>& shift_coeffs); \
|
||||
Vector<Coeff_t>& lower, Vector<Coeff_t>& diag, Vector<Coeff_t>& upper, Vector<Coeff_t>& shift_coeffs); \
|
||||
template void MobiusEOFAFermion<A>::MooeeInv(const FermionField& psi, FermionField& chi); \
|
||||
template void MobiusEOFAFermion<A>::MooeeInv_shift(const FermionField& psi, FermionField& chi); \
|
||||
template void MobiusEOFAFermion<A>::MooeeInvDag(const FermionField& psi, FermionField& chi); \
|
||||
|
||||
@ -37,7 +37,7 @@ NAMESPACE_BEGIN(Grid);
|
||||
|
||||
template<class Impl>
|
||||
void MobiusEOFAFermion<Impl>::M5D(const FermionField &psi_i, const FermionField &phi_i, FermionField &chi_i,
|
||||
std::vector<Coeff_t> &lower, std::vector<Coeff_t> &diag, std::vector<Coeff_t> &upper)
|
||||
Vector<Coeff_t> &lower, Vector<Coeff_t> &diag, Vector<Coeff_t> &upper)
|
||||
{
|
||||
chi_i.Checkerboard() = psi_i.Checkerboard();
|
||||
GridBase *grid = psi_i.Grid();
|
||||
@ -79,8 +79,8 @@ void MobiusEOFAFermion<Impl>::M5D(const FermionField &psi_i, const FermionField
|
||||
|
||||
template<class Impl>
|
||||
void MobiusEOFAFermion<Impl>::M5D_shift(const FermionField &psi_i, const FermionField &phi_i, FermionField &chi_i,
|
||||
std::vector<Coeff_t> &lower, std::vector<Coeff_t> &diag, std::vector<Coeff_t> &upper,
|
||||
std::vector<Coeff_t> &shift_coeffs)
|
||||
Vector<Coeff_t> &lower, Vector<Coeff_t> &diag, Vector<Coeff_t> &upper,
|
||||
Vector<Coeff_t> &shift_coeffs)
|
||||
{
|
||||
chi_i.Checkerboard() = psi_i.Checkerboard();
|
||||
GridBase *grid = psi_i.Grid();
|
||||
@ -127,7 +127,7 @@ void MobiusEOFAFermion<Impl>::M5D_shift(const FermionField &psi_i, const Fermion
|
||||
|
||||
template<class Impl>
|
||||
void MobiusEOFAFermion<Impl>::M5Ddag(const FermionField &psi_i, const FermionField &phi_i, FermionField &chi_i,
|
||||
std::vector<Coeff_t> &lower, std::vector<Coeff_t> &diag, std::vector<Coeff_t> &upper)
|
||||
Vector<Coeff_t> &lower, Vector<Coeff_t> &diag, Vector<Coeff_t> &upper)
|
||||
{
|
||||
chi_i.Checkerboard() = psi_i.Checkerboard();
|
||||
GridBase *grid = psi_i.Grid();
|
||||
@ -169,8 +169,8 @@ void MobiusEOFAFermion<Impl>::M5Ddag(const FermionField &psi_i, const FermionFie
|
||||
|
||||
template<class Impl>
|
||||
void MobiusEOFAFermion<Impl>::M5Ddag_shift(const FermionField &psi_i, const FermionField &phi_i, FermionField &chi_i,
|
||||
std::vector<Coeff_t> &lower, std::vector<Coeff_t> &diag, std::vector<Coeff_t> &upper,
|
||||
std::vector<Coeff_t> &shift_coeffs)
|
||||
Vector<Coeff_t> &lower, Vector<Coeff_t> &diag, Vector<Coeff_t> &upper,
|
||||
Vector<Coeff_t> &shift_coeffs)
|
||||
{
|
||||
chi_i.Checkerboard() = psi_i.Checkerboard();
|
||||
GridBase *grid = psi_i.Grid();
|
||||
|
||||
@ -40,7 +40,7 @@ NAMESPACE_BEGIN(Grid);
|
||||
// Pplus backwards
|
||||
template<class Impl>
|
||||
void MobiusEOFAFermion<Impl>::M5D(const FermionField& psi, const FermionField& phi,
|
||||
FermionField& chi, std::vector<Coeff_t>& lower, std::vector<Coeff_t>& diag, std::vector<Coeff_t>& upper)
|
||||
FermionField& chi, Vector<Coeff_t>& lower, Vector<Coeff_t>& diag, Vector<Coeff_t>& upper)
|
||||
{
|
||||
Coeff_t one(1.0);
|
||||
int Ls = this->Ls;
|
||||
@ -60,8 +60,8 @@ void MobiusEOFAFermion<Impl>::M5D(const FermionField& psi, const FermionField& p
|
||||
|
||||
template<class Impl>
|
||||
void MobiusEOFAFermion<Impl>::M5D_shift(const FermionField& psi, const FermionField& phi,
|
||||
FermionField& chi, std::vector<Coeff_t>& lower, std::vector<Coeff_t>& diag, std::vector<Coeff_t>& upper,
|
||||
std::vector<Coeff_t>& shift_coeffs)
|
||||
FermionField& chi, Vector<Coeff_t>& lower, Vector<Coeff_t>& diag, Vector<Coeff_t>& upper,
|
||||
Vector<Coeff_t>& shift_coeffs)
|
||||
{
|
||||
Coeff_t one(1.0);
|
||||
int Ls = this->Ls;
|
||||
@ -83,7 +83,7 @@ void MobiusEOFAFermion<Impl>::M5D_shift(const FermionField& psi, const FermionFi
|
||||
|
||||
template<class Impl>
|
||||
void MobiusEOFAFermion<Impl>::M5Ddag(const FermionField& psi, const FermionField& phi,
|
||||
FermionField& chi, std::vector<Coeff_t>& lower, std::vector<Coeff_t>& diag, std::vector<Coeff_t>& upper)
|
||||
FermionField& chi, Vector<Coeff_t>& lower, Vector<Coeff_t>& diag, Vector<Coeff_t>& upper)
|
||||
{
|
||||
Coeff_t one(1.0);
|
||||
int Ls = this->Ls;
|
||||
@ -103,8 +103,8 @@ void MobiusEOFAFermion<Impl>::M5Ddag(const FermionField& psi, const FermionField
|
||||
|
||||
template<class Impl>
|
||||
void MobiusEOFAFermion<Impl>::M5Ddag_shift(const FermionField& psi, const FermionField& phi,
|
||||
FermionField& chi, std::vector<Coeff_t>& lower, std::vector<Coeff_t>& diag, std::vector<Coeff_t>& upper,
|
||||
std::vector<Coeff_t>& shift_coeffs)
|
||||
FermionField& chi, Vector<Coeff_t>& lower, Vector<Coeff_t>& diag, Vector<Coeff_t>& upper,
|
||||
Vector<Coeff_t>& shift_coeffs)
|
||||
{
|
||||
Coeff_t one(1.0);
|
||||
int Ls = this->Ls;
|
||||
|
||||
@ -64,7 +64,7 @@ void MobiusEOFAFermion<Impl>::MooeeInvDag_shift(const FermionField& psi, Fermion
|
||||
|
||||
template<class Impl>
|
||||
void MobiusEOFAFermion<Impl>::M5D(const FermionField& psi_i, const FermionField& phi_i,FermionField& chi_i,
|
||||
std::vector<Coeff_t>& lower, std::vector<Coeff_t>& diag, std::vector<Coeff_t>& upper)
|
||||
Vector<Coeff_t>& lower, Vector<Coeff_t>& diag, Vector<Coeff_t>& upper)
|
||||
{
|
||||
chi_i.Checkerboard() = psi_i.Checkerboard();
|
||||
GridBase* grid = psi_i.Grid();
|
||||
@ -211,8 +211,8 @@ void MobiusEOFAFermion<Impl>::M5D(const FermionField& psi_i, const FermionField&
|
||||
|
||||
template<class Impl>
|
||||
void MobiusEOFAFermion<Impl>::M5D_shift(const FermionField& psi_i, const FermionField& phi_i,
|
||||
FermionField& chi_i, std::vector<Coeff_t>& lower, std::vector<Coeff_t>& diag, std::vector<Coeff_t>& upper,
|
||||
std::vector<Coeff_t>& shift_coeffs)
|
||||
FermionField& chi_i, Vector<Coeff_t>& lower, Vector<Coeff_t>& diag, Vector<Coeff_t>& upper,
|
||||
Vector<Coeff_t>& shift_coeffs)
|
||||
{
|
||||
#if 0
|
||||
auto & psi = psi_i;
|
||||
@ -397,7 +397,7 @@ void MobiusEOFAFermion<Impl>::M5D_shift(const FermionField& psi_i, const Fermion
|
||||
|
||||
template<class Impl>
|
||||
void MobiusEOFAFermion<Impl>::M5Ddag(const FermionField& psi_i, const FermionField& phi_i,FermionField& chi_i,
|
||||
std::vector<Coeff_t>& lower, std::vector<Coeff_t>& diag, std::vector<Coeff_t>& upper)
|
||||
Vector<Coeff_t>& lower, Vector<Coeff_t>& diag, Vector<Coeff_t>& upper)
|
||||
{
|
||||
chi_i.Checkerboard() = psi_i.Checkerboard();
|
||||
GridBase* grid = psi_i.Grid();
|
||||
@ -542,8 +542,8 @@ void MobiusEOFAFermion<Impl>::M5Ddag(const FermionField& psi_i, const FermionFie
|
||||
|
||||
template<class Impl>
|
||||
void MobiusEOFAFermion<Impl>::M5Ddag_shift(const FermionField& psi_i, const FermionField& phi_i, FermionField& chi_i,
|
||||
std::vector<Coeff_t>& lower, std::vector<Coeff_t>& diag, std::vector<Coeff_t>& upper,
|
||||
std::vector<Coeff_t>& shift_coeffs)
|
||||
Vector<Coeff_t>& lower, Vector<Coeff_t>& diag, Vector<Coeff_t>& upper,
|
||||
Vector<Coeff_t>& shift_coeffs)
|
||||
{
|
||||
#if 0
|
||||
auto & psi = psi_i;
|
||||
|
||||
@ -93,8 +93,8 @@ protected:
|
||||
RealD R;
|
||||
RealD amax;
|
||||
RealD scale;
|
||||
std::vector<double> p;
|
||||
std::vector<double> q;
|
||||
Vector<double> p;
|
||||
Vector<double> q;
|
||||
|
||||
};
|
||||
|
||||
|
||||
@ -36,7 +36,7 @@ template<class Matrix, class Field>
|
||||
class KappaSimilarityTransform {
|
||||
public:
|
||||
INHERIT_IMPL_TYPES(Matrix);
|
||||
std::vector<Coeff_t> kappa, kappaDag, kappaInv, kappaInvDag;
|
||||
Vector<Coeff_t> kappa, kappaDag, kappaInv, kappaInvDag;
|
||||
|
||||
KappaSimilarityTransform (Matrix &zmob) {
|
||||
for (int i=0;i<(int)zmob.bs.size();i++) {
|
||||
|
||||
@ -48,7 +48,7 @@ public:
|
||||
GridCartesian &FourDimGrid,
|
||||
GridRedBlackCartesian &FourDimRedBlackGrid,
|
||||
RealD _mass,RealD _M5,
|
||||
std::vector<ComplexD> &gamma, RealD b,RealD c,const ImplParams &p= ImplParams()) :
|
||||
Vector<ComplexD> &gamma, RealD b,RealD c,const ImplParams &p= ImplParams()) :
|
||||
|
||||
CayleyFermion5D<Impl>(_Umu,
|
||||
FiveDimGrid,
|
||||
@ -59,7 +59,7 @@ public:
|
||||
{
|
||||
// RealD eps = 1.0;
|
||||
std::cout<<GridLogMessage << "ZMobiusFermion (b="<<b<<",c="<<c<<") with Ls= "<<this->Ls<<" gamma passed in"<<std::endl;
|
||||
std::vector<Coeff_t> zgamma(this->Ls);
|
||||
Vector<Coeff_t> zgamma(this->Ls);
|
||||
for(int s=0;s<this->Ls;s++){
|
||||
zgamma[s] = gamma[s];
|
||||
}
|
||||
|
||||
@ -245,18 +245,18 @@ namespace Optimization {
|
||||
struct MultRealPart{
|
||||
accelerator_inline float4 operator()(float4 a, float4 b){
|
||||
float4 ymm0;
|
||||
ymm0.x = a.y;
|
||||
ymm0.y = a.y;
|
||||
ymm0.z = a.w;
|
||||
ymm0.w = a.w;
|
||||
ymm0.x = a.x;
|
||||
ymm0.y = a.x;
|
||||
ymm0.z = a.z;
|
||||
ymm0.w = a.z;
|
||||
return ymm0*b;
|
||||
// ymm0 = _mm_shuffle_ps(a,a,_MM_SELECT_FOUR_FOUR(2,2,0,0)); // ymm0 <- ar ar,
|
||||
// return _mm_mul_ps(ymm0,b); // ymm0 <- ar bi, ar br
|
||||
}
|
||||
accelerator_inline double2 operator()(double2 a, double2 b){
|
||||
double2 ymm0;
|
||||
ymm0.x = a.y;
|
||||
ymm0.y = a.y;
|
||||
ymm0.x = a.x;
|
||||
ymm0.y = a.x;
|
||||
return ymm0*b;
|
||||
// ymm0 = _mm_shuffle_pd(a,a,0x0); // ymm0 <- ar ar, ar,ar b'00,00
|
||||
// return _mm_mul_pd(ymm0,b); // ymm0 <- ar bi, ar br
|
||||
@ -265,17 +265,17 @@ namespace Optimization {
|
||||
struct MaddRealPart{
|
||||
accelerator_inline float4 operator()(float4 a, float4 b, float4 c){
|
||||
float4 ymm0; // = _mm_shuffle_ps(a,a,_MM_SELECT_FOUR_FOUR(2,2,0,0)); // ymm0 <- ar ar,
|
||||
ymm0.x = a.y;
|
||||
ymm0.y = a.y;
|
||||
ymm0.z = a.w;
|
||||
ymm0.w = a.w;
|
||||
ymm0.x = a.x;
|
||||
ymm0.y = a.x;
|
||||
ymm0.z = a.z;
|
||||
ymm0.w = a.z;
|
||||
return c+ymm0*b;
|
||||
}
|
||||
accelerator_inline double2 operator()(double2 a, double2 b, double2 c){
|
||||
// ymm0 = _mm_shuffle_pd( a, a, 0x0 );
|
||||
double2 ymm0;
|
||||
ymm0.x = a.y;
|
||||
ymm0.y = a.y;
|
||||
ymm0.x = a.x;
|
||||
ymm0.y = a.x;
|
||||
return c+ymm0*b;
|
||||
}
|
||||
};
|
||||
|
||||
@ -76,13 +76,20 @@ int main (int argc, char ** argv)
|
||||
std::cout << GridLogMessage<< "* Benchmarking DomainWallFermionR::Dhop "<<std::endl;
|
||||
std::cout << GridLogMessage<< "*********************************************************" <<std::endl;
|
||||
|
||||
GridParallelRNG RNG5(FGrid);
|
||||
GridParallelRNG RNG5(FGrid); RNG5.SeedFixedIntegers(seeds5);
|
||||
LatticeFermion src(FGrid); random(RNG5,src);
|
||||
LatticeFermion result(FGrid);
|
||||
|
||||
DomainWallFermionR Dw(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5);
|
||||
double t0,t1;
|
||||
|
||||
|
||||
typedef typename DomainWallFermionR::Coeff_t Coeff_t;
|
||||
Vector<Coeff_t> diag = Dw.bs;
|
||||
Vector<Coeff_t> upper= Dw.cs;
|
||||
Vector<Coeff_t> lower= Dw.cs;
|
||||
upper[Ls-1]=-Dw.mass*upper[Ls-1];
|
||||
lower[0] =-Dw.mass*lower[0];
|
||||
|
||||
LatticeFermion r_eo(FGrid);
|
||||
LatticeFermion src_e (FrbGrid);
|
||||
LatticeFermion src_o (FrbGrid);
|
||||
@ -99,13 +106,13 @@ int main (int argc, char ** argv)
|
||||
r_o = Zero();
|
||||
|
||||
|
||||
#define BENCH_DW(A,in,out) \
|
||||
Dw.CayleyZeroCounters(); \
|
||||
Dw. A (in,out); \
|
||||
#define BENCH_DW(A,...) \
|
||||
Dw. A (__VA_ARGS__); \
|
||||
FGrid->Barrier(); \
|
||||
Dw.CayleyZeroCounters(); \
|
||||
t0=usecond(); \
|
||||
for(int i=0;i<ncall;i++){ \
|
||||
Dw. A (in,out); \
|
||||
Dw. A (__VA_ARGS__); \
|
||||
} \
|
||||
t1=usecond(); \
|
||||
FGrid->Barrier(); \
|
||||
@ -114,9 +121,9 @@ int main (int argc, char ** argv)
|
||||
std::cout<<GridLogMessage << "******************"<<std::endl;
|
||||
|
||||
#define BENCH_ZDW(A,in,out) \
|
||||
zDw.CayleyZeroCounters(); \
|
||||
zDw. A (in,out); \
|
||||
FGrid->Barrier(); \
|
||||
zDw.CayleyZeroCounters(); \
|
||||
t0=usecond(); \
|
||||
for(int i=0;i<ncall;i++){ \
|
||||
zDw. A (in,out); \
|
||||
@ -128,9 +135,9 @@ int main (int argc, char ** argv)
|
||||
std::cout<<GridLogMessage << "******************"<<std::endl;
|
||||
|
||||
#define BENCH_DW_SSC(A,in,out) \
|
||||
Dw.CayleyZeroCounters(); \
|
||||
Dw. A (in,out); \
|
||||
FGrid->Barrier(); \
|
||||
Dw.CayleyZeroCounters(); \
|
||||
t0=usecond(); \
|
||||
for(int i=0;i<ncall;i++){ \
|
||||
__SSC_START ; \
|
||||
@ -143,23 +150,10 @@ int main (int argc, char ** argv)
|
||||
std::cout<<GridLogMessage << "Called " #A " "<< (t1-t0)/ncall<<" us"<<std::endl;\
|
||||
std::cout<<GridLogMessage << "******************"<<std::endl;
|
||||
|
||||
#define BENCH_DW_MEO(A,in,out) \
|
||||
Dw.CayleyZeroCounters(); \
|
||||
Dw. A (in,out,0); \
|
||||
FGrid->Barrier(); \
|
||||
t0=usecond(); \
|
||||
for(int i=0;i<ncall;i++){ \
|
||||
Dw. A (in,out,0); \
|
||||
} \
|
||||
t1=usecond(); \
|
||||
FGrid->Barrier(); \
|
||||
Dw.CayleyReport(); \
|
||||
std::cout<<GridLogMessage << "Called " #A " "<< (t1-t0)/ncall<<" us"<<std::endl;\
|
||||
std::cout<<GridLogMessage << "******************"<<std::endl;
|
||||
|
||||
BENCH_DW_MEO(Dhop ,src,result);
|
||||
BENCH_DW_MEO(DhopEO ,src_o,r_e);
|
||||
BENCH_DW(Dhop ,src,result,0);
|
||||
BENCH_DW(DhopEO ,src_o,r_e,0);
|
||||
BENCH_DW(Meooe ,src_o,r_e);
|
||||
BENCH_DW(M5D ,src_o,src_o,r_e,lower,diag,upper);
|
||||
BENCH_DW(Mooee ,src_o,r_o);
|
||||
BENCH_DW(MooeeInv,src_o,r_o);
|
||||
|
||||
@ -173,7 +167,7 @@ int main (int argc, char ** argv)
|
||||
std::cout << GridLogMessage<< "* Benchmarking DomainWallFermionVec5dR::Dhop "<<std::endl;
|
||||
std::cout << GridLogMessage<< "*********************************************************" <<std::endl;
|
||||
|
||||
GridParallelRNG RNG5(sFGrid);
|
||||
GridParallelRNG RNG5(sFGrid); RNG5.SeedFixedIntegers(seeds5);
|
||||
LatticeFermion src(sFGrid); random(RNG5,src);
|
||||
LatticeFermion sref(sFGrid);
|
||||
LatticeFermion result(sFGrid);
|
||||
@ -184,7 +178,7 @@ int main (int argc, char ** argv)
|
||||
|
||||
RealD b=1.5;// Scale factor b+c=2, b-c=1
|
||||
RealD c=0.5;
|
||||
std::vector<ComplexD> gamma(Ls,std::complex<double>(1.0,0.0));
|
||||
Vector<ComplexD> gamma(Ls,std::complex<double>(1.0,0.0));
|
||||
ZMobiusFermionVec5dR zDw(Umu,*sFGrid,*sFrbGrid,*sUGrid,*sUrbGrid,mass,M5,gamma,b,c);
|
||||
|
||||
std::cout<<GridLogMessage << "Calling Dhop "<<std::endl;
|
||||
@ -207,8 +201,8 @@ int main (int argc, char ** argv)
|
||||
r_e = Zero();
|
||||
r_o = Zero();
|
||||
|
||||
BENCH_DW_MEO(Dhop ,src,result);
|
||||
BENCH_DW_MEO(DhopEO ,src_o,r_e);
|
||||
BENCH_DW(Dhop ,src,result,0);
|
||||
BENCH_DW(DhopEO ,src_o,r_e,0);
|
||||
BENCH_DW_SSC(Meooe ,src_o,r_e);
|
||||
BENCH_DW(Mooee ,src_o,r_o);
|
||||
BENCH_DW(MooeeInv,src_o,r_o);
|
||||
|
||||
@ -474,7 +474,7 @@ esac
|
||||
case ${ac_COMMS} in
|
||||
*-auto)
|
||||
LX_FIND_MPI
|
||||
# if test "x$have_CXX_mpi" = 'xno'; then AC_MSG_ERROR(["The configure could not find the MPI compilation flags. N.B. The -auto mode is not supported by Cray wrappers. Use the non -auto version in this case."]); fi
|
||||
## if test "x$have_CXX_mpi" = 'xno'; then AC_MSG_ERROR(["The configure could not find the MPI compilation flags. N.B. The -auto mode is not supported by Cray wrappers. Use the non -auto version in this case."]); fi
|
||||
AM_CXXFLAGS="$MPI_CXXFLAGS $AM_CXXFLAGS"
|
||||
AM_CFLAGS="$MPI_CFLAGS $AM_CFLAGS"
|
||||
AM_LDFLAGS="`echo $MPI_CXXLDFLAGS | sed -E 's/-l@<:@^ @:>@+//g'` $AM_LDFLAGS"
|
||||
|
||||
Loading…
Reference in New Issue
Block a user