bullet3

/******************************************************************************
 * Copyright 2010 Duane Merrill
 * 
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may ob3ain 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. 
 * 
 * 
 * 
 * 
 * AUTHORS' REQUEST: 
 * 
 * 		If you use|reference|benchmark this code, please cite our Technical 
 * 		Report (http://www.cs.virginia.edu/~dgm4d/papers/RadixSortTR.pdf):
 * 
 *		@TechReport{ Merrill:Sorting:2010,
 *        	author = "Duane Merrill and Andrew Grimshaw",
 *        	title = "Revisiting Sorting for GPGPU Stream Architectures",
 *        	year = "2010",
 *        	institution = "University of Virginia, Department of Computer Science",
 *        	address = "Charlottesville, VA, USA",
 *        	number = "CS2010-03"
 *		}
 * 
 * For more information, see our Google Code project site: 
 * http://code.google.com/p/back40computing/
 * 
 * Thanks!
 ******************************************************************************/

/******************************************************************************
 * Simple test driver program for *large-problem* radix sorting.
 *
 * Useful for demonstrating how to integrate radix sorting into 
 * your application 
 ******************************************************************************/

/******************************************************************************
 * Converted from CUDA to OpenCL/DirectCompute by Erwin Coumans
 ******************************************************************************/
#ifdef _WIN32
#pragma warning(disable : 4996)
#endif
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>
#include <float.h>
#include <algorithm>
#include <string>

//#include <iostream>
#include <sstream>
/**********************
*
*/

#include "Bullet3OpenCL/ParallelPrimitives/b3RadixSort32CL.h"
#include "Bullet3OpenCL/Initialize/b3OpenCLUtils.h"
#include "../btgui/Bullet3AppSupport/b3Clock.h"

cl_context g_cxMainContext;
cl_device_id g_device;
cl_command_queue g_cqCommandQueue;

/***********************
*
*/

bool g_verbose;
///Preferred OpenCL device/platform. When < 0 then no preference is used.
///Note that b3OpenCLUtils might still use the preference of using a platform vendor that matches the SDK vendor used to build the application.
///Preferred device/platform take priority over this platform-vendor match
int gPreferredDeviceId = -1;
int gPreferredPlatformId = -1;

/******************************************************************************
 * Routines
 ******************************************************************************/

/**
 * Keys-only sorting.  Uses the GPU to sort the specified vector of elements for the given 
 * number of iterations, displaying runtime information.
 *
 * @param[in] 		num_elements 
 * 		Size in elements of the vector to sort
 * @param[in] 		h_keys 
 * 		Vector of keys to sort 
 * @param[in] 		iterations  
 * 		Number of times to invoke the GPU sorting primitive
  * @param[in] 		cfg 
 * 		Config
 */
template <typename K>
void TimedSort(
	unsigned int num_elements,
	K *h_keys,
	unsigned int iterations)
{
	printf("Keys only, %d iterations, %d elements\n", iterations, num_elements);

	int max_elements = num_elements;
	b3AlignedObjectArray<unsigned int> hostData;
	hostData.resize(num_elements);
	for (int i = 0; i < num_elements; i++)
	{
		hostData[i] = h_keys[i];
	}

	b3RadixSort32CL sorter(g_cxMainContext, g_device, g_cqCommandQueue);

	b3OpenCLArray<unsigned int> gpuData(g_cxMainContext, g_cqCommandQueue);
	gpuData.copyFromHost(hostData);
	//sorter.executeHost(gpuData);
	sorter.execute(gpuData);

	b3AlignedObjectArray<unsigned int> hostDataSorted;
	gpuData.copyToHost(hostDataSorted);

	clFinish(g_cqCommandQueue);

	{
		//printf("Key-values, %d iterations, %d elements", iterations, num_elements);

		// Create sorting enactor

		// Perform the timed number of sorting iterations
		double elapsed = 0;
		float duration = 0;
		b3Clock watch;

		//warm-start
		gpuData.copyFromHost(hostData);
		clFinish(g_cqCommandQueue);
		sorter.execute(gpuData);

		watch.reset();

		for (int i = 0; i < iterations; i++)
		{
			// Move a fresh copy of the problem into device storage
			gpuData.copyFromHost(hostData);
			clFinish(g_cqCommandQueue);

			// Start GPU timing record
			double startMs = watch.getTimeMicroseconds() / 1e3;

			// Call the sorting API routine
			sorter.execute(gpuData);

			clFinish(g_cqCommandQueue);

			double stopMs = watch.getTimeMicroseconds() / 1e3;

			duration = stopMs - startMs;

			// End GPU timing record
			elapsed += (double)duration;
			printf("duration = %f\n", duration);
		}

		// Display timing information
		double avg_runtime = elapsed / iterations;
		//	double throughput = ((double) num_elements) / avg_runtime / 1000.0 / 1000.0;
		//   printf(", %f GPU ms, %f x10^9 elts/sec\n", 	avg_runtime,	throughput);
		double throughput = ((double)num_elements) / avg_runtime / 1000.0;
		printf(", %f GPU ms, %f x10^6 elts/sec\n", avg_runtime, throughput);

		gpuData.copyToHost(hostData);
		for (int i = 0; i < num_elements; i++)
		{
			h_keys[i] = hostData[i];
		}
	}
}

/**
 * Key-value sorting.  Uses the GPU to sort the specified vector of elements for the given 
 * number of iterations, displaying runtime information.
 *
 * @param[in] 		num_elements 
 * 		Size in elements of the vector to sort
 * @param[in] 		h_keys 
 * 		Vector of keys to sort 
 * @param[in,out] 	h_values  
 * 		Vector of values to sort 
 * @param[in] 		iterations  
 * 		Number of times to invoke the GPU sorting primitive
  * @param[in] 		cfg 
 * 		Config
 */
template <typename K, typename V>
void TimedSort(
	unsigned int num_elements,
	K *h_keys,
	V *h_values,
	unsigned int iterations)
{
	printf("Key-values, %d iterations, %d elements\n", iterations, num_elements);

	int max_elements = num_elements;
	b3AlignedObjectArray<b3SortData> hostData;
	hostData.resize(num_elements);
	for (int i = 0; i < num_elements; i++)
	{
		hostData[i].m_key = h_keys[i];
		hostData[i].m_value = h_values[i];
	}

	b3RadixSort32CL sorter(g_cxMainContext, g_device, g_cqCommandQueue);

	b3OpenCLArray<b3SortData> gpuData(g_cxMainContext, g_cqCommandQueue);
	gpuData.copyFromHost(hostData);
	//sorter.executeHost(gpuData);
	sorter.execute(gpuData);

	b3AlignedObjectArray<b3SortData> hostDataSorted;
	gpuData.copyToHost(hostDataSorted);
#if 0
    for (int i=0;i<num_elements;i++)
	{
		printf("hostData[%d].m_key = %d\n",i, hostDataSorted[i].m_key);
        printf("hostData[%d].m_value = %d\n",i,hostDataSorted[i].m_value);
	}
#endif

	clFinish(g_cqCommandQueue);

	{
		//printf("Key-values, %d iterations, %d elements", iterations, num_elements);

		// Create sorting enactor

		// Perform the timed number of sorting iterations
		double elapsed = 0;
		float duration = 0;
		b3Clock watch;

		//warm-start
		gpuData.copyFromHost(hostData);
		sorter.execute(gpuData);
		clFinish(g_cqCommandQueue);

		watch.reset();

		for (int i = 0; i < iterations; i++)
		{
			// Move a fresh copy of the problem into device storage
			gpuData.copyFromHost(hostData);
			clFinish(g_cqCommandQueue);

			// Start GPU timing record
			double startMs = watch.getTimeMicroseconds() / 1e3;

			// Call the sorting API routine
			sorter.execute(gpuData);
			clFinish(g_cqCommandQueue);

			double stopMs = watch.getTimeMicroseconds() / 1e3;

			duration = stopMs - startMs;

			// End GPU timing record
			elapsed += (double)duration;
			printf("duration = %f\n", duration);
		}

		// Display timing information
		double avg_runtime = elapsed / iterations;
		//	double throughput = ((double) num_elements) / avg_runtime / 1000.0 / 1000.0;
		//   printf(", %f GPU ms, %f x10^9 elts/sec\n", 	avg_runtime,	throughput);
		double throughput = ((double)num_elements) / avg_runtime / 1000.0;
		printf(", %f GPU ms, %f x10^6 elts/sec\n", avg_runtime, throughput);

		gpuData.copyToHost(hostData);
		for (int i = 0; i < num_elements; i++)
		{
			h_keys[i] = hostData[i].m_key;
			h_values[i] = hostData[i].m_value;
		}
	}
}

/**
 * Generates random 32-bit keys.
 * 
 * We always take the second-order byte from rand() because the higher-order 
 * bits returned by rand() are commonly considered more uniformly distributed
 * than the lower-order bits.
 * 
 * We can decrease the entropy level of keys by adopting the technique 
 * of Thearling and Smith in which keys are computed from the bitwise AND of 
 * multiple random samples: 
 * 
 * entropy_reduction	| Effectively-unique bits per key
 * -----------------------------------------------------
 * -1					| 0
 * 0					| 32
 * 1					| 25.95
 * 2					| 17.41
 * 3					| 10.78
 * 4					| 6.42
 * ...					| ...
 * 
 */
template <typename K>
void RandomBits(K &key, int entropy_reduction = 0, int lower_key_bits = sizeof(K) * 8)
{
	const unsigned int NUM_UCHARS = (sizeof(K) + sizeof(unsigned char) - 1) / sizeof(unsigned char);
	unsigned char key_bits[NUM_UCHARS];

	do
	{
		for (int j = 0; j < NUM_UCHARS; j++)
		{
			unsigned char quarterword = 0xff;
			for (int i = 0; i <= entropy_reduction; i++)
			{
				quarterword &= (rand() >> 7);
			}
			key_bits[j] = quarterword;
		}

		if (lower_key_bits < sizeof(K) * 8)
		{
			unsigned long long base = 0;
			memcpy(&base, key_bits, sizeof(K));
			base &= (1 << lower_key_bits) - 1;
			memcpy(key_bits, &base, sizeof(K));
		}

		memcpy(&key, key_bits, sizeof(K));

	} while (key != key);  // avoids NaNs when generating random floating point numbers
}

/******************************************************************************
 * Templated routines for printing keys/values to the console 
 ******************************************************************************/

template <typename T>
void PrintValue(T val)
{
	printf("%d", val);
}

template <>
void PrintValue<float>(float val)
{
	printf("%f", val);
}

template <>
void PrintValue<double>(double val)
{
	printf("%f", val);
}

template <>
void PrintValue<unsigned char>(unsigned char val)
{
	printf("%u", val);
}

template <>
void PrintValue<unsigned short>(unsigned short val)
{
	printf("%u", val);
}

template <>
void PrintValue<unsigned int>(unsigned int val)
{
	printf("%u", val);
}

template <>
void PrintValue<long>(long val)
{
	printf("%ld", val);
}

template <>
void PrintValue<unsigned long>(unsigned long val)
{
	printf("%lu", val);
}

template <>
void PrintValue<long long>(long long val)
{
	printf("%lld", val);
}

template <>
void PrintValue<unsigned long long>(unsigned long long val)
{
	printf("%llu", val);
}

/**
 * Compares the equivalence of two arrays
 */
template <typename T, typename SizeT>
int CompareResults(T *computed, T *reference, SizeT len, bool verbose = true)
{
	printf("\n");
	for (SizeT i = 0; i < len; i++)
	{
		if (computed[i] != reference[i])
		{
			printf("INCORRECT: [%lu]: ", (unsigned long)i);
			PrintValue<T>(computed[i]);
			printf(" != ");
			PrintValue<T>(reference[i]);

			if (verbose)
			{
				printf("\nresult[...");
				for (size_t j = (i >= 5) ? i - 5 : 0; (j < i + 5) && (j < len); j++)
				{
					PrintValue<T>(computed[j]);
					printf(", ");
				}
				printf("...]");
				printf("\nreference[...");
				for (size_t j = (i >= 5) ? i - 5 : 0; (j < i + 5) && (j < len); j++)
				{
					PrintValue<T>(reference[j]);
					printf(", ");
				}
				printf("...]");
			}

			return 1;
		}
	}

	printf("CORRECT\n");
	return 0;
}

/**
 * Creates an example sorting problem whose keys is a vector of the specified 
 * number of K elements, values of V elements, and then dispatches the problem 
 * to the GPU for the given number of iterations, displaying runtime information.
 *
 * @param[in] 		iterations  
 * 		Number of times to invoke the GPU sorting primitive
 * @param[in] 		num_elements 
 * 		Size in elements of the vector to sort
 * @param[in] 		cfg 
 * 		Config
 */
template <typename K, typename V>
void TestSort(
	unsigned int iterations,
	int num_elements,
	bool keys_only)
{
	// Allocate the sorting problem on the host and fill the keys with random bytes

	K *h_keys = NULL;
	K *h_reference_keys = NULL;
	V *h_values = NULL;
	h_keys = (K *)malloc(num_elements * sizeof(K));
	h_reference_keys = (K *)malloc(num_elements * sizeof(K));
	if (!keys_only) h_values = (V *)malloc(num_elements * sizeof(V));

	// Use random bits
	for (unsigned int i = 0; i < num_elements; ++i)
	{
		RandomBits<K>(h_keys[i], 0);
		//h_keys[i] = num_elements-i;
		//h_keys[i] = 0xffffffffu-i;
		if (!keys_only)
			h_values[i] = h_keys[i];  //0xffffffffu-i;

		h_reference_keys[i] = h_keys[i];
	}

	// Run the timing test
	if (keys_only)
	{
		TimedSort<K>(num_elements, h_keys, iterations);
	}
	else
	{
		TimedSort<K, V>(num_elements, h_keys, h_values, iterations);
	}

	//	cudaThreadSynchronize();

	// Display sorted key data
	if (g_verbose)
	{
		printf("\n\nKeys:\n");
		for (int i = 0; i < num_elements; i++)
		{
			PrintValue<K>(h_keys[i]);
			printf(", ");
		}
		printf("\n\n");
	}

	// Verify solution
	std::sort(h_reference_keys, h_reference_keys + num_elements);
	CompareResults<K>(h_keys, h_reference_keys, num_elements, true);
	printf("\n");
	fflush(stdout);

	// Free our allocated host memory
	if (h_keys != NULL) free(h_keys);
	if (h_values != NULL) free(h_values);
}

/**
 * Displays the commandline usage for this tool
 */
void Usage()
{
	printf("\ntest_large_problem_sorting [--device=<device index>] [--v] [--i=<num-iterations>] [--n=<num-elements>] [--key-values] [--deviceId=<int>] [--platformId=<int>]\n");
	printf("\n");
	printf("\t--v\tDisplays sorted results to the console.\n");
	printf("\n");
	printf("\t--i\tPerforms the sorting operation <num-iterations> times\n");
	printf("\t\t\ton the device. Re-copies original input each time. Default = 1\n");
	printf("\n");
	printf("\t--n\tThe number of elements to comprise the sample problem\n");
	printf("\t\t\tDefault = 512\n");
	printf("\n");
	printf("\t--key-values\tSpecifies that keys are accommodated by value pairings\n");
	printf("\n");
}

/******************************************************************************
 * Command-line parsing
 ******************************************************************************/
#include <map>
#include <algorithm>
#include <string>

class b3CommandLineArgs
{
protected:
	std::map<std::string, std::string> pairs;

public:
	// Constructor
	b3CommandLineArgs(int argc, char **argv)
	{
		using namespace std;

		for (int i = 1; i < argc; i++)
		{
			string arg = argv[i];

			if ((arg[0] != '-') || (arg[1] != '-'))
			{
				continue;
			}

			string::size_type pos;
			string key, val;
			if ((pos = arg.find('=')) == string::npos)
			{
				key = string(arg, 2, arg.length() - 2);
				val = "";
			}
			else
			{
				key = string(arg, 2, pos - 2);
				val = string(arg, pos + 1, arg.length() - 1);
			}
			pairs[key] = val;
		}
	}

	bool CheckCmdLineFlag(const char *arg_name)
	{
		using namespace std;
		map<string, string>::iterator itr;
		if ((itr = pairs.find(arg_name)) != pairs.end())
		{
			return true;
		}
		return false;
	}

	template <typename T>
	void GetCmdLineArgument(const char *arg_name, T &val);

	int ParsedArgc()
	{
		return pairs.size();
	}
};

template <typename T>
void b3CommandLineArgs::GetCmdLineArgument(const char *arg_name, T &val)
{
	using namespace std;
	map<string, string>::iterator itr;
	if ((itr = pairs.find(arg_name)) != pairs.end())
	{
		istringstream strstream(itr->second);
		strstream >> val;
	}
}

template <>
void b3CommandLineArgs::GetCmdLineArgument<char *>(const char *arg_name, char *&val)
{
	using namespace std;
	map<string, string>::iterator itr;
	if ((itr = pairs.find(arg_name)) != pairs.end())
	{
		string s = itr->second;
		val = (char *)malloc(sizeof(char) * (s.length() + 1));
		strcpy(val, s.c_str());
	}
	else
	{
		val = NULL;
	}
}

/******************************************************************************
 * Main
 ******************************************************************************/

extern bool gDebugSkipLoadingBinary;

void myprintf(const char *msg)
{
	(void *)msg;
}
int main(int argc, char **argv)
{
	//gDebugSkipLoadingBinary = true;

	//	b3SetCustomPrintfFunc(myprintf);

	cl_int ciErrNum;
	b3CommandLineArgs args(argc, argv);

	args.GetCmdLineArgument("deviceId", gPreferredDeviceId);
	args.GetCmdLineArgument("platformId", gPreferredPlatformId);

	b3Printf("Initialize OpenCL using b3OpenCLUtils_createContextFromType\n");
	cl_platform_id platformId;
	//	g_cxMainContext = b3OpenCLUtils_createContextFromType(CL_DEVICE_TYPE_ALL, &ciErrNum, 0, 0,gPreferredDeviceId,gPreferredPlatformId,&platformId);
	g_cxMainContext = b3OpenCLUtils_createContextFromType(CL_DEVICE_TYPE_GPU, &ciErrNum, 0, 0, gPreferredDeviceId, gPreferredPlatformId, &platformId);
	//g_cxMainContext = b3OpenCLUtils_createContextFromType(CL_DEVICE_TYPE_CPU, &ciErrNum, 0, 0,gPreferredDeviceId,gPreferredPlatformId,&platformId);

	oclCHECKERROR(ciErrNum, CL_SUCCESS);

	int numDev = b3OpenCLUtils_getNumDevices(g_cxMainContext);

	if (!numDev)
	{
		b3Error("error: no OpenCL devices\n");
		exit(0);
	}
	int devId = 0;
	g_device = b3OpenCLUtils_getDevice(g_cxMainContext, devId);
	b3OpenCLUtils_printDeviceInfo(g_device);
	// create a command-queue
	g_cqCommandQueue = clCreateCommandQueue(g_cxMainContext, g_device, 0, &ciErrNum);
	oclCHECKERROR(ciErrNum, CL_SUCCESS);

	//srand(time(NULL));
	srand(0);  // presently deterministic

	unsigned int num_elements = 8 * 1024 * 1024;  //4*1024*1024;//4*1024*1024;//257;//8*524288;//2048;//512;//524288;
	unsigned int iterations = 10;
	bool keys_only = true;

	//
	// Check command line arguments
	//

	if (args.CheckCmdLineFlag("help"))
	{
		Usage();
		return 0;
	}

	args.GetCmdLineArgument("i", iterations);
	args.GetCmdLineArgument("n", num_elements);

	keys_only = !args.CheckCmdLineFlag("key-values");
	g_verbose = args.CheckCmdLineFlag("v");

	TestSort<unsigned int, unsigned int>(
		iterations,
		num_elements,
		keys_only);
}