sparsematrixongpu1.cpp

PDE simulator on GPU.

// SparseMatrixOnGPU1.cpp: implementation of the SparseMatrixOnGPU class.
//
//////////////////////////////////////////////////////////////////////

#include "stdafx.h"
#include "shadow.h"
#include "SparseMatrixOnGPU1.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

//////////////////////////////////////////////////////////////////////
// FOR sparse matrix operation --- from Bolze SIGGRAPH2003
// Class: CBolzMatrix
//////////////////////////////////////////////////////////////////////
CBolzMatrix::CBolzMatrix()
{
	AiiTexID = 0;       ///////Aii
	AijTexID = 0;

	RTexID = 0;
	CTexID = 0;

	_bInitialized = false;
	_bFloat = g_Float;
	_iChannelsNum = 1;
	
}

//////////////////////////////////////////////////////////////////////
// 
//////////////////////////////////////////////////////////////////////
CBolzMatrix::~CBolzMatrix()
{
	glDeleteTextures(1,&AiiTexID);
	glDeleteTextures(1,&AijTexID);
	glDeleteTextures(1,&RTexID);
	glDeleteTextures(1,&CTexID);
}

//////////////////////////////////////////////////////////////////////
/// Arrange the input data into textures for computations on GPU,
/// at the same time, create the connection information textures
/// A----DimSize*DimSize
/// Currently, just One Channel with TEX_RECTANGLE_NV supported
///
///                     THIS one DO NOT work
///
//////////////////////////////////////////////////////////////////////
void CBolzMatrix::SetDataAnotherWay(float *A, int DimSize)
{
	int i,j,k,m,n;
	//////////////////////////////////////////////////////////////////////
	//// The length of vector x and y
	_iDim = DimSize;

	GLuint iInternalFormat;
	GLuint iFormat;

	float *fInputMatrix = NULL;
	if( _iChannelsNum == 1)
	{
		//// Determine the format of texture according to specifications
		if(_bFloat == true )
		{
			iInternalFormat = GL_FLOAT_R32_NV;
			iFormat = GL_RED;
		}
		else
		{
			iInternalFormat = GL_LUMINANCE;
			iFormat = GL_RED;
		}		

	}
	else////4 channels
	{
		if(_bFloat == true )
		{
			iInternalFormat = GL_FLOAT_RGBA32_NV;
			iFormat = GL_RGBA;
		}
		else
		{
			iInternalFormat = GL_RGBA8;
			iFormat = GL_RGBA;
		}
	}
/*
	if(IsPowerOfTwo(_iWidth) && IsPowerOfTwo(_iHeight) && _bFloat == false)
	{
		_iTextureTarget = GL_TEXTURE_2D;
	}
	else
		_iTextureTarget = GL_TEXTURE_RECTANGLE_NV;
*/
	_iTextureTarget = g_TextureTarget;

	//////////////////////////////////////////////////////////////////////
	//// Texture--->AijTexID---->Matrix A non-diagonal elements aij
	////                    ---->ONLY Non-zero elements are stored
	//// Len = DimSize*DimSize
	//// Texture = XTexSize*XTexSize
	//// This step to Define the sizeof texture AijTexID
	_iMaxPassNum = 0; 
	int sum;
	int TotalNonZeroNum = 0;
	for(i=0;i<DimSize;i++)
	{
		sum = 0;
		for(j=0;j<DimSize;j++)
		{
			/////////Non-zero elements and Non-diagonal
			if(fabs(A[i*_iDim + j] > 1.0e-3) && i != j)
			{
				//////Compact the nonzero elements into AijTexID texture
				TotalNonZeroNum++;
				sum++;
			}
		}
		if(_iMaxPassNum < sum)
			_iMaxPassNum = sum;
		
	}

	///// Just ONE channel Here
	CalOptTexRSize(TotalNonZeroNum, _iIndexHeight, _iIndexWidth);
	CalOptTexRSize(_iDim, _iHeight, _iWidth);

	//////////////////////////////////////////////////////////////////////
	/////// Texture of Aii in 1 channel ----------_iHeight*_iWidth
	/////// Texture of Aij in 1 channel ----------_iIndexHeight*_iIndexWidth
	/////// Texture of C   in 3 channels----------_iIndexHeight*_iIndexWidth*3
	/////// Texture of R   in 3 channels----------_iHeight*_iWidth*3
	/////// Currently just in 1 Channel
	float *AiiData = new float[_iHeight*_iWidth];
	memset(AiiData,   0,sizeof(float)*_iHeight*_iWidth);

	float *ImageRData = new float[_iHeight*_iWidth*3];
	memset(ImageRData,0,sizeof(float)*_iHeight*_iWidth*3);
	
	float *AijData = new float[_iIndexHeight*_iIndexWidth];
	memset(AijData,   0,sizeof(float)*_iIndexHeight*_iIndexWidth);

	float *ImageCData = new float[_iIndexHeight*_iIndexWidth*3];
	memset(ImageCData,0,sizeof(float)*_iIndexHeight*_iIndexWidth*3);

	bool isFirstNonZero;
	k = m = n = 0;	
	for(i=0;i<DimSize;i++)
	{
		/////////Non-zero elements of each ROW
		isFirstNonZero = true;
		for(j=0;j<DimSize;j++)
		{
			AiiData[i] = A[i*_iDim + i];
			/////////Non-zero elements and Non-diagonal
			if(fabs(A[i*_iDim + j] > 1.0e-3) && i != j)
			{
				////// Just record No.1 non-zero of each ROW 
				////// Save its Index Position in ImageRData
				if(isFirstNonZero == true)
				{
					//////---RGB component save the texture coordinates of 
					//////---No.1 Non-zero element of each Row of Matrix A
					//////---Just G & B are used here for dependent texture
					////// ImageRData ---> AijTexID
					ImageRData[k++] = (float)(m%_iIndexWidth);    ///R--->U-->Aij
					ImageRData[k++] = (float)(m/_iIndexWidth);    ///G--->V-->
					ImageRData[k++] = 0;			              ///B
					isFirstNonZero = false;
				}
				
				////// Compact the nonzero elements into AijTexID texture
				AijData[m++] = A[i*_iDim + j];
				
				////// ImageRData ---> ImageCData ---> xj in accordance with Aij
				////// ImageCData ---> xj
				ImageCData[n++] = (float)(j%_iWidth);    ///R--->U--->Xj
				ImageCData[n++] = (float)(j/_iWidth);	 ///G			
				ImageCData[n++] = 0;                     ///B
			}
		}		
	}

	//////////////////////////////////////////////////////////////////////
	//////Generate the corresponding textures
	//////Texture for Aii (diagonal elements)
	glGenTextures(1,&AiiTexID);
	glBindTexture(_iTextureTarget,AiiTexID);
	glPixelStorei(GL_UNPACK_ALIGNMENT, 1);
	glTexEnvi(GL_TEXTURE_ENV, GL_TEXTURE_ENV_MODE, GL_REPLACE);
	glTexParameterf(_iTextureTarget, GL_TEXTURE_MAG_FILTER, GL_NEAREST);
	glTexParameterf(_iTextureTarget, GL_TEXTURE_MIN_FILTER, GL_NEAREST);	

    glTexParameterf(_iTextureTarget, GL_TEXTURE_WRAP_S, GL_CLAMP_TO_EDGE);
    glTexParameterf(_iTextureTarget, GL_TEXTURE_WRAP_T, GL_CLAMP_TO_EDGE);
	glTexImage2D( _iTextureTarget, 0, iInternalFormat, _iWidth, _iHeight, 0, iFormat, GL_FLOAT, AiiData );

	//////Texture for Aij (Non-diagonal & Non-zero elements)
	glGenTextures(1,&AijTexID);
	glBindTexture(_iTextureTarget,AijTexID);
	glPixelStorei(GL_UNPACK_ALIGNMENT, 1);
	glTexEnvi(GL_TEXTURE_ENV, GL_TEXTURE_ENV_MODE, GL_REPLACE);
	glTexParameterf(_iTextureTarget, GL_TEXTURE_MAG_FILTER, GL_NEAREST);
	glTexParameterf(_iTextureTarget, GL_TEXTURE_MIN_FILTER, GL_NEAREST);	

    glTexParameterf(_iTextureTarget, GL_TEXTURE_WRAP_S, GL_CLAMP_TO_EDGE);
    glTexParameterf(_iTextureTarget, GL_TEXTURE_WRAP_T, GL_CLAMP_TO_EDGE);
	glTexImage2D( _iTextureTarget, 0, iInternalFormat, _iIndexWidth, _iIndexHeight, 0, iFormat, GL_FLOAT, AijData );

	//////////////////////////////////////////////////////////////////////
	////Texture--->CTexID---->Indirectional Texture
	///            Same layout to Texture of Aij in RGB       
	glGenTextures(1,&CTexID);
	glBindTexture(_iTextureTarget,CTexID);
	glPixelStorei(GL_UNPACK_ALIGNMENT, 1);
	glTexEnvi(GL_TEXTURE_ENV, GL_TEXTURE_ENV_MODE, GL_REPLACE);
	glTexParameterf(_iTextureTarget, GL_TEXTURE_MAG_FILTER, GL_NEAREST);
	glTexParameterf(_iTextureTarget, GL_TEXTURE_MIN_FILTER, GL_NEAREST);	

    glTexParameterf(_iTextureTarget, GL_TEXTURE_WRAP_S, GL_CLAMP_TO_EDGE);
    glTexParameterf(_iTextureTarget, GL_TEXTURE_WRAP_T, GL_CLAMP_TO_EDGE);
	glTexImage2D( _iTextureTarget, 0, GL_FLOAT_RGBA32_NV, _iIndexWidth, _iIndexHeight, 0, GL_RGB, GL_FLOAT, ImageCData );

	//////////////////////////////////////////////////////////////////////
	////Texture--->RTexID---->Indirectional Texture 
	/// In Same Size with X texture (SIZE = XTexSize*XTexSize)
	glGenTextures(1,&RTexID);
	glBindTexture(_iTextureTarget,RTexID);
	glPixelStorei(GL_UNPACK_ALIGNMENT, 1);
	glTexEnvi(GL_TEXTURE_ENV, GL_TEXTURE_ENV_MODE, GL_REPLACE);
	glTexParameterf(_iTextureTarget, GL_TEXTURE_MAG_FILTER, GL_NEAREST);
	glTexParameterf(_iTextureTarget, GL_TEXTURE_MIN_FILTER, GL_NEAREST);	

    glTexParameterf(_iTextureTarget, GL_TEXTURE_WRAP_S, GL_CLAMP_TO_EDGE);
    glTexParameterf(_iTextureTarget, GL_TEXTURE_WRAP_T, GL_CLAMP_TO_EDGE);
	glTexImage2D( _iTextureTarget, 0, GL_FLOAT_RGBA32_NV, _iWidth, _iHeight, 0, GL_RGB, GL_FLOAT, ImageRData );

	
	delete[] AiiData;
	AiiData = NULL;
	delete[] AijData;
	AijData = NULL;
	delete[] ImageCData;
	ImageCData = NULL;
	delete[] ImageRData;
	ImageRData = NULL;

	_bInitialized = true;
	
	return;
}


void CBolzMatrix::SetData(float *A, int DimSize)
{
	int i,j,k,m,n;
	//////////////////////////////////////////////////////////////////////
	//// The length of vector x and y
	_iDim = DimSize;

	GLuint iInternalFormat;
	GLuint iFormat;

	float *fInputMatrix = NULL;
	if( _iChannelsNum == 1)
	{
		//// Determine the format of texture according to specifications
		if(_bFloat == true )
		{
			iInternalFormat = GL_FLOAT_R32_NV;
			iFormat = GL_RED;
		}
		else
		{
			iInternalFormat = GL_LUMINANCE;
			iFormat = GL_RED;
		}		

	}
	else////4 channels
	{
		if(_bFloat == true )
		{
			iInternalFormat = GL_FLOAT_RGBA32_NV;
			iFormat = GL_RGBA;
		}
		else
		{
			iInternalFormat = GL_RGBA8;
			iFormat = GL_RGBA;
		}
	}
/*
	if(IsPowerOfTwo(_iWidth) && IsPowerOfTwo(_iHeight) && _bFloat == false)
	{
		_iTextureTarget = GL_TEXTURE_2D;
	}
	else
		_iTextureTarget = GL_TEXTURE_RECTANGLE_NV;
*/
	_iTextureTarget = g_TextureTarget;

	//////////////////////////////////////////////////////////////////////
	//// Texture--->AijTexID---->Matrix A non-diagonal elements aij
	////                    ---->ONLY Non-zero elements are stored
	//// Len = DimSize*DimSize
	//// Texture = XTexSize*XTexSize
	//// This step to Define the sizeof texture AijTexID
	_iMaxPassNum = 0; 
	int sum;
	int TotalNonZeroNum = 0;
	int *len = new int[DimSize];
	for(i=0;i<DimSize;i++)
	{
		sum = 0;
		for(j=0;j<DimSize;j++)
		{
			/////////Non-zero elements and Non-diagonal
			if(fabs(A[i*_iDim + j] > 1.0e-3) && i != j)
			{