GemmSpecialCases.cpp

/* ************************************************************************
* Copyright 2015 Advanced Micro Devices, Inc.
*
* Licensed under the Apache License, Version 2.0 (the "License");
* you may not use this file except in compliance with the License.
* You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
* ************************************************************************/

#include "GemmSpecialCases.h"
#include "UserGemmKernelSources/UserGemmKernelSourceIncludes.h"
#include "UserGemmKernelSources/UserGemmClKernels.h"
#include "xgemm.h" //helper functions defined in xgemm.cpp
#include "AutoGemmIncludes/AutoGemmClKernels.h"
#include "AutoGemmIncludes/AutoGemmKernelSources.h"
#include "AutoGemmIncludes/AutoGemmKernelBinaries.h"

/******************************************************************************
* Check OpenCL Errors
*****************************************************************************/
#define CL_CHECK(RET) \
  if(RET != CL_SUCCESS) { \
    printf("OpenCL error %i on line %u\n", RET, __LINE__); \
    assert(false); \
    }
/*
template<typename precision>
clblasStatus SGEMM_SPLIT_CALLS(
	cl_kernel  ClKernel, clblasOrder order,
	unsigned int tile_size, unsigned int WG_size,
	unsigned int M_split_factor,
	unsigned int N_split_factor,
	unsigned int K_split_factor,
	clblasTranspose transA,
	clblasTranspose transB,
	cl_uint M, cl_uint N, cl_uint K,
	precision alpha,
	cl_mem A, cl_uint offA, cl_uint lda,
	cl_mem B, cl_uint offB, cl_uint ldb,
	precision beta,
	cl_mem C, cl_uint offC, cl_uint ldc,
	cl_uint numCommandQueues,
	cl_command_queue *commandQueues,
	cl_uint numEventsInWaitList,
	const cl_event *eventWaitList,
	cl_event *events);
*/
template<typename precision>
clblasStatus GEMM_SPLIT_CALLS(
	cl_kernel  ClKernel, clblasOrder order,
	unsigned int tile_size, unsigned int WG_size,
	unsigned int M_split_factor,
	unsigned int N_split_factor,
	unsigned int K_split_factor,
	clblasTranspose transA,
	clblasTranspose transB,
	cl_uint M, cl_uint N, cl_uint K,
	precision alpha,
	cl_mem A, cl_uint offA, cl_uint lda,
	cl_mem B, cl_uint offB, cl_uint ldb,
	precision beta,
	cl_mem C, cl_uint offC, cl_uint ldc,
	cl_uint numCommandQueues,
	cl_command_queue *commandQueues,
	cl_uint numEventsInWaitList,
	const cl_event *eventWaitList,
	cl_event *events)
{
	//for example, when M=N=K=8192 in GEMM col NT
	//we are gonna call 16 GEMMs 
	//each GEMM has M=N=K=4096
	//note are direct GEMM call has a 0.7 TFLOPS performance

	//     [ A11 | A12 | A13 | A14 ]      [ B11 | B12 | B13 | B14 ]      [ C11 | C12 ]
	// A = [ A21 | A22 | A23 | A24 ]  B = [ B21 | B22 | B23 | B24 ]  C = [ C21 | C22 ] 

	// 16 GEMMs are
	// #01: C11 = a*A11*B11 + b*C11
	// #02: C11 = a*A12*B12 + 1*C11
	// #03: C11 = a*A13*B13 + 1*C11
	// #04: C11 = a*A14*B14 + 1*C11 now we are done with C11

	// #05: C12 = a*A11*B21 + b*C12
	// #06: C12 = a*A12*B22 + 1*C12
	// #07: C12 = a*A12*B22 + 1*C12
	// #08: C12 = a*A12*B22 + 1*C12 now we are done with C12

	// #09: C21 = a*A21*B11 + b*C21
	// #10: C21 = a*A22*B12 + 1*C21
	// #11: C21 = a*A23*B13 + 1*C21
	// #12: C21 = a*A24*B14 + 1*C21 now we are done with C21

	// #13: C22 = a*A21*B21 + b*C22
	// #14: C22 = a*A22*B22 + 1*C22
	// #15: C22 = a*A23*B23 + 1*C22
	// #16: C22 = a*A24*B24 + 1*C22 now we are done with C22

	if (transA == clblasNoTrans && transB == clblasTrans)
	{
		unsigned int small_M = M / M_split_factor;
		unsigned int small_N = N / N_split_factor;
		unsigned int small_K = K / K_split_factor;

		size_t GlobalX = ((small_M - 1) / tile_size + 1) * WG_size;
		size_t GlobalY = ((small_N - 1) / tile_size + 1) * WG_size;
		size_t gs[2] = { GlobalX, GlobalY };
		size_t wgsize[2] = { WG_size, WG_size };
		cl_int error = 0;

		precision betaone = 1;

		error = clSetKernelArg(ClKernel, 5, sizeof(cl_uint), &small_M);
		assert(error == CL_SUCCESS);
		error = clSetKernelArg(ClKernel, 6, sizeof(cl_uint), &small_N);
		assert(error == CL_SUCCESS);
		error = clSetKernelArg(ClKernel, 7, sizeof(cl_uint), &small_K);
		assert(error == CL_SUCCESS);

		for (int M_split_index = 0; M_split_index < M_split_factor; M_split_index++)
		{
			for (int N_split_index = 0; N_split_index < N_split_factor; N_split_index++)
			{
				unsigned int offc_C = ldc*N / N_split_factor * N_split_index + M / M_split_factor * M_split_index + offC;
				error = clSetKernelArg(ClKernel, 13, sizeof(cl_uint), &offc_C);
				assert(error == CL_SUCCESS);

				for (int K_split_index = 0; K_split_index < K_split_factor; K_split_index++)
				{
					unsigned int offa_A = (M / M_split_factor * M_split_index) + (lda * K / K_split_factor * K_split_index) + offA;
					unsigned int offb_B = (N / N_split_factor * N_split_index) + (ldb * K / K_split_factor * K_split_index) + offB;
					error = clSetKernelArg(ClKernel, 11, sizeof(cl_uint), &offa_A);
					assert(error == CL_SUCCESS);
					error = clSetKernelArg(ClKernel, 12, sizeof(cl_uint), &offb_B);
					assert(error == CL_SUCCESS);

					if (K_split_index == 0)
					{
						error = clSetKernelArg(ClKernel, 4, sizeof(precision), &(beta));
						assert(error == CL_SUCCESS);

						if (M_split_index == 0 && N_split_index == 0)
						{
							//very first GEMM call
							if ((M_split_factor == 1) && (N_split_factor == 1) && (K_split_factor == 1))
							{
								//also very last GEMM call
								error = clEnqueueNDRangeKernel(commandQueues[0], ClKernel, 2, NULL,
									gs, wgsize, numEventsInWaitList, eventWaitList, &events[0]);
								assert(error == CL_SUCCESS);
							}
							else
							{
								error = clEnqueueNDRangeKernel(commandQueues[0], ClKernel, 2, NULL,
									gs, wgsize, numEventsInWaitList, eventWaitList, NULL);
								assert(error == CL_SUCCESS);
							}
						}
						else
						{
							error = clEnqueueNDRangeKernel(commandQueues[0], ClKernel, 2, NULL,
								gs, wgsize, 0, NULL, NULL);
							assert(error == CL_SUCCESS);
						}
					}
					else
					{
						error = clSetKernelArg(ClKernel, 4, sizeof(precision), &betaone);
						assert(error == CL_SUCCESS);

						if ((M_split_index == (M_split_factor - 1)) && (N_split_index == (N_split_factor - 1)) && (K_split_index == (K_split_factor - 1)))
						{
							//very last GEMM call
							error = clEnqueueNDRangeKernel(commandQueues[0], ClKernel, 2, NULL,
								gs, wgsize, 0, NULL, events);
							assert(error == CL_SUCCESS);
						}
						else
						{
							error = clEnqueueNDRangeKernel(commandQueues[0], ClKernel, 2, NULL,
								gs, wgsize, 0, NULL, NULL);
							assert(error == CL_SUCCESS);
						}
					}
				}
			}
		}

		return clblasSuccess;
	}
	return clblasNotImplemented;
}

clblasStatus SGEMM_mod1024(
	clblasTranspose transA,
	clblasTranspose transB,
	cl_uint M, cl_uint N, cl_uint K,
	float alpha,
	cl_mem A, cl_uint offA, cl_uint lda,
	cl_mem B, cl_uint offB, cl_uint ldb,
	float beta,
	cl_mem C, cl_uint offC, cl_uint ldc,
	cl_uint numCommandQueues,
	cl_command_queue *commandQueues,
	cl_uint numEventsInWaitList,
	const cl_event *eventWaitList,
	cl_event *events,
	bool &specialCaseHandled)
{
	const char *tileKernelSource = NULL;
	cl_kernel   tileClKernel = NULL;
	size_t tileKernelBinarySize = 0;
	cl_int err;


	const unsigned char *tileKernelBinary = NULL;

	clblasStatus status;


	//split the kernel calls to handle sgemm NT perf drop at big multiples of 1024
	if ((lda % 1024 == 0) && (ldb % 1024 == 0) && (K > lda / 4))
	{
		if ((lda == ldb) && (lda >= 4096) && (lda <= 8192)) // between 4096 and 8192 for now
		{
			if (lda != 6144)// 6144 is handled by 96 x 96 kernel
			{
				// we are going to call 16 GEMMs with M=M/2, N=N/2, K=K/4
				// each GEMM requires M%128 == 0, N%128 == 0, K%16 == 0
				if (M % 256 == 0 && N % 256 == 0 && K % 64 == 0)
				{
					if (!((transA == clblasNoTrans) && (transB == clblasTrans)))
						return clblasNotImplemented;

					specialCaseHandled = true;
					unsigned int M_split_factor;
					unsigned int N_split_factor;
					unsigned int K_split_factor;

					if (lda < 7168)
					{
						M_split_factor = 1;
						N_split_factor = 1;
						K_split_factor = 1;
					}
					else
					{
						//7168, 8192
						M_split_factor = 2;
						N_split_factor = 2;
						K_split_factor = 4;
					}

					tileKernelSource = sgemm_Col_NT_B1_MX128_NX128_KX16_src;
					tileClKernel = sgemm_Col_NT_B1_MX128_NX128_KX16_clKernel;
					tileKernelBinary = sgemm_Col_NT_B1_MX128_NX128_KX16_bin;
					tileKernelBinarySize = sgemm_Col_NT_B1_MX128_NX128_KX16_binSize;

					makeGemmKernel(&tileClKernel, commandQueues[0], tileKernelSource, User_srcBuildOptions, &tileKernelBinary, &tileKernelBinarySize, User_binBuildOptions);

					err = clSetKernelArg(tileClKernel, 0, sizeof(cl_mem), &A);
					CL_CHECK(err);
					err = clSetKernelArg(tileClKernel, 1, sizeof(cl_mem), &B);
					CL_CHECK(err);
					err = clSetKernelArg(tileClKernel, 2, sizeof(cl_mem), &C);
					CL_CHECK(err);
					err = clSetKernelArg(tileClKernel, 3, sizeof(cl_float), &alpha);
					CL_CHECK(err);
					err = clSetKernelArg(tileClKernel, 4, sizeof(cl_float), &beta);
					CL_CHECK(err);
					err = clSetKernelArg(tileClKernel, 5, sizeof(cl_uint), &M);
					CL_CHECK(err);
					err = clSetKernelArg(tileClKernel, 6, sizeof(cl_uint), &N);
					CL_CHECK(err);
					err = clSetKernelArg(tileClKernel, 7, sizeof(cl_uint), &K);
					CL_CHECK(err);
					err = clSetKernelArg(tileClKernel, 8, sizeof(cl_uint), &lda);
					CL_CHECK(err);
					err = clSetKernelArg(tileClKernel, 9, sizeof(cl_uint), &ldb);
					CL_CHECK(err);
					err = clSetKernelArg(tileClKernel, 10, sizeof(cl_uint), &ldc);
					CL_CHECK(err);
					err = clSetKernelArg(tileClKernel, 11, sizeof(cl_uint), &offA);
					CL_CHECK(err);
					err = clSetKernelArg(tileClKernel, 12, sizeof(cl_uint), &offB);
					CL_CHECK(err);
					err = clSetKernelArg(tileClKernel, 13, sizeof(cl_uint), &offC);
					CL_CHECK(err);

					status = GEMM_SPLIT_CALLS(
						tileClKernel, clblasColumnMajor,
						128, 16,
						M_split_factor,
						N_split_factor, K_split_factor,
						transA,
						transB,
						M, N, K,
						alpha,
						A, offA, lda,
						B, offB, ldb,
						beta,
						C, offC, ldc,
						numCommandQueues,
						commandQueues,
						numEventsInWaitList,
						eventWaitList,
						events);


					return status;
				}
			}
			else
			{
				// lda == ldb == 6144
				// we are going to call 4 GEMMs each with K = K/4
				if (M % 96 == 0 && N % 96 == 0 && K % 64 == 0)
				{
					if (!((transA == clblasNoTrans) && (transB == clblasTrans)))
						return clblasNotImplemented;

					specialCaseHandled = true;
					unsigned int M_split_factor = 1;
					unsigned int N_split_factor = 1;
					unsigned int K_split_factor = 4;



					tileKernelSource = sgemm_Col_NT_B1_MX096_NX096_KX16_src;
					tileClKernel = sgemm_Col_NT_B1_MX096_NX096_KX16_clKernel;
					tileKernelBinary = sgemm_Col_NT_B1_MX096_NX096_KX16_bin;
					tileKernelBinarySize = sgemm_Col_NT_B1_MX096_NX096_KX16_binSize;

					makeGemmKernel(&tileClKernel, commandQueues[0], tileKernelSource, User_srcBuildOptions, &tileKernelBinary, &tileKernelBinarySize, User_binBuildOptions);

					err = clSetKernelArg(tileClKernel, 0, sizeof(cl_mem), &A);
					CL_CHECK(err);
					err = clSetKernelArg(tileClKernel, 1, sizeof(cl_mem), &B);
					CL_CHECK(err);
					err = clSetKernelArg(tileClKernel, 2, sizeof(cl_mem), &C);
					CL_CHECK(err);
					err = clSetKernelArg(tileClKernel, 3, sizeof(cl_float), &alpha);
					CL_CHECK(err);
					err = clSetKernelArg(tileClKernel, 4, sizeof(cl_float), &beta);
					CL_CHECK(err);
					err = clSetKernelArg(tileClKernel, 5, sizeof(cl_uint), &M);
					CL_CHECK(err);
					err = clSetKernelArg(tileClKernel, 6, sizeof(cl_uint), &N);
					CL_CHECK(err);
					err = clSetKernelArg(tileClKernel, 7, sizeof(cl_uint), &K);
					CL_CHECK(err);
					err = clSetKernelArg(tileClKernel, 8, sizeof(cl_uint), &lda);
					CL_CHECK(err);
					err = clSetKernelArg(tileClKernel, 9, sizeof(cl_uint), &ldb);
					CL_CHECK(err);
					err = clSetKernelArg(tileClKernel, 10, sizeof(cl_uint), &ldc);
					CL_CHECK(err);
					err = clSetKernelArg(tileClKernel, 11, sizeof(cl_uint), &offA);
					CL_CHECK(err);
					err = clSetKernelArg(tileClKernel, 12, sizeof(cl_uint), &offB);
					CL_CHECK(err);
					err = clSetKernelArg(tileClKernel, 13, sizeof(cl_uint), &offC);
					CL_CHECK(err);


					status = GEMM_SPLIT_CALLS(
						tileClKernel, clblasColumnMajor,
						96, 16,
						M_split_factor,
						N_split_factor, K_split_factor,
						transA,
						transB,
						M, N, K,
						alpha,
						A, offA, lda,
						B, offB, ldb,
						beta,
						C, offC, ldc,
						numCommandQueues,
						commandQueues,
						numEventsInWaitList,
						eventWaitList,
						events);


					return status;
				}
			}
		}
	}

	return clblasNotImplemented;
}


clblasStatus SGEMM_SPLIT64_32(
	clblasTranspose transA,
	clblasTranspose transB,
	cl_uint M, cl_uint N, cl_uint K,
	float alpha,
	cl_mem A, cl_uint offA, cl_uint lda,
	cl_mem B, cl_uint offB, cl_uint ldb,
	float beta,
	cl_mem C, cl_uint offC, cl_uint ldc,
	cl_uint numCommandQueues,
	cl_command_queue *commandQueues,
	cl_uint numEventsInWaitList,
	const cl_event *eventWaitList,
	cl_event *events,
	bool &specialCaseHandled)
{
	//all the mod32 sizes that is not mod64 or mod96 ranging from 1184 to 3872
	//non mod32 cases are not implemented in this approach and are of less interest
	const char *tileKernelSource = NULL;
	const char *rowKernelSource = NULL;
	const char *columnKernelSource = NULL;
	const char *singleKernelSource = NULL;

	cl_kernel   tileClKernel = NULL;
	cl_kernel   rowClKernel = NULL;
	cl_kernel   columnClKernel = NULL;
	cl_kernel   singleClKernel = NULL;

	const unsigned char *tileKernelBinary = NULL;
	const unsigned char *rowKernelBinary = NULL;
	const unsigned char *columnKernelBinary = NULL;
	const unsigned char *singleKernelBinary = NULL;

	size_t tileKernelBinarySize = 0;
	size_t rowKernelBinarySize = 0;
	size_t columnKernelBinarySize = 0;
	size_t singleKernelBinarySize = 0;

	cl_int err;

	if ((M >= 1184 && N >= 1184) && (M <= 3872 && N <= 3872) && (M % 64 != 0 && N % 64 != 0) && (M % 96 != 0 && N % 96 != 0) && (K % 16 == 0))
	{
		if ((M % 32 == 0 && N % 32 == 0) && (transA == clblasNoTrans && transB == clblasTrans))
		{
			specialCaseHandled = true;
			//execute the kernels

			//GlobalX = ((Mvalue - 1) / 64) * 16
			//GlobalY = ((Nvalue - 1) / 64) * 16
			size_t GlobalX = ((M - 1) / 64) * 16;
			size_t GlobalY = ((N - 1) / 64) * 16;
			size_t gs[2] = { GlobalX, GlobalY };
			size_t wgsize[2] = { 16, 16 };

			tileKernelSource = sgemm_Col_NT_B1_MX064_NX064_KX16_src;
			tileClKernel = sgemm_Col_NT_B1_MX064_NX064_KX16_clKernel;
			tileKernelBinary = sgemm_Col_NT_B1_MX064_NX064_KX16_bin;
			tileKernelBinarySize = sgemm_Col_NT_B1_MX064_NX064_KX16_binSize;

			rowKernelSource = sgemm_Col_NT_B1_MX032_NX064_KX16_ROW_src;
			rowClKernel = sgemm_Col_NT_B1_MX032_NX064_KX16_ROW_clKernel;
			rowKernelBinary = sgemm_Col_NT_B1_MX032_NX064_KX16_ROW_bin;
			rowKernelBinarySize = sgemm_Col_NT_B1_MX032_NX064_KX16_ROW_binSize;

			columnKernelSource = sgemm_Col_NT_B1_MX064_NX032_KX16_COLUMN_src;
			columnClKernel = sgemm_Col_NT_B1_MX064_NX032_KX16_COLUMN_clKernel;
			columnKernelBinary = sgemm_Col_NT_B1_MX064_NX032_KX16_COLUMN_bin;
			columnKernelBinarySize = sgemm_Col_NT_B1_MX064_NX032_KX16_COLUMN_binSize;

			singleKernelSource = sgemm_Col_NT_B1_MX032_NX032_KX16_SINGLE_src;
			singleClKernel = sgemm_Col_NT_B1_MX032_NX032_KX16_SINGLE_clKernel;
			singleKernelBinary = sgemm_Col_NT_B1_MX032_NX032_KX16_SINGLE_bin;
			singleKernelBinarySize = sgemm_Col_NT_B1_MX032_NX032_KX16_SINGLE_binSize;

			cl_kernel   Kernels[4] = { tileClKernel, rowClKernel, columnClKernel, singleClKernel };


			makeGemmKernel(&tileClKernel,   commandQueues[0], tileKernelSource,   User_srcBuildOptions, &tileKernelBinary,   &tileKernelBinarySize,   User_binBuildOptions);
			makeGemmKernel(&rowClKernel,    commandQueues[0], rowKernelSource,    User_srcBuildOptions, &rowKernelBinary,    &rowKernelBinarySize,    User_binBuildOptions);
			makeGemmKernel(&columnClKernel, commandQueues[0], columnKernelSource, User_srcBuildOptions, &columnKernelBinary, &columnKernelBinarySize, User_binBuildOptions);
			makeGemmKernel(&singleClKernel, commandQueues[0], singleKernelSource, User_srcBuildOptions, &singleKernelBinary, &singleKernelBinarySize, User_binBuildOptions);

			for (int i = 0; i < 4; i++)
			{
				err = clSetKernelArg(Kernels[i], 0, sizeof(cl_mem), &A);
				CL_CHECK(err);
				err = clSetKernelArg(Kernels[i], 1, sizeof(cl_mem), &B);
				CL_CHECK(err);
				err = clSetKernelArg(Kernels[i], 2, sizeof(cl_mem), &C);
				CL_CHECK(err);
				err = clSetKernelArg(Kernels[i], 3, sizeof(cl_float), &alpha);
				CL_CHECK(err);
				err = clSetKernelArg(Kernels[i], 4, sizeof(cl_float), &beta);
				CL_CHECK(err);
				err = clSetKernelArg(Kernels[i], 5, sizeof(cl_uint), &M);
				CL_CHECK(err);
				err = clSetKernelArg(Kernels[i], 6, sizeof(cl_uint), &N);
				CL_CHECK(err);
				err = clSetKernelArg(Kernels[i], 7, sizeof(cl_uint), &K);
				CL_CHECK(err);
				err = clSetKernelArg(Kernels[i], 8, sizeof(cl_uint), &lda);
				CL_CHECK(err);
				err = clSetKernelArg(Kernels[i], 9, sizeof(cl_uint), &ldb);
				CL_CHECK(err);
				err = clSetKernelArg(Kernels[i], 10, sizeof(cl_uint), &ldc);
				CL_CHECK(err);
				err = clSetKernelArg(Kernels[i], 11, sizeof(cl_uint), &offA);
				CL_CHECK(err);
				err = clSetKernelArg(Kernels[i], 12, sizeof(cl_uint), &offB);
				CL_CHECK(err);
				err = clSetKernelArg(Kernels[i], 13, sizeof(cl_uint), &offC);
				CL_CHECK(err);
			}

			err = clEnqueueNDRangeKernel(commandQueues[0], Kernels[0], 2, NULL, gs, wgsize, numEventsInWaitList, eventWaitList, NULL);

			gs[0] = 16;
			err |= clEnqueueNDRangeKernel(commandQueues[0], Kernels[1], 2, NULL, gs, wgsize, 0, NULL, NULL);

			gs[1] = 16;
			gs[0] = GlobalX;
			err |= clEnqueueNDRangeKernel(commandQueues[0], Kernels[2], 2, NULL, gs, wgsize, 0, NULL, NULL);

			gs[0] = 16; gs[1] = 16;
			err |= clEnqueueNDRangeKernel(commandQueues[0], Kernels[3], 2, NULL, gs, wgsize, 0, NULL, events);

			if (err == 0)
				return clblasSuccess;

		}
	}

	return clblasNotImplemented;
}

clblasStatus SGEMM_BRANCH_32(
	clblasTranspose transA,
	clblasTranspose transB,
	cl_uint M, cl_uint N, cl_uint K,
	float alpha,
	cl_mem A, cl_uint offA, cl_uint lda,
	cl_mem B, cl_uint offB, cl_uint ldb,
	float beta,
	cl_mem C, cl_uint offC, cl_uint ldc,
	cl_uint numCommandQueues,
	cl_command_queue *commandQueues,
	cl_uint numEventsInWaitList,
	const cl_event *eventWaitList,
	cl_event *events,
	bool &specialCaseHandled)
{
	const char *tileKernelSource = NULL;
	cl_kernel   tileClKernel = NULL;
	size_t tileKernelBinarySize = 0;
	cl_int err;


	const unsigned char *tileKernelBinary = NULL;

	clblasStatus status;

	if ((M * N < 1080 * 1080) && (M % 32 != 0 || N % 32 != 0) && (K%16==0))
	{
		// ((Mvalue - 1) / 32 + 1) * 16
		size_t GlobalX = ((M - 1) / 32 + 1) * 16;
		size_t GlobalY = ((N - 1) / 32 + 1) * 16;
		size_t gs[2] = { GlobalX, GlobalY };
		size_t wgsize[2] = { 16, 16 };

		if (transA == clblasNoTrans && transB == clblasNoTrans)
		{
			specialCaseHandled = true;
			tileKernelSource = sgemm_Col_NN_B1_MX032_NX032_KX16_BRANCH_src;
			tileClKernel = sgemm_Col_NN_B1_MX032_NX032_KX16_BRANCH_clKernel;
			tileKernelBinary = sgemm_Col_NN_B1_MX032_NX032_KX16_BRANCH_bin;
			tileKernelBinarySize = sgemm_Col_NN_B1_MX032_NX032_KX16_BRANCH_binSize;

			makeGemmKernel(&tileClKernel, commandQueues[0], tileKernelSource, User_srcBuildOptions, &tileKernelBinary, &tileKernelBinarySize, User_binBuildOptions);

			err = clSetKernelArg(tileClKernel, 0, sizeof(cl_mem), &A);
			CL_CHECK(err);
			err = clSetKernelArg(tileClKernel, 1, sizeof(cl_mem), &B);
			CL_CHECK(err);
			err = clSetKernelArg(tileClKernel, 2, sizeof(cl_mem), &C);
			CL_CHECK(err);
			err = clSetKernelArg(tileClKernel, 3, sizeof(cl_float), &alpha);
			CL_CHECK(err);
			err = clSetKernelArg(tileClKernel, 4, sizeof(cl_float), &beta);
			CL_CHECK(err);
			err = clSetKernelArg(tileClKernel, 5, sizeof(cl_uint), &M);
			CL_CHECK(err);
			err = clSetKernelArg(tileClKernel, 6, sizeof(cl_uint), &N);
			CL_CHECK(err);
			err = clSetKernelArg(tileClKernel, 7, sizeof(cl_uint), &K);
			CL_CHECK(err);
			err = clSetKernelArg(tileClKernel, 8, sizeof(cl_uint), &lda);
			CL_CHECK(err);
			err = clSetKernelArg(tileClKernel, 9, sizeof(cl_uint), &ldb);
			CL_CHECK(err);
			err = clSetKernelArg(tileClKernel, 10, sizeof(cl_uint), &ldc);
			CL_CHECK(err);
			err = clSetKernelArg(tileClKernel, 11, sizeof(cl_uint), &offA);
			CL_CHECK(err);
			err = clSetKernelArg(tileClKernel, 12, sizeof(cl_uint), &offB);
			CL_CHECK(err);
			err = clSetKernelArg(tileClKernel, 13, sizeof(cl_uint), &offC);
			CL_CHECK(err);

			err = clEnqueueNDRangeKernel(commandQueues[0], tileClKernel, 2, NULL,
				gs, wgsize, numEventsInWaitList, eventWaitList, &events[0]);

			if (err == 0)
				return clblasSuccess;
		}
		if (transA == clblasNoTrans && transB == clblasTrans)
		{
			specialCaseHandled = true;
			tileKernelSource = sgemm_Col_NT_B1_MX032_NX032_KX16_BRANCH_src;
			tileClKernel = sgemm_Col_NT_B1_MX032_NX032_KX16_BRANCH_clKernel;
			tileKernelBinary = sgemm_Col_NT_B1_MX032_NX032_KX16_BRANCH_bin;
			tileKernelBinarySize = sgemm_Col_NT_B1_MX032_NX032_KX16_BRANCH_binSize;

			makeGemmKernel(&tileClKernel, commandQueues[0], tileKernelSource, User_srcBuildOptions, &tileKernelBinary, &tileKernelBinarySize, User_binBuildOptions);

			err = clSetKernelArg(tileClKernel, 0, sizeof(cl_mem), &A);
			CL_CHECK(err);
			err = clSetKernelArg(tileClKernel, 1, sizeof(cl_mem), &B);
			CL_CHECK(err);
			err = clSetKernelArg(tileClKernel, 2, sizeof(cl_mem), &C);
			CL_CHECK(err);
			err = clSetKernelArg(tileClKernel, 3, sizeof(cl_float), &alpha);
			CL_CHECK(err);
			err = clSetKernelArg(tileClKernel, 4, sizeof(cl_float), &beta);
			CL_CHECK(err);
			err = clSetKernelArg(tileClKernel, 5, sizeof(cl_uint), &M);
			CL_CHECK(err);
			err = clSetKernelArg(tileClKernel, 6, sizeof(cl_uint), &N);
			CL_CHECK(err);
			err = clSetKernelArg(tileClKernel, 7, sizeof(cl_uint), &K);
			CL_CHECK(err);
			err = clSetKernelArg(tileClKernel, 8, sizeof(cl_uint), &lda);
			CL_CHECK(err);
			err = clSetKernelArg(tileClKernel, 9, sizeof(cl_uint), &ldb);
			CL_CHECK(err);
			err = clSetKernelArg(tileClKernel, 10, sizeof(cl_uint), &ldc);
			CL_CHECK(err);
			err = clSetKernelArg(tileClKernel, 11, sizeof(cl_uint), &offA);
			CL_CHECK(err);
			err = clSetKernelArg(tileClKernel, 12, sizeof(cl_uint), &offB);
			CL_CHECK(err);
			err = clSetKernelArg(tileClKernel, 13, sizeof(cl_uint), &offC);
			CL_CHECK(err);

			err = clEnqueueNDRangeKernel(commandQueues[0], tileClKernel, 2, NULL,
				gs, wgsize, numEventsInWaitList, eventWaitList, &events[0]);

			if (err == 0)
				return clblasSuccess;
		}
		if (transA == clblasTrans && transB == clblasNoTrans)
		{
			specialCaseHandled = true;
			tileKernelSource = sgemm_Col_TN_B1_MX032_NX032_KX16_BRANCH_src;
			tileClKernel = sgemm_Col_TN_B1_MX032_NX032_KX16_BRANCH_clKernel;
			tileKernelBinary = sgemm_Col_TN_B1_MX032_NX032_KX16_BRANCH_bin;
			tileKernelBinarySize = sgemm_Col_TN_B1_MX032_NX032_KX16_BRANCH_binSize;

			makeGemmKernel(&tileClKernel, commandQueues[0], tileKernelSource, User_srcBuildOptions, &tileKernelBinary, &tileKernelBinarySize, User_binBuildOptions);

			err = clSetKernelArg(tileClKernel, 0, sizeof(cl_mem), &A);
			CL_CHECK(err);
			err = clSetKernelArg(tileClKernel, 1, sizeof(cl_mem), &B);
			CL_CHECK(err);
			err = clSetKernelArg(tileClKernel, 2, sizeof(cl_mem), &C);
			CL_CHECK(err);
			err = clSetKernelArg(tileClKernel, 3, sizeof(cl_float), &alpha);
			CL_CHECK(err);
			err = clSetKernelArg(tileClKernel, 4, sizeof(cl_float), &beta);
			CL_CHECK(err);
			err = clSetKernelArg(tileClKernel, 5, sizeof(cl_uint), &M);
			CL_CHECK(err);
			err = clSetKernelArg(tileClKernel, 6, sizeof(cl_uint), &N);
			CL_CHECK(err);
			err = clSetKernelArg(tileClKernel, 7, sizeof(cl_uint), &K);
			CL_CHECK(err);
			err = clSetKernelArg(tileClKernel, 8, sizeof(cl_uint), &lda);
			CL_CHECK(err);
			err = clSetKernelArg(tileClKernel, 9, sizeof(cl_uint), &ldb);
			CL_CHECK(err);
			err = clSetKernelArg(tileClKernel, 10, sizeof(cl_uint), &ldc);
			CL_CHECK(err);
			err = clSetKernelArg(tileClKernel, 11, sizeof(cl_uint), &offA);
			CL_CHECK(err);
			err = clSetKernelArg(tileClKernel, 12, sizeof(cl_uint), &offB);
			CL_CHECK(err);
			err = clSetKernelArg(tileClKernel, 13, sizeof(cl_uint), &offC);
			CL_CHECK(err);

			err = clEnqueueNDRangeKernel(commandQueues[0], tileClKernel, 2, NULL,
				gs, wgsize, numEventsInWaitList, eventWaitList, &events[0]);

			if (err == 0)
				return clblasSuccess;
		}
	}

	return clblasNotImplemented;
}

clblasStatus DGEMM_BIG_MOD48(
	clblasTranspose transA,
	clblasTranspose transB,
	cl_uint M, cl_uint N, cl_uint K,
	double alpha,
	cl_mem A, cl_uint offA, cl_uint lda,
	cl_mem B, cl_uint offB, cl_uint ldb,
	double beta,
	cl_mem C, cl_uint offC, cl_uint ldc,
	cl_uint numCommandQueues,
	cl_command_queue *commandQueues,
	cl_uint numEventsInWaitList,
	const cl_event *eventWaitList,
	cl_event *events,
	bool &specialCaseHandled)
{
	const char *tileKernelSource = NULL;
	cl_kernel   tileClKernel = NULL;
	size_t tileKernelBinarySize = 0;
	cl_int err;


	const unsigned char *tileKernelBinary = NULL;

	clblasStatus status;
	//split the kernel calls to handle dgemm NT perf drop when matrix sizes are big
	if ((lda == ldb) && (lda >= 18000) && (lda <= 36000)) // between 18000 and 36000 for now
	{
		if (!((transA == clblasNoTrans) && (transB == clblasTrans)))
			return clblasNotImplemented;

		unsigned int M_split_factor;
		unsigned int N_split_factor;
		unsigned int K_split_factor;
		if ((M % 192 == 0) && (N % 192 == 0) && (K % 192 == 0) && (K > lda / 4))
		{
			M_split_factor = 4;
			N_split_factor = 4;
			K_split_factor = 4;
		}
		else if ((M % 96 == 0) && (N % 96 == 0) && (K % 96 == 0) && (K > lda / 4))
		{
			M_split_factor = 2;
			N_split_factor = 2;
			K_split_factor = 2;
		}
		else
		{
			return clblasNotImplemented;
		}

		tileKernelSource = dgemm_Col_NT_B1_MX048_NX048_KX08_src;
		tileClKernel = dgemm_Col_NT_B1_MX048_NX048_KX08_clKernel;
		tileKernelBinary = dgemm_Col_NT_B1_MX048_NX048_KX08_bin;
		tileKernelBinarySize = dgemm_Col_NT_B1_MX048_NX048_KX08_binSize;

		makeGemmKernel(&tileClKernel, commandQueues[0], tileKernelSource, User_srcBuildOptions, &tileKernelBinary, &tileKernelBinarySize, User_binBuildOptions);

		err = clSetKernelArg(tileClKernel, 0, sizeof(cl_mem), &A);
		CL_CHECK(err);
		err = clSetKernelArg(tileClKernel, 1, sizeof(cl_mem), &B);
		CL_CHECK(err);
		err = clSetKernelArg(tileClKernel, 2, sizeof(cl_mem), &C);
		CL_CHECK(err);
		err = clSetKernelArg(tileClKernel, 3, sizeof(cl_double), &alpha);
		CL_CHECK(err);
		err = clSetKernelArg(tileClKernel, 4, sizeof(cl_double), &beta);
		CL_CHECK(err);
		err = clSetKernelArg(tileClKernel, 5, sizeof(cl_uint), &M);
		CL_CHECK(err);
		err = clSetKernelArg(tileClKernel, 6, sizeof(cl_uint), &N);
		CL_CHECK(err);
		err = clSetKernelArg(tileClKernel, 7, sizeof(cl_uint), &K);
		CL_CHECK(err);
		err = clSetKernelArg(tileClKernel, 8, sizeof(cl_uint), &lda);
		CL_CHECK(err);
		err = clSetKernelArg(tileClKernel, 9, sizeof(cl_uint), &ldb);
		CL_CHECK(err);
		err = clSetKernelArg(tileClKernel, 10, sizeof(cl_uint), &ldc);
		CL_CHECK(err);
		err = clSetKernelArg(tileClKernel, 11, sizeof(cl_uint), &offA);
		CL_CHECK(err);
		err = clSetKernelArg(tileClKernel, 12, sizeof(cl_uint), &offB);
		CL_CHECK(err);
		err = clSetKernelArg(tileClKernel, 13, sizeof(cl_uint), &offC);
		CL_CHECK(err);

		status = GEMM_SPLIT_CALLS(
			tileClKernel, clblasColumnMajor,
			48, 8,
			M_split_factor,
			N_split_factor, K_split_factor,
			transA,
			transB,
			M, N, K,
			alpha,
			A, offA, lda,
			B, offB, ldb,
			beta,
			C, offC, ldc,
			numCommandQueues,
			commandQueues,
			numEventsInWaitList,
			eventWaitList,
			events);
		if (status == clblasSuccess)
			specialCaseHandled = true;

		return status;
	}


	return clblasNotImplemented;
}

template<>
clblasStatus
GemmSpecialCases<float>(clblasOrder order,
clblasTranspose transA,
clblasTranspose transB,
cl_uint M, cl_uint N, cl_uint K,
float alpha,
cl_mem A, cl_uint offA, cl_uint lda,
cl_mem B, cl_uint offB, cl_uint ldb,
float beta,
cl_mem C, cl_uint offC, cl_uint ldc,
cl_uint numCommandQueues,
cl_command_queue *commandQueues,
cl_uint numEventsInWaitList,
const cl_event *eventWaitList,
cl_event *events,
bool &specialCaseHandled)
{

	if (order == clblasRowMajor)
		return clblasNotImplemented;

	clblasStatus status;

	//handles big multiples of 1024
	status = SGEMM_mod1024(transA,
		transB,
		M, N, K,
		alpha,
		A, offA, lda,
		B, offB, ldb,
		beta,
		C, offC, ldc,
		numCommandQueues,
		commandQueues,
		numEventsInWaitList,
		eventWaitList,
		events,
		specialCaseHandled);
	if (specialCaseHandled)
		return status;

	//handles mod32 but not mod64
	status = SGEMM_SPLIT64_32(transA,
		transB,
		M, N, K,
		alpha,
		A, offA, lda,
		B, offB, ldb,
		beta,
		C, offC, ldc,
		numCommandQueues,
		commandQueues,
		numEventsInWaitList,
		eventWaitList,
		events,
		specialCaseHandled);
	if (specialCaseHandled)
		return status;

	//handles middle range sgemm (M*N<1080*1080) that are not mod32 (M%32!=0 || N%32!=0)
	//use 32x32 micro tile kernels with branch statement within kernels
	status = SGEMM_BRANCH_32(transA,
		transB,
		M, N, K,
		alpha,
		A, offA, lda,
		B, offB, ldb,
		beta,
		C, offC, ldc,
		numCommandQueues,
		commandQueues,
		numEventsInWaitList,
		eventWaitList,
		events,
		specialCaseHandled);
	if (specialCaseHandled)
		return status;


	return clblasNotImplemented;
}

template<>
clblasStatus
GemmSpecialCases<double>(clblasOrder order,
clblasTranspose transA,
clblasTranspose transB,
cl_uint M, cl_uint N, cl_uint K,
double alpha,
cl_mem A, cl_uint offA, cl_uint lda,
cl_mem B, cl_uint offB, cl_uint ldb,
double beta,
cl_mem C, cl_uint offC, cl_uint ldc,
cl_uint numCommandQueues,
cl_command_queue *commandQueues,
cl_uint numEventsInWaitList,
const cl_event *eventWaitList,
cl_event *events,
bool &specialCaseHandled)
{
	if (order == clblasRowMajor)
		return clblasNotImplemented;

	clblasStatus status;
	
	status = DGEMM_BIG_MOD48(transA,
		transB,
		M, N, K,
		alpha,
		A, offA, lda,
		B, offB, ldb,
		beta,
		C, offC, ldc,
		numCommandQueues,
		commandQueues,
		numEventsInWaitList,
		eventWaitList,
		events,
		specialCaseHandled);
		
	if (specialCaseHandled)
		return status;

	return clblasNotImplemented;
}

template<>
clblasStatus
GemmSpecialCases<FloatComplex>(clblasOrder order,
clblasTranspose transA,
clblasTranspose transB,
cl_uint M, cl_uint N, cl_uint K,
FloatComplex alpha,
cl_mem A, cl_uint offA, cl_uint lda,
cl_mem B, cl_uint offB, cl_uint ldb,
FloatComplex beta,
cl_mem C, cl_uint offC, cl_uint ldc,
cl_uint numCommandQueues,
cl_command_queue *commandQueues,
cl_uint numEventsInWaitList,
const cl_event *eventWaitList,
cl_event *events,
bool &specialCaseHandled)
{
	return clblasNotImplemented;
}

template<>
clblasStatus
GemmSpecialCases<DoubleComplex>(clblasOrder order,
clblasTranspose transA,
clblasTranspose transB,
cl_uint M, cl_uint N, cl_uint K,
DoubleComplex alpha,
cl_mem A, cl_uint offA, cl_uint lda,
cl_mem B, cl_uint offB, cl_uint ldb,
DoubleComplex beta,
cl_mem C, cl_uint offC, cl_uint ldc,
cl_uint numCommandQueues,
cl_command_queue *commandQueues,
cl_uint numEventsInWaitList,
const cl_event *eventWaitList,
cl_event *events,
bool &specialCaseHandled)
{
	return clblasNotImplemented;
}