sparsematrixongpu2.cpp

PDE simulator on GPU.

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

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

//////////////////////////////////////////////////////////////////////
// From Brook System  2004/1/1
// Class: CBrookMatrix
//////////////////////////////////////////////////////////////////////
CBrookMatrix::CBrookMatrix()
{
	_iDataTexID   = 0;         // record the non zero elements
	_iPosTexID    = 0;         // record the accordingly position(Column J)of each non zero element
	_iStartTexID  = 0;         // record the first non zero element's position(Column J) of each row
	_iLenTexID    = 0;         // record the length of each row which sum non zero
	
	_bInitialized = false;
	_bFloat = g_Float;
	_iChannelsNum = 1;
}

//////////////////////////////////////////////////////////////////////
// 
//////////////////////////////////////////////////////////////////////
CBrookMatrix::~CBrookMatrix()
{
	glDeleteTextures(1,&_iDataTexID);
	glDeleteTextures(1,&_iPosTexID);
	glDeleteTextures(1,&_iStartTexID);
	glDeleteTextures(1,&_iLenTexID);
}

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

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

	float *fInputVector = 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
	{
		//// Determine the format of texture according to specifications
		if(_bFloat == true )
		{
			iInternalFormat = GL_FLOAT_RGBA32_NV;
			iFormat = GL_RGBA;
		}
		else
		{
			iInternalFormat = GL_RGBA8;
			iFormat = GL_RGBA;
		}
	}
/*	
	if(IsPowerOfTwo(_iWidth) && IsPowerOfTwo(_iHeight))
	{
		_iTextureTarget = GL_TEXTURE_2D;
	}
	else
		_iTextureTarget = GL_TEXTURE_RECTANGLE_NV;
*/
	_iTextureTarget = g_TextureTarget;

	int TotalNonZeroNum = 0;
	for(i=0; i<_iDim; i++)
	{
		for(j=0; j<_iDim; j++)
		{
			////////sum the Non-zero elements 
			if(fabs(M[i*_iDim + j]) > 1.0e-3)
			{
				TotalNonZeroNum++;
			}
		}
	}

	///// Just ONE channel Here
	CalOptTexRSize(TotalNonZeroNum, _iDataHeight, _iDataWidth);
	CalOptTexRSize(_iDim, _iIndexHeight, _iIndexWidth);

	float *elem = new float[_iDataHeight*_iDataWidth];
	float *ipos = new float[_iDataHeight*_iDataWidth*3];
	float *start= new float[_iIndexHeight*_iIndexWidth];
	float *len  = new float[_iIndexHeight*_iIndexWidth];
	memset(elem,  0 , sizeof(float)*_iDataHeight*_iDataWidth);
	memset(ipos,  0 , sizeof(float)*_iDataHeight*_iDataWidth*3);
	memset(start, 0 , sizeof(float)*_iIndexHeight*_iIndexWidth);
	memset(len,   0 , sizeof(float)*_iIndexHeight*_iIndexWidth);

	int sum, first;
	l = k = 0;
	_iMaxPassNum = 0; 
	for(i=0; i<_iDim; i++)
	{
		sum   = 0;
		first = _iDataHeight*_iDataWidth;
		for(j=0; j<_iDim; j++)
		{
			////////sum the Non-zero elements 
			if(fabs(M[i*_iDim + j]) > 1.0e-3)
			{
				elem[k] = M[i*_iDim + j];
				ipos[l++] = (float)(j%_iIndexWidth); /// x
				ipos[l++] = (float)(j/_iIndexWidth); /// y
				ipos[l++] = 0.0f;
				if(first > k)
					first = k;							
				k++;
				sum++;
	
			}			
		}
		if(_iMaxPassNum < sum)
			_iMaxPassNum = sum;
		len[i]   = sum;
		start[i] = first;
	}


	////Generate Element Textures
    glGenTextures(1,&_iDataTexID);
	glBindTexture(_iTextureTarget, _iDataTexID);
	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, _iDataWidth, _iDataHeight, 0, iFormat, GL_FLOAT, elem );

	delete[] elem;
	elem = NULL;

	
	////Generate Position of j Textures
	glGenTextures(1, &_iPosTexID);
	glBindTexture(_iTextureTarget, _iPosTexID);
	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, _iDataWidth, _iDataHeight, 0, iFormat, GL_FLOAT, ipos );
	glTexImage2D(_iTextureTarget, 0, GL_FLOAT_RGB32_NV, _iDataWidth, _iDataHeight, 0, GL_RGB, GL_FLOAT, ipos );

	delete[] ipos;
	ipos = NULL;

	////Generate Start Textures to record the first Column ID of NonZero in each row
	glGenTextures(1,&_iStartTexID);
	glBindTexture(_iTextureTarget, _iStartTexID);
	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, start );
	
	delete[] start;
	start = NULL;


	////Generate Length Textures
	glGenTextures(1,&_iLenTexID);
	glBindTexture(_iTextureTarget, _iLenTexID);
	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, len );

	delete[] len;
	len = NULL;	

	_bInitialized = true;

	return;
}

//////////////////////////////////////////////////////////////////////////
/// Save the same layout data into ONE texture by filling with RGBA
/// Just support GL_TEXTURE_RECTANGLE_NV
/// ELEMENT data are stored in ONE channel
///	elem= new float[TotalNonZeroNum];	 --|
/// ipos = new float[TotalNonZeroNum*3]; --|---same layout, so in ONE texture
///	start= new float[_iDim]; ---|
///	len  = new float[_iDim]; ---|---same layout, so in ONE texture
//////////////////////////////////////////////////////////////////////////
void CBrookMatrix::SetMultiData(float *M, int DimSize)
{
	int i,j,k,l;

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

	float *fInputVector = NULL;
	if( _iChannelsNum != 1 || _bFloat == false)
	{
		AfxMessageBox("Not support in this method(CBrookMatrix::SetMultiData)!");
		return;
	}

	_iTextureTarget = GL_TEXTURE_RECTANGLE_NV;

	int TotalNonZeroNum = 0;
	for(i=0; i<_iDim; i++)
	{
		for(j=0; j<_iDim; j++)
		{
			////////sum the Non-zero elements 
			if(fabs(M[i*_iDim + j]) > 1.0e-3)
			{
				TotalNonZeroNum++;
			}
		}
	}

	///// Just ONE channel Here
	CalOptTexRSize(TotalNonZeroNum, _iDataHeight, _iDataWidth);
	CalOptTexRSize(_iDim, _iIndexHeight, _iIndexWidth);

	float *elem_ipos = new float[_iDataHeight*_iDataWidth*3];
	memset(elem_ipos, 0 , sizeof(float)*_iDataHeight*_iDataWidth*3);
	float *start_len= new float[_iIndexHeight*_iIndexWidth*3];
	memset(start_len, 0 , sizeof(float)*_iIndexHeight*_iIndexWidth*3);
	int sum, first;
	l = k = 0;
	_iMaxPassNum = 0; 
	for(i=0; i<_iDim; i++)
	{
		sum   = 0;
		first = _iDataHeight*_iDataWidth;
		for(j=0; j<_iDim; j++)
		{
			////////sum the Non-zero elements 
			if(fabs(M[i*_iDim + j]) > 1.0e-3)
			{
				elem_ipos[k*3] = M[i*_iDim + j];
				elem_ipos[k*3+1] = (float)(j%_iIndexWidth); /// x
				elem_ipos[k*3+2] = (float)(j/_iIndexWidth); /// y
				if(first > k)
					first = k;							
				k++;
				sum++;
	
			}			
		}
		if(_iMaxPassNum < sum)
			_iMaxPassNum = sum;
		start_len[i*3]   = sum;
		start_len[i*3+1] = first;
		start_len[i*3+2] = 0;
	}


	////Generate Element Textures
    glGenTextures(1,&_iDataTexID);
	glBindTexture(_iTextureTarget, _iDataTexID);
	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_RGB32_NV, _iDataWidth, _iDataHeight, 0, GL_RGB, GL_FLOAT, elem_ipos );

	delete[] elem_ipos;
	elem_ipos = NULL;
	
	////Generate Start Textures to record the first Column ID of NonZero in each row
	glGenTextures(1,&_iStartTexID);
	glBindTexture(_iTextureTarget, _iStartTexID);
	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_RGB32_NV, _iIndexWidth, _iIndexHeight, 0, GL_RGB, GL_FLOAT, start_len );
	
	delete[] start_len;
	start_len = NULL;

	_bInitialized = true;

	return;
}