Skip to content

Commit

Permalink
minor cleaning of the contraction.cuh, helper structures are now move…
Browse files Browse the repository at this point in the history
…d into gamma.cuh
  • Loading branch information
alexstrel committed Dec 4, 2023
1 parent f49efe2 commit 4d76007
Show file tree
Hide file tree
Showing 2 changed files with 144 additions and 202 deletions.
139 changes: 139 additions & 0 deletions include/gamma.cuh
Original file line number Diff line number Diff line change
Expand Up @@ -290,5 +290,144 @@ namespace quda {
//Returns the type of Gamma matrix
inline constexpr int Dir() const { return dir; }
};

// list of specialized structures used in the contraction kernels:
static constexpr int nspin = 4;

constexpr array<array<int, nspin>, nspin*nspin> get_dr_gm_i() {
return {{
// VECTORS
// G_idx = 1: \gamma_1
{3, 2, 1, 0},

// G_idx = 2: \gamma_2
{3, 2, 1, 0},

// G_idx = 3: \gamma_3
{2, 3, 0, 1},

// G_idx = 4: \gamma_4
{2, 3, 0, 1},

// PSEUDO-VECTORS
// G_idx = 6: \gamma_5\gamma_1
{3, 2, 1, 0},

// G_idx = 7: \gamma_5\gamma_2
{3, 2, 1, 0},

// G_idx = 8: \gamma_5\gamma_3
{2, 3, 0, 1},

// G_idx = 9: \gamma_5\gamma_4
{2, 3, 0, 1},

// SCALAR
// G_idx = 0: I
{0, 1, 2, 3},

// PSEUDO-SCALAR
// G_idx = 5: \gamma_5
{0, 1, 2, 3},

// TENSORS
// G_idx = 10: (i/2) * [\gamma_1, \gamma_2]
{0, 1, 2, 3},

// G_idx = 11: (i/2) * [\gamma_1, \gamma_3]. this matrix was corrected
{1, 0, 3, 2},

// G_idx = 12: (i/2) * [\gamma_1, \gamma_4]
{1, 0, 3, 2},

// G_idx = 13: (i/2) * [\gamma_2, \gamma_3]
{1, 0, 3, 2},

// G_idx = 14: (i/2) * [\gamma_2, \gamma_4]
{1, 0, 3, 2},

// G_idx = 15: (i/2) * [\gamma_3, \gamma_4]. this matrix was corrected
{0, 1, 2, 3}
}};
}

template<typename T>
constexpr array<array<complex<T>, nspin>, nspin*nspin> get_dr_g5gm_z() {

constexpr complex<T> i = complex<T>(0., 1.);
constexpr complex<T> one = complex<T>(1., 0.);

return {{
// VECTORS
// G_idx = 1: \gamma_1
{i, i, i, i},

// G_idx = 2: \gamma_2
{-one, one, -one, one},

// G_idx = 3: \gamma_3
{i, -i, i, -i},

// G_idx = 4: \gamma_4
{one, one, -one, -one},

// PSEUDO-VECTORS
// G_idx = 6: \gamma_5\gamma_1
{i, i, -i, -i},

// G_idx = 7: \gamma_5\gamma_2
{-one, one, one, -one},

// G_idx = 8: \gamma_5\gamma_3
{i, -i, -i, i},

// G_idx = 9: \gamma_5\gamma_4
{one, one, one, one},

// SCALAR
// G_idx = 0: I
{one, one, -one, -one},

// PSEUDO-SCALAR
// G_idx = 5: \gamma_5
{one, one, one, one},

// TENSORS
// G_idx = 10: (i/2) * [\gamma_1, \gamma_2]
{one, -one, -one, one},

// G_idx = 11: (i/2) * [\gamma_1, \gamma_3]. this matrix was corrected
{-i, i, i, -i},

// G_idx = 12: (i/2) * [\gamma_1, \gamma_4]
{-one, -one, -one, -one},

// G_idx = 13: (i/2) * [\gamma_2, \gamma_3]
{one, one, -one, -one},

// G_idx = 14: (i/2) * [\gamma_2, \gamma_4]
{-i, i, -i, i},

// G_idx = 15: (i/2) * [\gamma_3, \gamma_4]. this matrix was corrected
{-one, one, -one, one}
}};
}

template <typename real, int nSpin> struct DRGamma { };

template <> struct DRGamma<double, 4> {
static constexpr int nSpin = 4;
//
const array<array<int, nSpin>, nSpin*nSpin> gm_i = get_dr_gm_i();
const array<array<complex<double>, nSpin>, nSpin*nSpin> g5gm_z = get_dr_g5gm_z<double>();
};

template <> struct DRGamma<float, 4> {
static constexpr int nSpin = 4;
//
const array<array<int, nSpin>, nSpin*nSpin> gm_i = get_dr_gm_i();
const array<array<complex<float>, nSpin>, nSpin*nSpin> g5gm_z = get_dr_g5gm_z<float>();
};


} // namespace quda
207 changes: 5 additions & 202 deletions include/kernels/contraction.cuh
Original file line number Diff line number Diff line change
Expand Up @@ -5,214 +5,15 @@
#include <quda_matrix.h>
#include <matrix_field.h>
#include <kernel.h>
#include <gamma.cuh>

namespace quda
{
static constexpr int max_contract_results = 16; // sized for nSpin**2 = 16


using spinor_array = array<array<double, 2>, max_contract_results>;
using staggered_spinor_array = array<double, 2>;

template <typename real, int nSpin = 4> class DRGammaMatrix {
public:
// Stores gamma matrix column index for non-zero complex value.
// This is shared by g5gm, gmg5.
int gm_i[nSpin*nSpin][nSpin] = {};

// Stores gamma matrix non-zero complex value for the corresponding g5gm_i
complex<real> g5gm_z[nSpin*nSpin][nSpin] = {};

// use tr[Gamma*Prop*Gamma*g5*conj(Prop)*g5] = tr[g5*Gamma*Prop*g5*Gamma*(-1)^{?}*conj(Prop)].
//the possible minus sign will be taken care of in the main function
//! Constructor
DRGammaMatrix() {
//if constexpr (nSpin == 4) {
if (nSpin == 4) {
const complex<real> i(0., 1.);
// VECTORS
// G_idx = 1: \gamma_1
gm_i[0][0] = 3;
gm_i[0][1] = 2;
gm_i[0][2] = 1;
gm_i[0][3] = 0;

g5gm_z[0][0] = i;
g5gm_z[0][1] = i;
g5gm_z[0][2] = i;
g5gm_z[0][3] = i;

// G_idx = 2: \gamma_2
gm_i[1][0] = 3;
gm_i[1][1] = 2;
gm_i[1][2] = 1;
gm_i[1][3] = 0;

g5gm_z[1][0] = -1.;
g5gm_z[1][1] = 1.;
g5gm_z[1][2] = -1.;
g5gm_z[1][3] = 1.;

// G_idx = 3: \gamma_3
gm_i[2][0] = 2;
gm_i[2][1] = 3;
gm_i[2][2] = 0;
gm_i[2][3] = 1;

g5gm_z[2][0] = i;
g5gm_z[2][1] = -i;
g5gm_z[2][2] = i;
g5gm_z[2][3] = -i;

// G_idx = 4: \gamma_4
gm_i[3][0] = 2;
gm_i[3][1] = 3;
gm_i[3][2] = 0;
gm_i[3][3] = 1;

g5gm_z[3][0] = 1.;
g5gm_z[3][1] = 1.;
g5gm_z[3][2] = -1.;
g5gm_z[3][3] = -1.;


// PSEUDO-VECTORS
// G_idx = 6: \gamma_5\gamma_1
gm_i[4][0] = 3;
gm_i[4][1] = 2;
gm_i[4][2] = 1;
gm_i[4][3] = 0;

g5gm_z[4][0] = i;
g5gm_z[4][1] = i;
g5gm_z[4][2] = -i;
g5gm_z[4][3] = -i;

// G_idx = 7: \gamma_5\gamma_2
gm_i[5][0] = 3;
gm_i[5][1] = 2;
gm_i[5][2] = 1;
gm_i[5][3] = 0;

g5gm_z[5][0] = -1.;
g5gm_z[5][1] = 1.;
g5gm_z[5][2] = 1.;
g5gm_z[5][3] = -1.;

// G_idx = 8: \gamma_5\gamma_3
gm_i[6][0] = 2;
gm_i[6][1] = 3;
gm_i[6][2] = 0;
gm_i[6][3] = 1;

g5gm_z[6][0] = i;
g5gm_z[6][1] = -i;
g5gm_z[6][2] = -i;
g5gm_z[6][3] = i;

// G_idx = 9: \gamma_5\gamma_4
gm_i[7][0] = 2;
gm_i[7][1] = 3;
gm_i[7][2] = 0;
gm_i[7][3] = 1;

g5gm_z[7][0] = 1.;
g5gm_z[7][1] = 1.;
g5gm_z[7][2] = 1.;
g5gm_z[7][3] = 1.;

// SCALAR
// G_idx = 0: I
gm_i[8][0] = 0;
gm_i[8][1] = 1;
gm_i[8][2] = 2;
gm_i[8][3] = 3;

g5gm_z[8][0] = 1.;
g5gm_z[8][1] = 1.;
g5gm_z[8][2] = -1.;
g5gm_z[8][3] = -1.;


// PSEUDO-SCALAR
// G_idx = 5: \gamma_5
gm_i[9][0] = 0;
gm_i[9][1] = 1;
gm_i[9][2] = 2;
gm_i[9][3] = 3;

g5gm_z[9][0] = 1.;
g5gm_z[9][1] = 1.;
g5gm_z[9][2] = 1.;
g5gm_z[9][3] = 1.;

// TENSORS
// G_idx = 10: (i/2) * [\gamma_1, \gamma_2]
gm_i[10][0] = 0;
gm_i[10][1] = 1;
gm_i[10][2] = 2;
gm_i[10][3] = 3;

g5gm_z[10][0] = 1.;
g5gm_z[10][1] = -1.;
g5gm_z[10][2] = -1.;
g5gm_z[10][3] = 1.;

// G_idx = 11: (i/2) * [\gamma_1, \gamma_3]. this matrix was corrected
gm_i[11][0] = 1;
gm_i[11][1] = 0;
gm_i[11][2] = 3;
gm_i[11][3] = 2;

g5gm_z[11][0] = -i;
g5gm_z[11][1] = i;
g5gm_z[11][2] = i;
g5gm_z[11][3] = -i;

// G_idx = 12: (i/2) * [\gamma_1, \gamma_4]
gm_i[12][0] = 1;
gm_i[12][1] = 0;
gm_i[12][2] = 3;
gm_i[12][3] = 2;

g5gm_z[12][0] = -1.;
g5gm_z[12][1] = -1.;
g5gm_z[12][2] = -1.;
g5gm_z[12][3] = -1.;

// G_idx = 13: (i/2) * [\gamma_2, \gamma_3]
gm_i[13][0] = 1;
gm_i[13][1] = 0;
gm_i[13][2] = 3;
gm_i[13][3] = 2;

g5gm_z[13][0] = 1.;
g5gm_z[13][1] = 1.;
g5gm_z[13][2] = -1.;
g5gm_z[13][3] = -1.;
// G_idx = 14: (i/2) * [\gamma_2, \gamma_4]
gm_i[14][0] = 1;
gm_i[14][1] = 0;
gm_i[14][2] = 3;
gm_i[14][3] = 2;

g5gm_z[14][0] = -i;
g5gm_z[14][1] = i;
g5gm_z[14][2] = -i;
g5gm_z[14][3] = i;

// G_idx = 15: (i/2) * [\gamma_3, \gamma_4]. this matrix was corrected
gm_i[15][0] = 0;
gm_i[15][1] = 1;
gm_i[15][2] = 2;
gm_i[15][3] = 3;

g5gm_z[15][0] = -1.;
g5gm_z[15][1] = 1.;
g5gm_z[15][2] = -1.;
g5gm_z[15][3] = 1.;
} // end if constexpr
};
};

template <int reduction_dim, class T> __device__ void sink_from_t_xyz(int sink[4], int t, int xyz, T X[4])
{
Expand Down Expand Up @@ -256,15 +57,17 @@ namespace quda
static constexpr bool spin_project = nSpin_ == 1 ? false : true;
static constexpr bool spinor_direct_load = false; // false means texture load

const DRGamma<real, nSpin> Gamma;//empty for nSpin != 4

typedef typename colorspinor_mapper<Float, nSpin, nColor, spin_project, spinor_direct_load>::type F;

F x;
F y;
int s1, b1;
int mom_mode[4];
QudaFFTSymmType fft_type[4];
int source_position[4];
int NxNyNzNt[4];
DRGammaMatrix<real> Gamma;
int t_offset;
int offsets[4];

Expand Down

0 comments on commit 4d76007

Please sign in to comment.