testequationssolver.cpp

PDE simulator on GPU.

// TestEquationsSlover.cpp: implementation of the solver
//
//////////////////////////////////////////////////////////////////////

#include "stdafx.h"
#include "shadow.h"
#include "TestEquationsSolver.h"
#include "LinearEquationsOnGPU.h"
#include "CPUcomputation.h"

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

//////////////////////////////////////////////////////////////////////
// Test Complex Linear Equations Iterative
//////////////////////////////////////////////////////////////////////
void PerformanceIterativeSolverTestOnGPU()
{
	int i,j;
	int Dim = 4;
	float M[4*4];
/*
	float **M = new float*[Dim];
	for(i=0;i<Dim;i++)
	{
		M[i] = new float[Dim];
	}
*/
    float *x = new float[Dim];
	float *b = new float[Dim];
	float *result = new float[Dim];

	////////////////////////
/*
	for(i=0;i<Dim;i++)
	{
		for(j=0;j<Dim;j++)
			M[i][j] = (float)rand()/RAND_MAX;
		b[i] = 5.0*(float)rand()/RAND_MAX;
		x[i] = (float)rand()/RAND_MAX;
	}
*/
	M[0] = 1;
	M[1] = 0;
	M[2] = 0;
	M[3] = 0;

	M[4] = 0;
	M[5] = 1;
	M[6] = 0;
	M[7] = 0;

	M[8] = 0;
	M[9] = 0;
	M[10] = 1;
	M[11] = 0;

	M[12] = 0;
	M[13] = 0;
	M[14] = 0;
	M[15] = 1;

	b[0] = 0.05;
	b[1] = 0.10;
	b[2] = 0.15;
	b[3] = 0.20;

	x[0] = 1;
	x[1] = 1;
	x[2] = 1;
	x[3] = 1;

/*
	/////////////////	On CPU	/////////////////
	////	GaussSeidelSolverOnCPU
	QueryPerformanceCounter(&startcounter);
	GaussSeidelSolverOnCPU(&M[0], x, b, Dim);
	QueryPerformanceCounter(&endcounter);
	time_cost = ( (double)(endcounter.QuadPart - startcounter.QuadPart) /(double)frequency.QuadPart ) *1000;
	fprintf(timefp, "Time cost : %f ", time_cost);

	////	JacobiSolverOnCPU
	QueryPerformanceCounter(&startcounter);
	JacobiSolverOnCPU(&M[0], x, b, Dim);
	QueryPerformanceCounter(&endcounter);
	time_cost = ( (double)(endcounter.QuadPart - startcounter.QuadPart) /(double)frequency.QuadPart ) *1000;
	fprintf(timefp, "Time cost : %f ", time_cost);
*/

/*
	////	SucOverRelaxSolverOnCPU
	QueryPerformanceCounter(&startcounter);
	SucOverRelaxSolverOnCPU(&M[0], x, b, Dim);
	QueryPerformanceCounter(&endcounter);
	time_cost = ( (double)(endcounter.QuadPart - startcounter.QuadPart) /(double)frequency.QuadPart ) *1000;
	fprintf(timefp, "Time cost : %f ", time_cost);
	
	////	SteepestDescentSolverOnCPU
	QueryPerformanceCounter(&startcounter);
	SteepestDescentSolverOnCPU(&M[0], x, b, Dim);
	QueryPerformanceCounter(&endcounter);
	time_cost = ( (double)(endcounter.QuadPart - startcounter.QuadPart) /(double)frequency.QuadPart ) *1000;
	fprintf(timefp, "Time cost : %f ", time_cost);
*/

/*	////	ConjugateGradientSolverOnCPU
	ConjugateGradientSolverOnCPU(&M[0], x, b, Dim);
*/
	////////////////	On GPU	////////////////
	////	ConjugateGradientSolverOnGPU
	ConjugateGradientSolverOnGPU(&M[0], x, b, Dim);

	////	GaussSeidelSolverOnGPU
	//GaussSeidelSolverOnGPU(&M[0], x, b, Dim);

	////	JacobiSolverOnGPU
	//JacobiSolverOnGPU(&M[0], x, b, Dim);


	for(i=0; i<Dim; i++)
	{
		result[i] = 0.0;
		for(j=0; j<Dim; j++)
			result[i] += M[i*Dim+j]*x[j];
		result[i] = result[i] - b[i];
	}

	FILE *fp = fopen("..\\linearsolver.txt","w");
	fprintf(fp,"Input Matrix\n");
	for(i=0; i<Dim; i++)
	{
		for(j=0; j<Dim; j++)
			fprintf(fp,"%4.3f  ",M[i*Dim+j]);
		fprintf(fp,"\n");
	}
	fprintf(fp,"Result\n");
	for(i=0; i<Dim; i++)
	{
		fprintf(fp,"x=%4.3f b=%4.3f error=%4.3f\n",x[i], b[i], result[i]);
	}
	fclose(fp);

	return;
}