cpucomputation.cpp

PDE simulator on GPU.

// CPUcomputation.cpp: implementation of the CCPUcomputation class.
//
//////////////////////////////////////////////////////////////////////

#include "stdafx.h"
#include "shadow.h"
#include <math.h>
#include "CPUcomputation.h"

#ifdef _DEBUG
#undef THIS_FILE
static char THIS_FILE[]=__FILE__;
#define new DEBUG_NEW
#endif


//////////////////////////////////////////////////////////////////////////////
//// Operations on CPU VS on GPU
//////////////////////////////////////////////////////////////////////////////
void	on_cpu_clVecOp(OperationOnGPUEnum op,  float a, float b, float* x, float* y, float* result, int Dim_vector)
{
	int i;
	switch (op)
	{
	case CL_VECT_ADD_VECT_RECT:
		//////////////////////////////////////////////////////////////////
		////	Vector + Vector
		for(i=0;i<Dim_vector;i++)
		{
			result[i] = x[i]*a + y[i]*b;
		}
		break;
	case CL_VECT_SUB_VECT_RECT:
		/////////////////////////////////////////////////
		////	Vector - Vector
		for(i=0;i<Dim_vector;i++)
		{
			result[i] = x[i]*a - y[i]*b;
		}
		break;
	case CL_VECT_MULT_VECT_RECT:
		/////////////////////////////////////////////////
		////	Vector * Vector
		for(i=0;i<Dim_vector;i++)
		{
			result[i] = x[i]*a * y[i]*b;
		}
		break;
	}

	return;
}

void	on_cpu_clMatVec(float* A,  float* x, float* result, int Dim_Vector)
{
	for(int i=0; i<Dim_Vector; i++)
	{
		result[i] = 0;
		for(int j=0; j<Dim_Vector; j++)
		{
			result[i] += A[i*Dim_Vector + j] * x[j];	
		}
	}

	return;
}

float	on_cpu_clVecReduce(OperationOnGPUEnum cmb, float* x, int Dim_Vector)
{
	int i;
	float result = 0;
	switch(cmb)
	{		
	case CL_VECT_NORMAL_RECT:////sum the elements of Vector x
		/////////////////////////////////////////////////
		////	sum the elements of Vector x
		result = 0.0;
		for(i=0; i<Dim_Vector; i++)
		{
			result += x[i];
		}
		break;
	case CL_VECT_MAX_RECT:////retrieve the max element of Vector x
		/////////////////////////////////////////////////
		////	the max element of Vector x
		result = x[0];
		for(i=1; i<Dim_Vector; i++)
		{
			if (x[i] > result)
			{
				result = x[i];
			}
		}
		break;
	case CL_VECT_MIN_RECT:////retrieve the min element of Vector x
		/////////////////////////////////////////////////
		////	the min element of Vector x
		result = x[0];
		for(i=1; i<Dim_Vector; i++)
		{
			if (x[i] > result)
			{
				result = x[i];
			}
		}
		break;
	default:
		AfxMessageBox("Error Operation Type in clVecReduceFP!");
		exit(0);
	}

	return result;
}

void	on_cpu_clMatMat(float* A, float* B, float* result, int Dim_Vector)
{
	/////////////////////////////////////////////////
	////	Matrix * Matrix
	/////////////////////////////////////////////////
	////	first I transform A and B to a 2D matrix 
	////	why I do this is that 
	////	I think in loops, many mul operation should be used 
	////	to locate the correct element 
	////	while after transformation, these locating mul can be avoid
	int i,j,k;
	for(i=0;i<Dim_Vector;i++)
	{
		for(j=0; j<Dim_Vector; j++)
		{
			result[i*Dim_Vector + j] = 0;
			for (k=0; k<Dim_Vector; k++)
			{
				result[i*Dim_Vector + j] += A[i*Dim_Vector + k] * B[k*Dim_Vector + j];
			}
		}
	}

/*
	float **tempA, **tempB, **tempresult;
	tempA = new float* [Dim_Vector];
	tempB = new float* [Dim_Vector];
	tempresult = new float* [Dim_Vector];

	/////////////////////////////////////////////////
	////	transform A to 2D matrix
	for (i=0; i<Dim_Vector; i++)
	{
		tempA[i] = new float [Dim_Vector];
		for (j=0; j<Dim_Vector; j++)
		{
			tempA[i][j] = A[Dim_Vector*i + j];
		}

		tempresult[i] = new float [Dim_Vector];
	}

	/////////////////////////////////////////////////
	////	transform B to 2D matrix
	for (i=0; i<Dim_Vector; i++)
	{
		tempB[i] = new float [Dim_Vector];
		for ( j=0; j<Dim_Vector; j++)
		{
			tempB[i][j] = B[Dim_Vector*i + j];
		}
	}

	for(i=0;i<Dim_Vector;i++)
	{
		for(j=0; j<Dim_Vector; j++)
		{
			tempresult[i][j] = 0;
			for (k=0; k<Dim_Vector; k++)
			{
				tempresult[i][j] += tempA[i][k] * tempB[k][j];
			}
		}
	}

	for(i=0;i<Dim_Vector;i++)
	{
		delete[] tempA[i];
		delete[] tempB[i];
		delete[] tempresult[i];
	}
	delete[] tempA;
	delete[] tempB;
	delete[] tempresult;
*/
	return;
}

float	on_cpu_clInnerProduct(float* x, float* y, int Dim_Vector)
{
	/////////////////////////////////////////////////
	////	Vector . Vector
	float result = 0;

	for(int i=0; i<Dim_Vector; i++)
	{	
		result += x[i]*y[i];
		
	}
	return result;
}



//////////////////////////////////////////////////////////////////////
/// Total Object
/// Linear Equations Solver on CPU
/// Ax = b ----> solve x
/// Matrix----M[Dim][Dim]
/// Output Vector---x[Dim]
/// All the following methods are Iterative methods
/// Jacobi Solver
//////////////////////////////////////////////////////////////////////
void JacobiSolverOnCPU(float *M, float *x, float *b, int Dim)
{
	///// Previous x vector
	float *prev_x = new float[Dim];
	///// Tolerance of Iterations Error
	float ToL = 1.0e-5;
	///// Maximum number of Iterations 
	int MaxIterNum  = 1000;
	int i,j,k;
	float Tolerance, sum;
	for(k=0; k<MaxIterNum; k++)
	{
		Tolerance = 0.0f;
		memcpy(prev_x, x, sizeof(float)*Dim);
		for(i=0; i<Dim; i++)
		{
			///////
			sum = 0.0f;
			for(j=0; j<i; j++)
				sum += M[i*Dim+j]*prev_x[j];
			for(j=i+1; j<Dim; j++)
				sum += M[i*Dim+j]*prev_x[j];

			x[i] = (b[i] - sum)/M[i*Dim+i];

			Tolerance += (x[i] - prev_x[i])*(x[i] - prev_x[i]);

		}
		TRACE("Tolerance = %f\n",Tolerance);
		if(Tolerance < ToL)
		{
			break;
		}
	}
	TRACE("Max Iterations Num = %d\n", k);
	delete[] prev_x;

	return;
}

//////////////////////////////////////////////////////////////////////
/// Gauss-Seidel Solver
//////////////////////////////////////////////////////////////////////
void GaussSeidelSolverOnCPU(float *M, float *x, float *b, int Dim)
{
	///// Tolerance of Iterations Error
	float ToL = 1.0e-5;
	///// Maximum number of Iterations 
	int MaxIterNum  = 1000;
	int i,j,k;
	float Tolerance, sum, prev;
	for(k=0; k<MaxIterNum; k++)
	{
		Tolerance = 0.0f;
		for(i=0; i<Dim; i++)
		{
			///////
			prev = x[i];
			sum = 0.0f;
			for(j=0; j<i; j++)
				sum += M[i*Dim+j]*x[j];
			for(j=i+1; j<Dim; j++)
				sum += M[i*Dim+j]*x[j];
			x[i] = (b[i] - sum)/M[i*Dim+i];
			Tolerance += (x[i] - prev)*(x[i] - prev);
		}
		TRACE("Tolerance = %f\n",Tolerance);
		if(Tolerance < ToL)
		{
			break;
		}
	}
	TRACE("Max Iterations Num = %d\n", k);

	return;
}

//////////////////////////////////////////////////////////////////////
/// Successive Over-Relaxation Solover
//////////////////////////////////////////////////////////////////////
void SucOverRelaxSolverOnCPU(float *M, float *x, float *b, int Dim)
{
	///// Previous x vector
	float *prev_x = new float[Dim];
	///// Tolerance of Iterations Error
	float ToL = 1.0e-5;
	///// Relaxation Factor
	float omiga = 0.5;
	///// Maximum number of Iterations 
	int MaxIterNum  = 1000;
	int i,j,k;
	float Tolerance, sum;
	for(k=0; k<MaxIterNum; k++)
	{
		Tolerance = 0.0f;
		memcpy(prev_x, x, sizeof(float)*Dim);
		for(i=0; i<Dim; i++)
		{
			///////
			sum = 0.0f;
			for(j=0; j<Dim; j++)
				sum += M[i*Dim+j]*prev_x[j];

			x[i] = x[i] + omiga*(b[i] - sum)/M[i*Dim+i];

			Tolerance += (x[i] - prev_x[i])*(x[i] - prev_x[i]);

		}
		TRACE("Tolerance = %f\n",Tolerance);
		if(Tolerance < ToL)
		{
			break;
		}
	}
	TRACE("Max Iterations Num = %d\n", k);
	delete[] prev_x;

	return;
}

//////////////////////////////////////////////////////////////////////
/// Steepest Descent Solover
/// Matrix M must be positive definite
//////////////////////////////////////////////////////////////////////
void SteepestDescentSolverOnCPU(float *M, float *x, float *b, int Dim)
{
	///// the search direction vector for x
	float *direction = new float[Dim];

	///// Tolerance of Iterations Error
	float ToL = 1.0e-5;
	///// Maximum number of Iterations 
	int MaxIterNum  = 1000;
	int i,j,k;
	float Tolerance, sum, tk, tt;
	for(k=0; k<MaxIterNum; k++)
	{
		///// calculate the search direction
		for(i=0; i<Dim; i++)
		{
			//// Calculate the residual vector
			sum = 0.0f;
			for(j=0;j<Dim;j++)
				sum += M[i*Dim+j]*x[j];
			direction[i] = b[i] - sum;
		}
		
		////<v,Av>--->tk
		////<v,v>---->tt
		tt = tk = 0;
		for(i=0; i<Dim; i++)
		{
			/// Av
			sum = 0.0f;
			for(j=0; j<Dim; j++)
			{
				sum += M[i*Dim+j]*direction[j];
			}
			tk += sum*direction[i];
			tt += direction[i]*direction[i];
		}
		//// tk = <v,v>/<v,Av>
		tk = tt/tk;

		Tolerance = 0.0f;
		//// x(k) = x(k-1) + tk*v
		for(i=0; i<Dim; i++)
		{
			/////// Record the previous one, then sum the error
			tt = x[i];			
			x[i] = tt + tk*direction[i];

			Tolerance += (x[i] - tt)*(x[i] - tt);
		}

		TRACE("Tolerance = %f\n",Tolerance);
		if(Tolerance < ToL)
		{
			break;
		}
	}
	TRACE("Max Iterations Num = %d\n", k);
	delete[] direction;

	return;

}

//////////////////////////////////////////////////////////////////////
/// Unpreconditioned Conjugate Gradient Solover
/// Matrix M must be positive definite
//////////////////////////////////////////////////////////////////////
void ConjugateGradientSolverOnCPU(float *M, float *x, float *b, int Dim)
{
	///// the search direction vector for x
	float *direction = new float[Dim];
	float *residua = new float[Dim];

	///// Tolerance of Iterations Error
	float ToL = 1.0e-5;
	///// Maximum number of Iterations 
	int MaxIterNum  = 1000;
	int i,j,k;
	float Tolerance, sum, tk, tt, rr, sk;

	//// calculate r0,v1
	for(i=0;i<Dim;i++)
	{
		sum = 0.0f;
		for(j=0; j<Dim; j++)
			sum += M[i*Dim+j]*x[j];
		residua[i]   = b[i] - sum;
		direction[i] = residua[i];
	}

	for(k=0; k<MaxIterNum; k++)
	{
		///// calculcate tk
		rr = tk = 0.0f;
		for(i=0; i<Dim; i++)
		{
			//// <r,r>
			rr += residua[i]*residua[i];

			//// <v,Av>
			sum = 0.0f;
			for(j=0; j<Dim; j++)
				sum += M[i*Dim+j]*direction[j];
			 tk += sum*direction[i];
		}
		tk = rr/tk;

		Tolerance = 0.0f;
		sk = 0.0f;
		//// x(k) = x(k-1) + tk*v
		for(i=0; i<Dim; i++)
		{
			/////// Record the previous one, then sum the error
			tt = x[i];			
			x[i] = tt + tk*direction[i];

			//// <v,Av>
			sum = 0.0f;
			for(j=0; j<Dim; j++)
				sum += M[i*Dim+j]*direction[j];
			residua[i] = residua[i] - tk*sum;
			
			sk += residua[i]*residua[i];

			Tolerance += (x[i] - tt)*(x[i] - tt);
		}
		
		//// recalculate direction v
		sk = sk/rr;
		for(i=0; i<Dim; i++)
		{
			direction[i] = residua[i] + sk*direction[i];
		}
		TRACE("Tolerance = %f\n",Tolerance);
		if(Tolerance < ToL)
		{
			break;
		}
	}
	TRACE("Max Iterations Num = %d\n", k);
	delete[] direction;
	delete[] residua;

	return;

}