gpgpu/Linear/Jacobi.cpp

#include "LinearTests.h"
#include <iostream>

void Jacobi::generateLinEq()
{
	Jx_c[0] = new float[Jn];
	Jx_c[1] = new float[Jn];
	Jx_g[0] = new float[Jn];
	Jx_g[1] = new float[Jn];

	for (int i = 0; i < Jn; ++i) {
		Jx_c[0][i] = 0.0f;
		Jx_c[1][i] = 0.0f;
		Jx_g[0][i] = 0.0f;
		Jx_g[1][i] = 0.0f;
	}

	JA_c = new float[Jn * Jn];
	JA_g = new float[Jn * Jn];
	for (int i = 0; i < Jn; ++i) {
		for (int j = 0; j < Jn; ++j) {
			float v = 0.0f;
			if (i == j) {
				v = 0.5f;
			}
			JA_c[i + j * Jn] = v;
			JA_g[i + j * Jn] = v;
		}
	}

	Jb_c = new float[Jn];
	Jb_g = new float[Jn];
	for (int i = 0; i < Jn; ++i) {
		Jb_c[i] = 1.0f;
		Jb_g[i] = 1.0f;
	}
}

void Jacobi::cpuScalarMV(int n, int m, float* y, const float* A, const float* x, const float* b)
{
	for (int i = 0; i < n; ++i) {
		float yi = b[i];
		for (int j = 0; j < m; ++j) {
			yi += A[i * m + j] * x[j];
		}
		y[i] = yi;
	}
}

void Jacobi::printMatrix(int n, int m, float* A)
{
	for (int i = 0; i < n; ++i) {
		for (int j = 0; j < m; ++j) {
			std::cout << A[j + i * n];
			if (j < m - 1) std::cout << ", ";
		}
		std::cout << std::endl;
	}
}

MatrixVectorMultiplier* Jacobi::MethodFactory(MVType type, cl::Context* context, cl::CommandQueue* queue, cl::Program* program)
{
	if (type == MVType::SimpleMV) {
		return new Simple(context, queue, program);
	}
	else if (type == MVType::ReduceMV) {
		return new Reduce(context, queue, program);
	}
	else if (type == MVType::LargeMV) {
		return new Large(context, queue, program);
	}
	else {
		return NULL;
	}
}

Jacobi::Jacobi(int n, MVType _type)
{
	type = _type;
	Jn = n;
	generateLinEq();
}

void Jacobi::collect_results(cl::CommandQueue* queue)
{
}

void Jacobi::gpu_compute(cl::Context* context, cl::CommandQueue* queue, cl::Program* program, cl::Event* Event)
{
	MatrixVectorMultiplier* MVMultiplier = MethodFactory(type, context, queue, program);

	cl_ulong time = 0;

	if (MVMultiplier != NULL) {
		int inputBuffer = 0;
		const int iterations = 20;
		for (int i = 0; i < iterations; ++i) {
			time += MVMultiplier->dewIt(Jn, Jn, Jx_g[(inputBuffer + 1) % 2], JA_g, Jx_g[inputBuffer], Jb_g);
			inputBuffer = (inputBuffer + 1) % 2;
			//printMatrix(1, Jn, Jx_g[inputBuffer]);
		}
	}
	else {
		std::cout << "Invalid factory parameter" << std::endl;
		exit(-1);
	}
	std::cout << "Jacobi>" << MVMultiplier->getName() << "@" << Jn  << ": " << time / 1e+06 << " ms" << std::endl;
	
}

void Jacobi::cpu_compute()
{
	int inputBuffer = 0;
	const int iterations = 20;
	for (int i = 0; i < iterations; ++i) {
		cpuScalarMV(Jn, Jn, Jx_c[(inputBuffer + 1) % 2], JA_c, Jx_c[inputBuffer], Jb_c);
		//printMatrix(1, Jn, Jx_c[inputBuffer + 1]);
		inputBuffer = (inputBuffer + 1) % 2;
	}

}

bool Jacobi::validate_results()
{
	bool result = true;
	// Actual validation
	//printMatrix(Jn, Jn, JA_c);
	//printMatrix(1, Jn, Jx_c[0]);
	//printMatrix(1, Jn, Jx_c[1]);

	// Cleanup
	if (Jx_c[0] == 0) delete[] Jx_c[0];
	if (Jx_c[1] == 0) delete[] Jx_c[1];
	if (Jx_g[0] == 0) delete[] Jx_g[0];
	if (Jx_g[1] == 0) delete[] Jx_g[1];
	
	if (JA_c == 0) delete[] JA_c;
	if (Jb_c == 0) delete[] Jb_c;
	if (JA_g == 0) delete[] JA_g;
	if (Jb_g == 0) delete[] Jb_g;

	return result;
}

std::string Jacobi::description()
{
	return std::string();
}