From b47f73c22231f5e7b4356bd43cee2dfdc979b3ee Mon Sep 17 00:00:00 2001 From: Peter Boyle Date: Tue, 4 Jun 2019 21:30:39 +0100 Subject: [PATCH] GPU happy --- Grid/qcd/action/fermion/WilsonImpl.h | 254 +++++++++++++++++++++++++++ 1 file changed, 254 insertions(+) create mode 100644 Grid/qcd/action/fermion/WilsonImpl.h diff --git a/Grid/qcd/action/fermion/WilsonImpl.h b/Grid/qcd/action/fermion/WilsonImpl.h new file mode 100644 index 00000000..c40bdb2b --- /dev/null +++ b/Grid/qcd/action/fermion/WilsonImpl.h @@ -0,0 +1,254 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: ./lib/qcd/action/fermion/FermionOperatorImpl.h + +Copyright (C) 2015 + +Author: Peter Boyle + +This program is free software; you can redistribute it and/or modify +it under the terms of the GNU General Public License as published by +the Free Software Foundation; either version 2 of the License, or +(at your option) any later version. + +This program is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +You should have received a copy of the GNU General Public License along +with this program; if not, write to the Free Software Foundation, Inc., +51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. + +See the full license in the file "LICENSE" in the top level distribution +directory +*************************************************************************************/ + /* END LEGAL */ +#pragma once + +NAMESPACE_BEGIN(Grid); + + +///////////////////////////////////////////////////////////////////////////// +// Single flavour four spinors with colour index +///////////////////////////////////////////////////////////////////////////// +template +class WilsonImpl : public PeriodicGaugeImpl > { +public: + + static const int Dimension = Representation::Dimension; + static const bool isFundamental = Representation::isFundamental; + static const bool LsVectorised=false; + static const int Nhcs = Options::Nhcs; + + typedef PeriodicGaugeImpl > Gimpl; + INHERIT_GIMPL_TYPES(Gimpl); + + //Necessary? + constexpr bool is_fundamental() const{return Dimension == Nc ? 1 : 0;} + + typedef typename Options::_Coeff_t Coeff_t; + typedef typename Options::template PrecisionMapper::LowerPrecVector SimdL; + + template using iImplSpinor = iScalar, Ns> >; + template using iImplPropagator = iScalar, Ns> >; + template using iImplHalfSpinor = iScalar, Nhs> >; + template using iImplHalfCommSpinor = iScalar, Nhcs> >; + template using iImplDoubledGaugeField = iVector >, Nds>; + + typedef iImplSpinor SiteSpinor; + typedef iImplPropagator SitePropagator; + typedef iImplHalfSpinor SiteHalfSpinor; + typedef iImplHalfCommSpinor SiteHalfCommSpinor; + typedef iImplDoubledGaugeField SiteDoubledGaugeField; + + typedef Lattice FermionField; + typedef Lattice PropagatorField; + typedef Lattice DoubledGaugeField; + + typedef WilsonCompressor Compressor; + typedef WilsonImplParams ImplParams; + typedef WilsonStencil StencilImpl; + typedef typename StencilImpl::View_type StencilView; + + ImplParams Params; + + WilsonImpl(const ImplParams &p = ImplParams()) : Params(p){ + assert(Params.boundary_phases.size() == Nd); + }; + + template + static accelerator_inline void multLink(_Spinor &phi, + const SiteDoubledGaugeField &U, + const _Spinor &chi, + int mu, + StencilEntry *SE, + StencilView &St) + { + auto UU = coalescedRead(U(mu)); + mult(&phi(), &UU, &chi()); + } + +#ifdef GPU_VEC + static accelerator_inline void copyLinkGpu(int lane, + SiteDoubledGaugeField & UU, + const SiteDoubledGaugeField &U) + { + auto U_l = extractLane(lane,U); + insertLane(lane,UU,U_l); + } + static accelerator_inline void multLinkGpu(int lane, + typename SiteHalfSpinor::scalar_object &phi, + const SiteDoubledGaugeField &U, + const typename SiteHalfSpinor::scalar_object &chi, + int mu) + { + auto U_l = extractLane(lane,U(mu)); + phi() = U_l * chi(); + } +#else + static accelerator_inline void multLinkGpu(int lane, + SiteHalfSpinor &phi, + const SiteDoubledGaugeField &U, + const SiteHalfSpinor &chi, + int mu) + { + auto U_l = U(mu); + phi() = U_l * chi(); + } +#endif + + static accelerator_inline void multLinkProp(SitePropagator &phi, + const SiteDoubledGaugeField &U, + const SitePropagator &chi, + int mu) + { + mult(&phi(), &U(mu), &chi()); + } + + template + static accelerator_inline void loadLinkElement(Simd ®, ref &memory) + { + reg = memory; + } + + inline void DoubleStore(GridBase *GaugeGrid, + DoubledGaugeField &Uds, + const GaugeField &Umu) + { + typedef typename Simd::scalar_type scalar_type; + + conformable(Uds.Grid(), GaugeGrid); + conformable(Umu.Grid(), GaugeGrid); + + GaugeLinkField U(GaugeGrid); + GaugeLinkField tmp(GaugeGrid); + + Lattice > coor(GaugeGrid); + //////////////////////////////////////////////////// + // apply any boundary phase or twists + //////////////////////////////////////////////////// + for (int mu = 0; mu < Nd; mu++) { + + ////////// boundary phase ///////////// + auto pha = Params.boundary_phases[mu]; + scalar_type phase( real(pha),imag(pha) ); + + int L = GaugeGrid->GlobalDimensions()[mu]; + int Lmu = L - 1; + + LatticeCoordinate(coor, mu); + + U = PeekIndex(Umu, mu); + + // apply any twists + RealD theta = Params.twist_n_2pi_L[mu] * 2*M_PI / L; + if ( theta != 0.0) { + scalar_type twphase(::cos(theta),::sin(theta)); + U = twphase*U; + std::cout << GridLogMessage << " Twist ["<(Uds, tmp, mu); + + U = adj(Cshift(U, mu, -1)); + U = where(coor == 0, conjugate(phase) * U, U); + PokeIndex(Uds, U, mu + 4); + } + } + + inline void InsertForce4D(GaugeField &mat, FermionField &Btilde, FermionField &A,int mu){ + GaugeLinkField link(mat.Grid()); + link = TraceIndex(outerProduct(Btilde,A)); + PokeIndex(mat,link,mu); + } + + inline void outerProductImpl(PropagatorField &mat, const FermionField &B, const FermionField &A){ + mat = outerProduct(B,A); + } + + inline void TraceSpinImpl(GaugeLinkField &mat, PropagatorField&P) { + mat = TraceIndex(P); + } + + inline void extractLinkField(std::vector &mat, DoubledGaugeField &Uds){ + for (int mu = 0; mu < Nd; mu++) + mat[mu] = PeekIndex(Uds, mu); + } + + + inline void InsertForce5D(GaugeField &mat, FermionField &Btilde, FermionField Ã,int mu){ + + int Ls=Btilde.Grid()->_fdimensions[0]; + GaugeLinkField tmp(mat.Grid()); + tmp = Zero(); + auto tmp_v = tmp.View(); + auto Btilde_v = Btilde.View(); + auto Atilde_v = Atilde.View(); + thread_loop( (int sss=0;sssoSites();sss++),{ + int sU=sss; + for(int s=0;s(outerProduct(Btilde_v[sF],Atilde_v[sF])); // ordering here + } + }); + PokeIndex(mat,tmp,mu); + + } +}; + + +typedef WilsonImpl WilsonImplR; // Real.. whichever prec +typedef WilsonImpl WilsonImplF; // Float +typedef WilsonImpl WilsonImplD; // Double + +typedef WilsonImpl WilsonImplRL; // Real.. whichever prec +typedef WilsonImpl WilsonImplFH; // Float +typedef WilsonImpl WilsonImplDF; // Double + +typedef WilsonImpl ZWilsonImplR; // Real.. whichever prec +typedef WilsonImpl ZWilsonImplF; // Float +typedef WilsonImpl ZWilsonImplD; // Double + +typedef WilsonImpl ZWilsonImplRL; // Real.. whichever prec +typedef WilsonImpl ZWilsonImplFH; // Float +typedef WilsonImpl ZWilsonImplDF; // Double + +typedef WilsonImpl WilsonAdjImplR; // Real.. whichever prec +typedef WilsonImpl WilsonAdjImplF; // Float +typedef WilsonImpl WilsonAdjImplD; // Double + +typedef WilsonImpl WilsonTwoIndexSymmetricImplR; // Real.. whichever prec +typedef WilsonImpl WilsonTwoIndexSymmetricImplF; // Float +typedef WilsonImpl WilsonTwoIndexSymmetricImplD; // Double + +typedef WilsonImpl WilsonTwoIndexAntiSymmetricImplR; // Real.. whichever prec +typedef WilsonImpl WilsonTwoIndexAntiSymmetricImplF; // Float +typedef WilsonImpl WilsonTwoIndexAntiSymmetricImplD; // Double + + +NAMESPACE_END(Grid); +