Skip to content

Commit

Permalink
WIP: very inefficient implementation of two-flavour tmclover together
Browse files Browse the repository at this point in the history 
with some performance improvements to single-flavour tmclover through
the usage of the CloverBlock term in AChiMinusBDPsi.

The two-flavour tmclover operator does not yet work correctly for
non-zero mass-splittings. It seems that the twisted mass part
works correctly.
  • Loading branch information
@kostrzewa
kostrzewa committed Aug 2, 2017
1 parent 24aafbc commit 67b9056
Show file tree
Hide file tree
Showing 9 changed files with 215 additions and 279 deletions.
56 changes: 26 additions & 30 deletions codegen/dslash_common.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1589,21 +1589,26 @@ void achiResult(InstVector &ivector,
}
}

// Apply clover term, and store result in out spinor.
// This is only on the AChi - bDPsi op (achimbdpsi = true)
// This is only in body kernel (face = false)
if (twisted_mass == TwistedMassVariant::none) {
clover_term(ivector, chi_spinor, false);
} else if (twisted_mass == TwistedMassVariant::degenerate) {
full_clover_term(ivector, chi_spinor, false);
} else if (twisted_mass == TwistedMassVariant::non_degenerate) {
// TODO Here something new for the ND case has to be
// implemented.
// Currently this is just copied from the degenerate case.
full_clover_term(ivector, chi_spinor, false);
} else {
unsupported_twisted_mass_variant();
}
clover_term(ivector, chi_spinor, false);

// BaKo WIP: for AChi - bDPsi, we always use CloverBlock, the twisted mass
// (or epsilon) is added using linear algebra in the EO operator
//
//// Apply clover term, and store result in out spinor.
//// This is only on the AChi - bDPsi op (achimbdpsi = true)
//// This is only in body kernel (face = false)
//if (twisted_mass == TwistedMassVariant::none) {
// clover_term(ivector, chi_spinor, false);
//} else if (twisted_mass == TwistedMassVariant::degenerate) {
// full_clover_term(ivector, chi_spinor, false);
//} else if (twisted_mass == TwistedMassVariant::non_degenerate) {
// // TODO Here something new for the ND case has to be
// // implemented.
// // Currently this is just copied from the degenerate case.
// full_clover_term(ivector, chi_spinor, false);
//} else {
// unsupported_twisted_mass_variant();
//}
} else { // NO CLOVER
if (twisted_mass == TwistedMassVariant::none) {
for (int col = 0; col < 3; col++) {
Expand Down Expand Up @@ -1641,9 +1646,8 @@ void achiResult(InstVector &ivector,
// for pure twisted mass and stores result in out_spinor
twisted_term(ivector, isPlus);
} else if (twisted_mass == TwistedMassVariant::non_degenerate) {
// TODO Here something new for the ND case has to be
// implemented.
// Currently this is just copied from the degenerate case.
// the full ND twisted mass term with mass splitting is implemented
// as a linalg functor in QPhiX, this is basically never used

// Loads spinor elements from chiBase, multiplies with A
// for pure twisted mass and stores result in out_spinor
Expand Down Expand Up @@ -1853,18 +1857,10 @@ void dslash_achimbdpsi_body(InstVector &ivector,
}

if (clover) {
if (twisted_mass == TwistedMassVariant::none) {
declare_clover(ivector);
} else if (twisted_mass == TwistedMassVariant::degenerate) {
declare_full_clover(ivector);
} else if (twisted_mass == TwistedMassVariant::non_degenerate) {
// TODO Here something new for the ND case has to be
// implemented.
// Currently this is just copied from the degenerate case.
declare_full_clover(ivector);
} else {
unsupported_twisted_mass_variant();
}
// whether we are using twisted mass or not, the clover term in AChiMinusBDPsi
// is always of the Wilson form. The twisted mass is added via the preconditioning
// mass below.
declare_clover(ivector);
}

// Fill result with a*chi
Expand Down
2 changes: 1 addition & 1 deletion codegen/jinja/tm_clov_dslash_decl.h.j2
View file
Original file line number Diff line number Diff line change
Expand Up @@ -62,7 +62,7 @@ void tm_clov_dslash_achimbdpsi_{{ isign }}_vec(
typename Types<FT, veclen, soalen, compress12>::FourSpinorBlock *oBase,
const typename Types<FT, veclen, soalen, compress12>::SU3MatrixBlock
*gBase,
const typename Types<FT, veclen, soalen, compress12>::FullCloverBlock
const typename Types<FT, veclen, soalen, compress12>::CloverBlock
*clBase,
const int xbOffs[veclen],
const int xfOffs[veclen],
Expand Down
2 changes: 1 addition & 1 deletion codegen/jinja/tm_clov_dslash_tbc_spec.h.j2
View file
Original file line number Diff line number Diff line change
Expand Up @@ -42,7 +42,7 @@ void tm_clov_dslash_achimbdpsi_{{ isign }}_vec<FPTYPE, VEC, SOA, COMPRESS12, {{
const Types<FPTYPE, VEC, SOA, COMPRESS12>::FourSpinorBlock *chiBase,
Types<FPTYPE, VEC, SOA, COMPRESS12>::FourSpinorBlock *oBase,
const Types<FPTYPE, VEC, SOA, COMPRESS12>::SU3MatrixBlock *gBase,
const Types<FPTYPE, VEC, SOA, COMPRESS12>::FullCloverBlock *clBase,
const Types<FPTYPE, VEC, SOA, COMPRESS12>::CloverBlock *clBase,
const int xbOffs[VEC],
const int xfOffs[VEC],
const int ybOffs[VEC],
Expand Down
2 changes: 1 addition & 1 deletion include/qphix/CMakeLists.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -42,7 +42,7 @@ if( ${twisted_mass} )
endif()

if( ${tm_clover} )
list(APPEND header_files twisted_mass_clover.h tm_clov_dslash_def.h tm_clov_dslash_body.h tm_clov_dslash_face.h)
list(APPEND header_files ndtm_reuse_operator_clover.h twisted_mass_clover.h tm_clov_dslash_def.h tm_clov_dslash_body.h tm_clov_dslash_face.h)
endif()

list(APPEND header_files ${PROJECT_BINARY_DIR}/include/qphix/qphix_config_internal.h)
Expand Down
2 changes: 0 additions & 2 deletions include/qphix/ndtm_reuse_operator.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -14,8 +14,6 @@ class EvenOddNDTMWilsonReuseOperator
: public TwoFlavEvenOddLinearOperator<FT, veclen, soalen, compress12>
{
public:
typedef
typename Geometry<FT, veclen, soalen, compress12>::FullCloverBlock FullCloverBlock;
typedef
typename Geometry<FT, veclen, soalen, compress12>::FourSpinorBlock FourSpinorBlock;
typedef
Expand Down
117 changes: 68 additions & 49 deletions include/qphix/ndtm_reuse_operator_clover.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,6 +2,7 @@

#include "qphix/linearOp.h"
#include "qphix/tm_clov_dslash_def.h"
#include "qphix/clover_dslash_def.h"

#include <memory>

Expand All @@ -26,108 +27,126 @@ class EvenOddNDTMCloverReuseOperator
public:
typedef typename Geometry<FT, veclen, soalen, compress12>::FullCloverBlock
FullCloverBlock;
typedef typename Geometry<FT, veclen, soalen, compress12>::CloverBlock
CloverBlock;
typedef typename Geometry<FT, veclen, soalen, compress12>::FourSpinorBlock
FourSpinorBlock;
typedef typename Geometry<FT, veclen, soalen, compress12>::SU3MatrixBlock
SU3MatrixBlock;

/**
\param[in] u Gauge fields, one for each checkerboard index.
\param[in] clov Two flavor parts of the odd-odd term. The first index is
the flavor index, the second index is for the hermitian conjugation, just
like \p clov of \ref dslashAChiMinusBDPsi or \p invclov of \ref
two_flav_dslash.
\param[in] invclov The four blocks of \f$ (A^{-1})_{ff'} \f$. The first
(slowest) index iterates through the flavor index \f$ f \f$, it coincides
with the resulting spinor flavor index. The second index is \f$ f' \f$
which coincides with the input spinor flavor indices. The last (fastest)
index is the same as for the \ref dslash function: It is the hermitian
conjugate of the inverse clover block.
\param[in] u_ Gauge fields, one for each checkerboard index.
\param[in] clov_ Single Wilson-clover style clover term which will be applied
in two_flav_AChiMinusBDPsi to all flavour components with the twisted quark
mass and the mass splitting applied subsequently.
\param[in] invclov_odiag_ The flavour off-diagonal inverse clover term
including the mass splitting.
\param[in] invclov_ FullCloverBlock inverse clover term for the flavour-diagonal
contribution as applied by the Twisted Clover Dslash.
*/
EvenOddNDTMCloverReuseOperator(SU3MatrixBlock *u_[2],
FullCloverBlock *clov_[2][2],
FullCloverBlock *invclov_[2][2][2],
double epsilon,
EvenOddNDTMCloverReuseOperator(double TwistedMass_,
double Epsilon_,
SU3MatrixBlock *u_[2],
CloverBlock *clov_,
CloverBlock *invclov_odiag_,
FullCloverBlock *invclov_[2],
Geometry<FT, veclen, soalen, compress12> *geom_,
double t_boundary,
double aniso_coeff_s,
double aniso_coeff_t,
bool use_tbc_[4] = nullptr,
double tbc_phases_[4][2] = nullptr)
: EvenOddNDTMCloverReuseOperator(
TwistedMass_, Epsilon_,
geom_, t_boundary, aniso_coeff_s, aniso_coeff_t, use_tbc_, tbc_phases_)
{
this->epsilon = epsilon;
setFields(u_, clov_, invclov_);
setFields(u_, clov_, invclov_odiag_, invclov_);
}

EvenOddNDTMCloverReuseOperator(Geometry<FT, veclen, soalen, compress12> *geom_,
EvenOddNDTMCloverReuseOperator(double TwistedMass_,
double Epsilon_,
Geometry<FT, veclen, soalen, compress12> *geom_,
double t_boundary,
double aniso_coeff_s,
double aniso_coeff_t,
bool use_tbc_[4] = nullptr,
double tbc_phases_[4][2] = nullptr)
: D(new TMClovDslash<FT, veclen, soalen, compress12>(
geom_, t_boundary, aniso_coeff_s, aniso_coeff_t, use_tbc_, tbc_phases_)),
tmp{D->getGeometry().allocCBFourSpinor(),
D->getGeometry().allocCBFourSpinor()}
: Dtmcl(new TMClovDslash<FT, veclen, soalen, compress12>(
geom_, t_boundary, aniso_coeff_s, aniso_coeff_t, use_tbc_, tbc_phases_,
TwistedMass_)),
Dwcl(new ClovDslash<FT, veclen, soalen, compress12>(
geom_, t_boundary, aniso_coeff_s, aniso_coeff_t, use_tbc_, tbc_phases_)),
tmp{Dtmcl->getGeometry().allocCBFourSpinor(),
Dtmcl->getGeometry().allocCBFourSpinor()},
odiag_tmp{Dtmcl->getGeometry().allocCBFourSpinor(),
Dtmcl->getGeometry().allocCBFourSpinor()},
Epsilon(Epsilon_)
{
}

~EvenOddNDTMCloverReuseOperator()
{
Geometry<FT, veclen, soalen, compress12> &geom = D->getGeometry();
for (int f : {0, 1}) {
geom.free(tmp[f]);
Geometry<FT, veclen, soalen, compress12> &geom = Dtmcl->getGeometry();
for (int fl : {0, 1}) {
geom.free(tmp[fl]);
geom.free(odiag_tmp[fl]);
}
}

void setFields(SU3MatrixBlock *u_[2],
FullCloverBlock *clov_[2][2],
FullCloverBlock *invclov_[2][2][2])
CloverBlock *clov_,
CloverBlock *invclov_odiag_,
FullCloverBlock *invclov_[2])
{
for (int pm : {0, 1}) {
u[pm] = u_[pm];
for (int f1 : {0, 1}) {
clov[pm][f1] = clov_[pm][f1];
for (int f2 : {0, 1}) {
invclov[pm][f1][f2] = invclov_[pm][f1][f2];
}
}
for (int cb : {0, 1}) {
u[cb] = u_[cb];
}
clov = clov_;
invclov_odiag = invclov_odiag_;
for (int fl : {0, 1}) {
invclov[fl] = invclov_[fl];
}
}

void operator()(FourSpinorBlock *res[2],
void operator()(FourSpinorBlock *const res[2],
const FourSpinorBlock *const in[2],
int isign, int target_cb) override
int isign,
int target_cb = 1) override
{
double beta = 0.25;
const int other_cb = 1-target_cb;

D->two_flav_dslash(tmp, in, u[other_cb], invclov, isign, other_cb);
D->two_flav_achimbdpsi(res, tmp, in, u[target_cb], clov, beta, epsilon, isign, target_cb);

for(int fl : {0, 1}){
Dtmcl->dslash(tmp[fl], in[fl], u[other_cb], invclov, isign, other_cb, fl);
Dwcl->dslash(odiag_tmp[fl], in[1-fl], u[other_cb], invclov_odiag, isign, other_cb);
}
axpy<FT, veclen, soalen, compress12, 2>(1.0, odiag_tmp, tmp, getGeometry(), getGeometry().getNSIMT());

Dtmcl->two_flav_AChiMinusBDPsi(res, tmp, in, u[target_cb], clov, beta, Epsilon, isign, target_cb);
}

Geometry<FT, veclen, soalen, compress12> &getGeometry()
{
return D->getGeometry();
return Dtmcl->getGeometry();
}

private:
double Mass;
double epsilon;
double Epsilon;

// XXX Use a simple member, then there is no need for a virtual function
// call within `operator()`.
std::unique_ptr<TMClovDslash<FT, veclen, soalen, compress12>> D;
std::unique_ptr<TMClovDslash<FT, veclen, soalen, compress12>> Dtmcl;
std::unique_ptr<ClovDslash<FT, veclen, soalen, compress12>> Dwcl;

SU3MatrixBlock const *u[2];
FullCloverBlock const *clov[2][2];
FullCloverBlock const *invclov[2][2][2];
CloverBlock const *clov;
CloverBlock const *invclov_odiag;
FullCloverBlock const *invclov[2];

FourSpinorBlock *tmp[2];
FourSpinorBlock *odiag_tmp[2];

}; // Class
}; // Namespace
Loading

0 comments on commit 67b9056

Please sign in to comment.