sparsematrixongpu3.cpp

PDE simulator on GPU.

// SparseMatrixOnGPU3.cpp: implementation of the SparseMatrixOnGPU3 class.
//
//////////////////////////////////////////////////////////////////////

#include "stdafx.h"
#include "shadow.h"
#include "SparseMatrixOnGPU3.h"

#include <math.h>
#include "CommonHW.h"
#include "InitOpenGL.h"

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

///////////////////////////////////////////////////////////////////////////////////////
/// Class: CKrugerMatrix
/// Method from Jens Kruger in SIGGRPAPH2003
/// Linear Algebra Operators for GPU Implementations of Numerical Algorithms
///////////////////////////////////////////////////////////////////////////////////////
CKrugerMatrix::CKrugerMatrix()
{
	ColorArray   = NULL;
	IndexArray   = NULL;
	VertexArray  = NULL;
	TextureArray = NULL;

	TotalNonZeroNumOfEachColumn = NULL;

	_bInitialized = false;
}

//////////////////////////////////////////////////////////////////////
// 
//////////////////////////////////////////////////////////////////////
CKrugerMatrix::~CKrugerMatrix()
{
	/////Release all on Main Memory
	int i;
	if(ColorArray != NULL)
	{
		for(i=0; i<_iDim; i++)
		  delete[] ColorArray[i];
		delete[] ColorArray;
		ColorArray = NULL;
	}
	
	if(IndexArray != NULL)
	{
		for(i=0; i<_iDim; i++)
			delete[] IndexArray[i];
		delete[] IndexArray;
		IndexArray = NULL;
	}

	if(VertexArray != NULL)
	{
		for(i=0; i<_iDim; i++)
			delete[] VertexArray[i];
		delete[] VertexArray;
		VertexArray = NULL;
	}

	if(TextureArray != NULL)
	{
		for(i=0; i<_iDim; i++)
			delete[] TextureArray[i];
		delete[] TextureArray;
		TextureArray = NULL;
	}

}

////////////////////////////////////////////////////////////////////////////////////////////
/// Arrange the Sparse matrix M into vertex arrays
/// M should be square
/// M[DimSize][DimSize] 
////////////////////////////////////////////////////////////////////////////////////////////
void CKrugerMatrix::SetData(float *M, int DimSize)
{
	int i,j,k;

	//// M has N*N elements, here DimSize means total elements number
	_iDim = DimSize;

	//// Previous Method: Texture resolution = res*res (NOTE: square)
	//// Last Modified Date: 2003/10/23
	//// ------To get appropriate Texture Size to reduce reduction
	//// Texture size = Height*Width
	CalOptTexRSize(_iDim, _iHeight, _iWidth);
/*	
	if(IsPowerOfTwo(_iWidth) && IsPowerOfTwo(_iHeight))
	{
		_iTextureTarget = GL_TEXTURE_2D;
	}
	else
		_iTextureTarget = GL_TEXTURE_RECTANGLE_NV;
*/
	_iTextureTarget = g_TextureTarget;

	TotalNonZeroNumOfEachColumn = new int[_iDim];
	if(TotalNonZeroNumOfEachColumn == NULL)
	{
		AfxMessageBox("Error to Allocate Memory for TotalNonZeroNumOfEachColumn in CKrugerMatrix::SetData!");
		return;
	}

	/// Save the non-zero elements of each column into ONE vertex array
	/// First in Column
	//int MaxNonZeroNum = 0;
	for(j=0; j<_iDim; j++)
	{
		TotalNonZeroNumOfEachColumn[j] = 0;
		for(i=0; i<_iDim; i++)
		{
			////////sum the Non-zero elements 
			if(fabs(M[i*_iDim + j]) > 1.0e-3)
			{
				TotalNonZeroNumOfEachColumn[j]++;
			}
			//if(MaxNonZeroNum < TotalNonZeroNumOfEachColumn[j])
			//	MaxNonZeroNum = TotalNonZeroNumOfEachColumn[j];

		}
	}

	//// Allocate the memory for the arrays
	VertexArray  = new int*[_iDim];
	ColorArray   = new float*[_iDim];
	TextureArray = new float*[_iDim];
	IndexArray   = new unsigned int*[_iDim];
	for(i=0; i<_iDim; i++)
	{
		VertexArray[i]  = new int[TotalNonZeroNumOfEachColumn[i]*2];
		ColorArray[i]   = new float[TotalNonZeroNumOfEachColumn[i]*3];
		TextureArray[i] = new float[TotalNonZeroNumOfEachColumn[i]*2];
		IndexArray[i]   = new unsigned int[TotalNonZeroNumOfEachColumn[i]];

	}
	if(VertexArray == NULL || ColorArray == NULL || TextureArray == NULL || IndexArray == NULL)
	{
		AfxMessageBox("Error to allocate memory in CKrugerMatrix::SetData!");
		return;
	}

	////Generate the Vertex Array according to the non-zero elements
	/// First in Column
	for(j=0; j<_iDim; j++)
	{
		k = 0;
		for(i=0; i<_iDim; i++)
		{
			////////Record the Non-zero elements 
			if(fabs(M[i*_iDim + j]) > 1.0e-3)
			{
				//// Vertex Position
				VertexArray[j][k*2]   = i%_iWidth ;  //x
				VertexArray[j][k*2+1] = i/_iWidth ;  //y
				
				if(_iTextureTarget == GL_TEXTURE_2D)
				{
					TextureArray[j][k*2]   = (float)(j%_iWidth)/(_iWidth-1);
					TextureArray[j][k*2+1] = (float)(j/_iWidth)/(_iWidth-1);
				}
				else ///GL_TEXTURE_RECTANGLE
				{
					TextureArray[j][k*2]   = (float)(j%_iWidth);
					TextureArray[j][k*2+1] = (float)(j/_iWidth);
				}

				ColorArray[j][k*3]   = M[i*_iDim + j];
				ColorArray[j][k*3+1] = 0.0f;
				ColorArray[j][k*3+2] = 0.0f;

				IndexArray[j][k] = k;
				k++;
			}
		}
	}

	_bInitialized = true;
/*
	FILE *fp = fopen("testttt.txt", "w");
	for(i=0;i<_iDim;i++)
	{
		fprintf(fp,"Column%d NonZeroNum=%d\n",i, TotalNonZeroNumOfEachColumn[i]);
		for(j=0;j<TotalNonZeroNumOfEachColumn[i];j++)
		{
			fprintf(fp,"Index%d-Pos=%d,%d, Color=%f, Tex Coord=%f,%f\n",
				IndexArray[i][j],
				VertexArray[i][j*2],VertexArray[i][j*2+1],
				ColorArray[i][j*3],
				TextureArray[i][j*2], TextureArray[i][j*2+1]);
		}

	}
	fclose(fp);
*/
	return;
}
//////////////////////////////////////////////////////////////////////////////