haraldwolff
/
CLBlast
mirror of https://github.com/CNugteren/CLBlast.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

Separated the tuning parameters of the new direct GEMM kernel from the indirect version

pull/108/head
Cedric Nugteren 2016-09-25 13:52:08 +02:00
parent 140dc12854
commit 669f43aed6
5 changed files with 299 additions and 151 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
src/database/database.cpp
View File

@ -21,6 +21,7 @@
#include "database/kernels/xgemv_fast_rot.hpp"
#include "database/kernels/xger.hpp"
#include "database/kernels/xgemm.hpp"
#include "database/kernels/xgemm_direct.hpp"
#include "database/kernels/copy.hpp"
#include "database/kernels/pad.hpp"
#include "database/kernels/transpose.hpp"
@ -38,6 +39,7 @@ const std::vector<Database::DatabaseEntry> Database::database = {
XgemvFastRotHalf, XgemvFastRotSingle, XgemvFastRotDouble, XgemvFastRotComplexSingle, XgemvFastRotComplexDouble,
XgerHalf, XgerSingle, XgerDouble, XgerComplexSingle, XgerComplexDouble,
XgemmHalf, XgemmSingle, XgemmDouble, XgemmComplexSingle, XgemmComplexDouble,
XgemmDirectHalf, XgemmDirectSingle, XgemmDirectDouble, XgemmDirectComplexSingle, XgemmDirectComplexDouble,
CopyHalf, CopySingle, CopyDouble, CopyComplexSingle, CopyComplexDouble,
PadHalf, PadSingle, PadDouble, PadComplexSingle, PadComplexDouble,
TransposeHalf, TransposeSingle, TransposeDouble, TransposeComplexSingle, TransposeComplexDouble,

1
src/database/database.hpp
View File

@ -75,6 +75,7 @@ class Database {
static const DatabaseEntry XgemvFastRotHalf, XgemvFastRotSingle, XgemvFastRotDouble, XgemvFastRotComplexSingle, XgemvFastRotComplexDouble;
static const DatabaseEntry XgerHalf, XgerSingle, XgerDouble, XgerComplexSingle, XgerComplexDouble;
static const DatabaseEntry XgemmHalf, XgemmSingle, XgemmDouble, XgemmComplexSingle, XgemmComplexDouble;
static const DatabaseEntry XgemmDirectHalf, XgemmDirectSingle, XgemmDirectDouble, XgemmDirectComplexSingle, XgemmDirectComplexDouble;
static const DatabaseEntry CopyHalf, CopySingle, CopyDouble, CopyComplexSingle, CopyComplexDouble;
static const DatabaseEntry PadHalf, PadSingle, PadDouble, PadComplexSingle, PadComplexDouble;
static const DatabaseEntry TransposeHalf, TransposeSingle, TransposeDouble, TransposeComplexSingle, TransposeComplexDouble;

76
src/database/kernels/xgemm_direct.hpp 100644
View File

@ -0,0 +1,76 @@
// =================================================================================================
// This file is part of the CLBlast project. The project is licensed under Apache Version 2.0. This
// project loosely follows the Google C++ styleguide and uses a tab-size of two spaces and a max-
// width of 100 characters per line.
//
// Author(s):
// Database generator <database.py>
//
// This file populates the database with best-found tuning parameters for the 'Xgemm' kernels.
//
// =================================================================================================
namespace clblast {
// =================================================================================================
const Database::DatabaseEntry Database::XgemmDirectHalf = {
"XgemmDirect", Precision::kHalf, {
{ // Default
kDeviceTypeAll, "default", {
{ "default", { {"KWGD",32}, {"KWID",2}, {"MDIMAD",8}, {"MDIMCD",8}, {"MWGD",32}, {"NDIMBD",8}, {"NDIMCD",8}, {"NWGD",32}, {"VWMD",1}, {"VWND",1} } },
}
},
}
};
// =================================================================================================
const Database::DatabaseEntry Database::XgemmDirectSingle = {
"XgemmDirect", Precision::kSingle, {
{ // Default
kDeviceTypeAll, "default", {
{ "default", { {"KWGD",32}, {"KWID",2}, {"MDIMAD",8}, {"MDIMCD",8}, {"MWGD",32}, {"NDIMBD",8}, {"NDIMCD",8}, {"NWGD",32}, {"VWMD",1}, {"VWND",1} } },
}
},
}
};
// =================================================================================================
const Database::DatabaseEntry Database::XgemmDirectComplexSingle = {
"XgemmDirect", Precision::kComplexSingle, {
{ // Default
kDeviceTypeAll, "default", {
{ "default", { {"KWGD",32}, {"KWID",2}, {"MDIMAD",8}, {"MDIMCD",8}, {"MWGD",32}, {"NDIMBD",8}, {"NDIMCD",8}, {"NWGD",32}, {"VWMD",1}, {"VWND",1} } },
}
},
}
};
// =================================================================================================
const Database::DatabaseEntry Database::XgemmDirectDouble = {
"XgemmDirect", Precision::kDouble, {
{ // Default
kDeviceTypeAll, "default", {
{ "default", { {"KWGD",32}, {"KWID",2}, {"MDIMAD",8}, {"MDIMCD",8}, {"MWGD",32}, {"NDIMBD",8}, {"NDIMCD",8}, {"NWGD",32}, {"VWMD",1}, {"VWND",1} } },
}
},
}
};
// =================================================================================================
const Database::DatabaseEntry Database::XgemmDirectComplexDouble = {
"XgemmDirect", Precision::kComplexDouble, {
{ // Default
kDeviceTypeAll, "default", {
{ "default", { {"KWGD",32}, {"KWID",2}, {"MDIMAD",8}, {"MDIMCD",8}, {"MWGD",32}, {"NDIMBD",8}, {"NDIMCD",8}, {"NWGD",32}, {"VWMD",1}, {"VWND",1} } },
}
},
}
};
// =================================================================================================
} // namespace clblast

358
src/kernels/level3/xgemm_direct.opencl
View File

@ -16,68 +16,140 @@
// literal). Comment-out this line for syntax-highlighting when developing.
R"(
// Parameters set by the tuner or by the database. Here they are given a basic default value in case
// this kernel file is used outside of the CLBlast library. Note that all parameters here have a
// suffix 'D' to denote that they are for the 'direct' version of the GEMM kernel.
#ifndef MWGD
#define MWGD 8 // Tile-size in dimension M (e.g. 64, 128)
#endif
#ifndef NWGD
#define NWGD 8 // Tile-size in dimension N (e.g. 64, 128)
#endif
#ifndef KWGD
#define KWGD 8 // Tile-size in dimension K (e.g. 8, 16)
#endif
#ifndef MDIMCD
#define MDIMCD 8 // Threads per workgroup in M-dimension (e.g. 8, 16, 32)
#endif
#ifndef NDIMCD
#define NDIMCD 8 // Threads per workgroup in N-dimension (e.g. 8, 16, 32)
#endif
#ifndef MDIMAD
#define MDIMAD 8 // Re-shaped tile dimension of matrix A: KDIMAD * MDIMAD
#endif
#ifndef NDIMBD
#define NDIMBD 8 // Re-shaped tile dimension of matrix B: KDIMBD * NDIMBD
#endif
#ifndef KWID
#define KWID 1 // Unroll factor of the KWGD loop (smaller or equal than KWGD)
#endif
#ifndef VWMD
#define VWMD 1 // Vector width of matrices A and C
#endif
#ifndef VWND
#define VWND 1 // Vector width of matrix B
#endif
// Helper parameters based on the above tuning parameters
#define MWID (MWGD/MDIMCD) // Work per work-item (M-dimension)
#define NWID (NWGD/NDIMCD) // Work per work-item (N-dimension)
#define KDIMAD ((MDIMCD*NDIMCD)/(MDIMAD)) // Re-shaped tile dimension of matrix A: KDIMAD * MDIMAD
#define KDIMBD ((MDIMCD*NDIMCD)/(NDIMBD)) // Re-shaped tile dimension of matrix B: KDIMBD * NDIMBD
#define MWAD (MWGD/MDIMAD) // Amount of loads-per-thread for matrix A (M-dimension)
#define KWAD (KWGD/KDIMAD) // Amount of loads-per-thread for matrix A (K-dimension)
#define KWBD (KWGD/KDIMBD) // Amount of loads-per-thread for matrix B (K-dimension)
#define NWBD (NWGD/NDIMBD) // Amount of loads-per-thread for matrix B (N-dimension)
// =================================================================================================
// Data-widths in dimension M
#if VWMD == 1
typedef real realMD;
#elif VWMD == 2
typedef real2 realMD;
#elif VWMD == 4
typedef real4 realMD;
#elif VWMD == 8
typedef real8 realMD;
#elif VWMD == 16
typedef real16 realMD;
#endif
// Data-widths in dimension N
#if VWND == 1
typedef real realND;
#elif VWND == 2
typedef real2 realND;
#elif VWND == 4
typedef real4 realND;
#elif VWND == 8
typedef real8 realND;
#elif VWND == 16
typedef real16 realND;
#endif
// =================================================================================================
// Caches global off-chip memory into local (shared) memory on-chip. This function is specific for
// caching the A input matrix.
inline void GlobalToLocalDirectA(const __global realM* restrict agm, __local real* alm,
inline void GlobalToLocalDirectA(const __global realMD* restrict agm, __local real* alm,
const int a_ld, const int a_offset, const int tid, const int kwg,
const int a_transpose, const int a_conjugate) {
const int la0 = tid % MDIMA;
const int la1 = tid / MDIMA;
const int la0 = tid % MDIMAD;
const int la1 = tid / MDIMAD;
#pragma unroll
for (int mia=0; mia<MWA/VWM; ++mia) {
for (int mia=0; mia<MWAD/VWMD; ++mia) {
#pragma unroll
for (int kia=0; kia<KWA; ++kia) {
for (int kia=0; kia<KWAD; ++kia) {
// Computes the indices for the global memory
int mg = mia + la0*(MWA/VWM);
int kg = kia + la1*KWA;
int idm = (a_transpose) ? mg + kwg/VWM : mg + GetGroupID0()*(MWG/VWM);
int idk = (a_transpose) ? kg + GetGroupID0()*MWG : kg + kwg;
int mg = mia + la0*(MWAD/VWMD);
int kg = kia + la1*KWAD;
int idm = (a_transpose) ? mg + kwg/VWMD : mg + GetGroupID0()*(MWGD/VWMD);
int idk = (a_transpose) ? kg + GetGroupID0()*MWGD : kg + kwg;
// Loads the data from global memory into the local memory
const realM avec = agm[idk*(a_ld/VWM) + idm + a_offset];
#if VWM == 1
alm[kg*MWG + mg] = avec;
#elif VWM == 2
alm[kg*MWG + mg*VWM + 0] = avec.x;
alm[kg*MWG + mg*VWM + 1] = avec.y;
#elif VWM == 4
alm[kg*MWG + mg*VWM + 0] = avec.x;
alm[kg*MWG + mg*VWM + 1] = avec.y;
alm[kg*MWG + mg*VWM + 2] = avec.z;
alm[kg*MWG + mg*VWM + 3] = avec.w;
#elif VWM == 8
alm[kg*MWG + mg*VWM + 0] = avec.s0;
alm[kg*MWG + mg*VWM + 1] = avec.s1;
alm[kg*MWG + mg*VWM + 2] = avec.s2;
alm[kg*MWG + mg*VWM + 3] = avec.s3;
alm[kg*MWG + mg*VWM + 4] = avec.s4;
alm[kg*MWG + mg*VWM + 5] = avec.s5;
alm[kg*MWG + mg*VWM + 6] = avec.s6;
alm[kg*MWG + mg*VWM + 7] = avec.s7;
#elif VWM == 16
alm[kg*MWG + mg*VWM + 0] = avec.s0;
alm[kg*MWG + mg*VWM + 1] = avec.s1;
alm[kg*MWG + mg*VWM + 2] = avec.s2;
alm[kg*MWG + mg*VWM + 3] = avec.s3;
alm[kg*MWG + mg*VWM + 4] = avec.s4;
alm[kg*MWG + mg*VWM + 5] = avec.s5;
alm[kg*MWG + mg*VWM + 6] = avec.s6;
alm[kg*MWG + mg*VWM + 7] = avec.s7;
alm[kg*MWG + mg*VWM + 8] = avec.s8;
alm[kg*MWG + mg*VWM + 9] = avec.s9;
alm[kg*MWG + mg*VWM + 10] = avec.sA;
alm[kg*MWG + mg*VWM + 11] = avec.sB;
alm[kg*MWG + mg*VWM + 12] = avec.sC;
alm[kg*MWG + mg*VWM + 13] = avec.sD;
alm[kg*MWG + mg*VWM + 14] = avec.sE;
alm[kg*MWG + mg*VWM + 15] = avec.sF;
const realMD avec = agm[idk*(a_ld/VWMD) + idm + a_offset];
#if VWMD == 1
alm[kg*MWGD + mg] = avec;
#elif VWMD == 2
alm[kg*MWGD + mg*VWMD + 0] = avec.x;
alm[kg*MWGD + mg*VWMD + 1] = avec.y;
#elif VWMD == 4
alm[kg*MWGD + mg*VWMD + 0] = avec.x;
alm[kg*MWGD + mg*VWMD + 1] = avec.y;
alm[kg*MWGD + mg*VWMD + 2] = avec.z;
alm[kg*MWGD + mg*VWMD + 3] = avec.w;
#elif VWMD == 8
alm[kg*MWGD + mg*VWMD + 0] = avec.s0;
alm[kg*MWGD + mg*VWMD + 1] = avec.s1;
alm[kg*MWGD + mg*VWMD + 2] = avec.s2;
alm[kg*MWGD + mg*VWMD + 3] = avec.s3;
alm[kg*MWGD + mg*VWMD + 4] = avec.s4;
alm[kg*MWGD + mg*VWMD + 5] = avec.s5;
alm[kg*MWGD + mg*VWMD + 6] = avec.s6;
alm[kg*MWGD + mg*VWMD + 7] = avec.s7;
#elif VWMD == 16
alm[kg*MWGD + mg*VWMD + 0] = avec.s0;
alm[kg*MWGD + mg*VWMD + 1] = avec.s1;
alm[kg*MWGD + mg*VWMD + 2] = avec.s2;
alm[kg*MWGD + mg*VWMD + 3] = avec.s3;
alm[kg*MWGD + mg*VWMD + 4] = avec.s4;
alm[kg*MWGD + mg*VWMD + 5] = avec.s5;
alm[kg*MWGD + mg*VWMD + 6] = avec.s6;
alm[kg*MWGD + mg*VWMD + 7] = avec.s7;
alm[kg*MWGD + mg*VWMD + 8] = avec.s8;
alm[kg*MWGD + mg*VWMD + 9] = avec.s9;
alm[kg*MWGD + mg*VWMD + 10] = avec.sA;
alm[kg*MWGD + mg*VWMD + 11] = avec.sB;
alm[kg*MWGD + mg*VWMD + 12] = avec.sC;
alm[kg*MWGD + mg*VWMD + 13] = avec.sD;
alm[kg*MWGD + mg*VWMD + 14] = avec.sE;
alm[kg*MWGD + mg*VWMD + 15] = avec.sF;
#endif
if (a_conjugate) {
for (int vm=0; vm<VWM; ++vm) {
COMPLEX_CONJUGATE(alm[kg*MWG + mg*VWM + vm]);
for (int vm=0; vm<VWMD; ++vm) {
COMPLEX_CONJUGATE(alm[kg*MWGD + mg*VWMD + vm]);
}
}
}
@ -85,64 +157,64 @@ inline void GlobalToLocalDirectA(const __global realM* restrict agm, __local rea
}
// Same as above, but now for the B input matrix
inline void GlobalToLocalDirectB(const __global realN* restrict bgm, __local real* blm,
inline void GlobalToLocalDirectB(const __global realND* restrict bgm, __local real* blm,
const int b_ld, const int b_offset, const int tid, const int kwg,
const int b_transpose, const int b_conjugate) {
const int lb0 = tid % NDIMB;
const int lb1 = tid / NDIMB;
const int lb0 = tid % NDIMBD;
const int lb1 = tid / NDIMBD;
#pragma unroll
for (int kib=0; kib<KWB; ++kib) {
for (int kib=0; kib<KWBD; ++kib) {
#pragma unroll
for (int nib=0; nib<NWB/VWN; ++nib) {
for (int nib=0; nib<NWBD/VWND; ++nib) {
// Computes the indices for the global memory
int ng = nib + lb0*(NWB/VWN);
int kg = kib + lb1*KWB;
int idn = (b_transpose) ? ng + kwg/VWN : ng + GetGroupID1()*(NWG/VWN);
int idk = (b_transpose) ? kg + GetGroupID1()*NWG : kg + kwg;
int ng = nib + lb0*(NWBD/VWND);
int kg = kib + lb1*KWBD;
int idn = (b_transpose) ? ng + kwg/VWND : ng + GetGroupID1()*(NWGD/VWND);
int idk = (b_transpose) ? kg + GetGroupID1()*NWGD : kg + kwg;
// Loads the data from global memory into the local memory
const realM bvec = bgm[idk*(b_ld/VWN) + idn + b_offset];
#if VWN == 1
blm[kg*NWG + ng] = bvec;
#elif VWN == 2
blm[kg*NWG + ng*VWN + 0] = bvec.x;
blm[kg*NWG + ng*VWN + 1] = bvec.y;
#elif VWN == 4
blm[kg*NWG + ng*VWN + 0] = bvec.x;
blm[kg*NWG + ng*VWN + 1] = bvec.y;
blm[kg*NWG + ng*VWN + 2] = bvec.z;
blm[kg*NWG + ng*VWN + 3] = bvec.w;
#elif VWN == 8
blm[kg*NWG + ng*VWN + 0] = bvec.s0;
blm[kg*NWG + ng*VWN + 1] = bvec.s1;
blm[kg*NWG + ng*VWN + 2] = bvec.s2;
blm[kg*NWG + ng*VWN + 3] = bvec.s3;
blm[kg*NWG + ng*VWN + 4] = bvec.s4;
blm[kg*NWG + ng*VWN + 5] = bvec.s5;
blm[kg*NWG + ng*VWN + 6] = bvec.s6;
blm[kg*NWG + ng*VWN + 7] = bvec.s7;
#elif VWN == 16
blm[kg*NWG + ng*VWN + 0] = bvec.s0;
blm[kg*NWG + ng*VWN + 1] = bvec.s1;
blm[kg*NWG + ng*VWN + 2] = bvec.s2;
blm[kg*NWG + ng*VWN + 3] = bvec.s3;
blm[kg*NWG + ng*VWN + 4] = bvec.s4;
blm[kg*NWG + ng*VWN + 5] = bvec.s5;
blm[kg*NWG + ng*VWN + 6] = bvec.s6;
blm[kg*NWG + ng*VWN + 7] = bvec.s7;
blm[kg*NWG + ng*VWN + 8] = bvec.s8;
blm[kg*NWG + ng*VWN + 9] = bvec.s9;
blm[kg*NWG + ng*VWN + 10] = bvec.sA;
blm[kg*NWG + ng*VWN + 11] = bvec.sB;
blm[kg*NWG + ng*VWN + 12] = bvec.sC;
blm[kg*NWG + ng*VWN + 13] = bvec.sD;
blm[kg*NWG + ng*VWN + 14] = bvec.sE;
blm[kg*NWG + ng*VWN + 15] = bvec.sF;
const realMD bvec = bgm[idk*(b_ld/VWND) + idn + b_offset];
#if VWND == 1
blm[kg*NWGD + ng] = bvec;
#elif VWND == 2
blm[kg*NWGD + ng*VWND + 0] = bvec.x;
blm[kg*NWGD + ng*VWND + 1] = bvec.y;
#elif VWND == 4
blm[kg*NWGD + ng*VWND + 0] = bvec.x;
blm[kg*NWGD + ng*VWND + 1] = bvec.y;
blm[kg*NWGD + ng*VWND + 2] = bvec.z;
blm[kg*NWGD + ng*VWND + 3] = bvec.w;
#elif VWND == 8
blm[kg*NWGD + ng*VWND + 0] = bvec.s0;
blm[kg*NWGD + ng*VWND + 1] = bvec.s1;
blm[kg*NWGD + ng*VWND + 2] = bvec.s2;
blm[kg*NWGD + ng*VWND + 3] = bvec.s3;
blm[kg*NWGD + ng*VWND + 4] = bvec.s4;
blm[kg*NWGD + ng*VWND + 5] = bvec.s5;
blm[kg*NWGD + ng*VWND + 6] = bvec.s6;
blm[kg*NWGD + ng*VWND + 7] = bvec.s7;
#elif VWND == 16
blm[kg*NWGD + ng*VWND + 0] = bvec.s0;
blm[kg*NWGD + ng*VWND + 1] = bvec.s1;
blm[kg*NWGD + ng*VWND + 2] = bvec.s2;
blm[kg*NWGD + ng*VWND + 3] = bvec.s3;
blm[kg*NWGD + ng*VWND + 4] = bvec.s4;
blm[kg*NWGD + ng*VWND + 5] = bvec.s5;
blm[kg*NWGD + ng*VWND + 6] = bvec.s6;
blm[kg*NWGD + ng*VWND + 7] = bvec.s7;
blm[kg*NWGD + ng*VWND + 8] = bvec.s8;
blm[kg*NWGD + ng*VWND + 9] = bvec.s9;
blm[kg*NWGD + ng*VWND + 10] = bvec.sA;
blm[kg*NWGD + ng*VWND + 11] = bvec.sB;
blm[kg*NWGD + ng*VWND + 12] = bvec.sC;
blm[kg*NWGD + ng*VWND + 13] = bvec.sD;
blm[kg*NWGD + ng*VWND + 14] = bvec.sE;
blm[kg*NWGD + ng*VWND + 15] = bvec.sF;
#endif
if (b_conjugate) {
for (int vn=0; vn<VWN; ++vn) {
COMPLEX_CONJUGATE(blm[kg*NWG + ng*VWN + vn]);
for (int vn=0; vn<VWND; ++vn) {
COMPLEX_CONJUGATE(blm[kg*NWGD + ng*VWND + vn]);
}
}
}
@ -153,23 +225,23 @@ inline void GlobalToLocalDirectB(const __global realN* restrict bgm, __local rea
// Caches on-chip local memory into per-thread private memory (registers). This function is specific
// for caching the A input matrix.
inline void LocalToPrivateDirectA(__local real* alm, real apm[MWI], const int kg,
inline void LocalToPrivateDirectA(__local real* alm, real apm[MWID], const int kg,
const int a_transpose) {
#pragma unroll
for (int mi=0; mi<MWI; ++mi) {
const int mg = mi + get_local_id(0)*MWI;
const int index = (a_transpose) ? mg*KWG + kg : kg*MWG + mg;
for (int mi=0; mi<MWID; ++mi) {
const int mg = mi + get_local_id(0)*MWID;
const int index = (a_transpose) ? mg*KWGD + kg : kg*MWGD + mg;
apm[mi] = alm[index];
}
}
// Same as above, but now for the B input matrix
inline void LocalToPrivateDirectB(__local real* blm, real bpm[NWI], const int kg,
inline void LocalToPrivateDirectB(__local real* blm, real bpm[NWID], const int kg,
const int b_transpose) {
#pragma unroll
for (int ni=0; ni<NWI; ++ni) {
const int ng = ni + get_local_id(1)*NWI;
const int index = (b_transpose) ? ng*KWG + kg : kg*NWG + ng;
for (int ni=0; ni<NWID; ++ni) {
const int ng = ni + get_local_id(1)*NWID;
const int index = (b_transpose) ? ng*KWGD + kg : kg*NWGD + ng;
bpm[ni] = blm[index];
}
}
@ -177,11 +249,11 @@ inline void LocalToPrivateDirectB(__local real* blm, real bpm[NWI], const int kg
// =================================================================================================
// Initializes the accumulation registers to zero
inline void InitAccRegistersDirect(real cpm[NWI][MWI]) {
inline void InitAccRegistersDirect(real cpm[NWID][MWID]) {
#pragma unroll
for (int mi=0; mi<MWI; ++mi) {
for (int mi=0; mi<MWID; ++mi) {
#pragma unroll
for (int ni=0; ni<NWI; ++ni) {
for (int ni=0; ni<NWID; ++ni) {
SetToZero(cpm[ni][mi]);
}
}
@ -190,11 +262,11 @@ inline void InitAccRegistersDirect(real cpm[NWI][MWI]) {
// =================================================================================================
// Performs the actual computation: Cpm += Apm * Bpm
inline void MultiplyAccumulateDirect(real cpm[NWI][MWI], real apm[MWI], real bpm[NWI]) {
inline void MultiplyAccumulateDirect(real cpm[NWID][MWID], real apm[MWID], real bpm[NWID]) {
#pragma unroll
for (int ni=0; ni<NWI; ++ni) {
for (int ni=0; ni<NWID; ++ni) {
#pragma unroll
for (int mi=0; mi<MWI; ++mi) {
for (int mi=0; mi<MWID; ++mi) {
MultiplyAdd(cpm[ni][mi], apm[mi], bpm[ni]);
}
}
@ -204,18 +276,18 @@ inline void MultiplyAccumulateDirect(real cpm[NWI][MWI], real apm[MWI], real bpm
// Merges the results in Cpm with the global array in Cgm. This also performs the multiplication
// with the constants: Cgm = alpha*A*B + beta*Cgm = alpha*Cpm + beta*Cgm
inline void StoreResultsDirect(__global real* cgm, real cpm[NWI][MWI],
inline void StoreResultsDirect(__global real* cgm, real cpm[NWID][MWID],
const int kSizeM, const int kSizeN,
const real alpha, const real beta,
const int c_ld, const int c_offset, const int c_transpose) {
#pragma unroll
for (int ni=0; ni<NWI; ++ni) {
for (int ni=0; ni<NWID; ++ni) {
#pragma unroll
for (int mi=0; mi<MWI; ++mi) {
int mg = mi + get_local_id(0)*MWI;
int ng = ni + get_local_id(1)*NWI;
int idm = mg + GetGroupID0() * MWG;
int idn = ng + GetGroupID1() * NWG;
for (int mi=0; mi<MWID; ++mi) {
int mg = mi + get_local_id(0)*MWID;
int ng = ni + get_local_id(1)*NWID;
int idm = mg + GetGroupID0() * MWGD;
int idn = ng + GetGroupID1() * NWGD;
// Determines the destination index
const int c_index = (c_transpose) ? idm*c_ld + idn : idn*c_ld + idm;
@ -231,12 +303,12 @@ inline void StoreResultsDirect(__global real* cgm, real cpm[NWI][MWI],
// =================================================================================================
// Main entry point of the kernel. This is the direct version without restrictions.
__attribute__((reqd_work_group_size(MDIMC, NDIMC, 1)))
__attribute__((reqd_work_group_size(MDIMCD, NDIMCD, 1)))
__kernel void XgemmDirect(const int kSizeM, const int kSizeN, const int kSizeK,
const real_arg arg_alpha,
const real_arg arg_beta,
const __global realM* restrict agm, const int a_offset, const int a_ld,
const __global realN* restrict bgm, const int b_offset, const int b_ld,
const __global realMD* restrict agm, const int a_offset, const int a_ld,
const __global realND* restrict bgm, const int b_offset, const int b_ld,
__global real* cgm, const int c_offset, const int c_ld,
const int a_transpose, const int b_transpose, const int c_transpose,
const int a_conjugate, const int b_conjugate) {
@ -248,40 +320,40 @@ __kernel void XgemmDirect(const int kSizeM, const int kSizeN, const int kSizeK,
const __global real* restrict bgms = (const __global real* restrict) bgm;
// Allocates workgroup-private memory (local memory)
__local real alm[KWG * MWG];
__local real blm[KWG * NWG];
__local real alm[KWGD * MWGD];
__local real blm[KWGD * NWGD];
// Combined thread identifier (volatile to disable caching)
volatile int tid = get_local_id(0) + MDIMC*get_local_id(1);
volatile int tid = get_local_id(0) + MDIMCD*get_local_id(1);
// Allocates workitem-private memory (registers)
real apm[MWI];
real bpm[NWI];
real cpm[NWI][MWI];
real apm[MWID];
real bpm[NWID];
real cpm[NWID][MWID];
// Initializes the accumulation registers
InitAccRegistersDirect(cpm);
// The faster version of GEMM is not allowed on the (incomplete) borders. Therefore, this section
// processes only the main parts: output blocks of MWG by NWG.
const int idm = get_local_id(0) * MWI + GetGroupID0() * MWG;
const int idn = get_local_id(1) * NWI + GetGroupID1() * NWG;
if ((idm < (kSizeM/MWG)*MWG) && (idn < (kSizeN/NWG)*NWG) &&
(a_ld % VWM == 0) && (b_ld % VWN == 0)) {
// processes only the main parts: output blocks of MWGD by NWGD.
const int idm = get_local_id(0) * MWID + GetGroupID0() * MWGD;
const int idn = get_local_id(1) * NWID + GetGroupID1() * NWGD;
if ((idm < (kSizeM/MWGD)*MWGD) && (idn < (kSizeN/NWGD)*NWGD) &&
(a_ld % VWMD == 0) && (b_ld % VWND == 0)) {
// Loops over all complete workgroup tiles
int kwg = 0;
for (; kwg < (kSizeK/KWG) * KWG; kwg+=KWG) {
for (; kwg < (kSizeK/KWGD) * KWGD; kwg+=KWGD) {
// Loads data: off-chip --> local (matrix A and B)
GlobalToLocalDirectA(agm, alm, a_ld, a_offset, tid, kwg, a_transpose, a_conjugate);
GlobalToLocalDirectB(bgm, blm, b_ld, b_offset, tid, kwg, b_transpose, b_conjugate);
barrier(CLK_LOCAL_MEM_FENCE);
// Loops over all workitem tiles, unrolled by a factor KWI
for (int pwi=0; pwi<KWG; pwi+=KWI) {
// Loops over all workitem tiles, unrolled by a factor KWID
for (int pwi=0; pwi<KWGD; pwi+=KWID) {
#pragma unroll
for (int pit=0; pit<KWI; ++pit) {
for (int pit=0; pit<KWID; ++pit) {
int kg = pwi + pit;
// Loads data: local --> private (matrix A)
@ -303,7 +375,7 @@ __kernel void XgemmDirect(const int kSizeM, const int kSizeN, const int kSizeK,
// Loads A into register memory
#pragma unroll
for (int mi=0; mi<MWI; ++mi) {
for (int mi=0; mi<MWID; ++mi) {
const int a_index = (a_transpose) ? (idm + mi)*a_ld + idk : idk*a_ld + (idm + mi);
apm[mi] = agms[a_index + a_offset];
if (a_conjugate) { COMPLEX_CONJUGATE(apm[mi]); }
@ -311,7 +383,7 @@ __kernel void XgemmDirect(const int kSizeM, const int kSizeN, const int kSizeK,
// Loads B into register memory
#pragma unroll
for (int ni=0; ni<NWI; ++ni) {
for (int ni=0; ni<NWID; ++ni) {
const int b_index = (b_transpose) ? (idn + ni)*b_ld + idk : idk*b_ld + (idn + ni);
bpm[ni] = bgms[b_index + b_offset];
if (b_conjugate) { COMPLEX_CONJUGATE(bpm[ni]); }
@ -321,10 +393,6 @@ __kernel void XgemmDirect(const int kSizeM, const int kSizeN, const int kSizeK,
MultiplyAccumulateDirect(cpm, apm, bpm);
}
#if GLOBAL_MEM_FENCE == 1
barrier(CLK_GLOBAL_MEM_FENCE);
#endif
// Stores a tile of results and performs the multiplication with alpha and beta
StoreResultsDirect(cgm, cpm, kSizeM, kSizeN, alpha, beta, c_ld, c_offset, c_transpose);
}
@ -337,7 +405,7 @@ __kernel void XgemmDirect(const int kSizeM, const int kSizeN, const int kSizeK,
// Loads A into register memory
#pragma unroll
for (int mi=0; mi<MWI; ++mi) {
for (int mi=0; mi<MWID; ++mi) {
if (idm + mi < kSizeM) {
const int a_index = (a_transpose) ? (idm + mi)*a_ld + idk : idk*a_ld + (idm + mi);
apm[mi] = agms[a_index + a_offset];
@ -350,7 +418,7 @@ __kernel void XgemmDirect(const int kSizeM, const int kSizeN, const int kSizeK,
// Loads B into register memory
#pragma unroll
for (int ni=0; ni<NWI; ++ni) {
for (int ni=0; ni<NWID; ++ni) {
if (idn + ni < kSizeN) {
const int b_index = (b_transpose) ? (idn + ni)*b_ld + idk : idk*b_ld + (idn + ni);
bpm[ni] = bgms[b_index + b_offset];
@ -367,9 +435,9 @@ __kernel void XgemmDirect(const int kSizeM, const int kSizeN, const int kSizeK,
// Stores the results
#pragma unroll
for (int ni=0; ni<NWI; ++ni) {
for (int ni=0; ni<NWID; ++ni) {
#pragma unroll
for (int mi=0; mi<MWI; ++mi) {
for (int mi=0; mi<MWID; ++mi) {
if ((idm + mi) < kSizeM && (idn + ni) < kSizeN) {
// Determines the destination index

13
src/routines/level3/xgemm.cpp
View File

@ -22,7 +22,8 @@ namespace clblast {
// Constructor: forwards to base class constructor
template <typename T>
Xgemm<T>::Xgemm(Queue &queue, EventPointer event, const std::string &name):
Routine(queue, event, name, {"Copy","Pad","Transpose","Padtranspose","Xgemm"}, PrecisionValue<T>()) {
Routine(queue, event, name, {"Copy","Pad","Transpose","Padtranspose","Xgemm", "XgemmDirect"},
PrecisionValue<T>()) {
source_string_ =
#include "../../kernels/level3/level3.opencl"
#include "../../kernels/level3/copy_fast.opencl"
@ -299,13 +300,13 @@ StatusCode Xgemm<T>::GemmDirect(const size_t m, const size_t n, const size_t k,
kernel.SetArgument(18, static_cast<int>(b_conjugate));
// Computes the global and local thread sizes
const auto m_ceiled = Ceil(m, db_["MWG"]);
const auto n_ceiled = Ceil(n, db_["NWG"]);
const auto m_ceiled = Ceil(m, db_["MWGD"]);
const auto n_ceiled = Ceil(n, db_["NWGD"]);
const auto global = std::vector<size_t>{
(m_ceiled * db_["MDIMC"]) / db_["MWG"],
(n_ceiled * db_["NDIMC"]) / db_["NWG"]
(m_ceiled * db_["MDIMCD"]) / db_["MWGD"],
(n_ceiled * db_["NDIMCD"]) / db_["NWGD"]
};
const auto local = std::vector<size_t>{db_["MDIMC"], db_["NDIMC"]};
const auto local = std::vector<size_t>{db_["MDIMCD"], db_["NDIMCD"]};
// Launches the kernel
auto status = RunKernel(kernel, queue_, device_, global, local, event_);