diff options
-rw-r--r-- | src/database/database.cpp | 2 | ||||
-rw-r--r-- | src/database/database.hpp | 1 | ||||
-rw-r--r-- | src/database/kernels/xgemm_direct.hpp | 76 | ||||
-rw-r--r-- | src/kernels/level3/xgemm_direct.opencl | 358 | ||||
-rw-r--r-- | src/routines/level3/xgemm.cpp | 13 |
5 files changed, 299 insertions, 151 deletions
diff --git a/src/database/database.cpp b/src/database/database.cpp index 34c44a29..2696fb9b 100644 --- a/src/database/database.cpp +++ b/src/database/database.cpp @@ -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, diff --git a/src/database/database.hpp b/src/database/database.hpp index a6ab49c5..7c0afb46 100644 --- a/src/database/database.hpp +++ b/src/database/database.hpp @@ -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; diff --git a/src/database/kernels/xgemm_direct.hpp b/src/database/kernels/xgemm_direct.hpp new file mode 100644 index 00000000..76055ef2 --- /dev/null +++ b/src/database/kernels/xgemm_direct.hpp @@ -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 diff --git a/src/kernels/level3/xgemm_direct.opencl b/src/kernels/level3/xgemm_direct.opencl index fb5972ba..801887dd 100644 --- a/src/kernels/level3/xgemm_direct.opencl +++ b/src/kernels/level3/xgemm_direct.opencl @@ -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 diff --git a/src/routines/level3/xgemm.cpp b/src/routines/level3/xgemm.cpp index 90e43fe4..2fb9f1fd 100644 --- a/src/routines/level3/xgemm.cpp +++ b/src/routines/level3/xgemm.cpp @@ -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_); |