sparsematrixongpu1.cpp

PDE simulator on GPU.

				//////Compact the nonzero elements into AijTexID texture
				TotalNonZeroNum++;
				sum++;
			}
		}
		if(_iMaxPassNum < sum)
			_iMaxPassNum = sum;
		len[i] = 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 R & G are used here for dependent texture
					////// ImageRData ---> AijTexID
					ImageRData[k++] = m;         ///R--->U,V-->Aij
					ImageRData[k++] = len[i];    ///G--->record the number nonzero elements
					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 Aij (Non-diagonal & Non-zero 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_RGB32_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_RGB32_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;

	delete[] len;
	len = NULL;

	_bInitialized = true;
	
	return;
}

//////////////////////////////////////////////////////////////////////////
/// Arrange the same layout data into ONE texture to reduce instructions
/// to fetch textures, so performance improved
//////////////////////////////////////////////////////////////////////////
void CBolzMatrix::SetMultiData(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)
			{
				//////Compact the nonzero elements into AijTexID texture
				TotalNonZeroNum++;
				sum++;
			}
		}
		if(_iMaxPassNum < sum)
			_iMaxPassNum = sum;
		len[i] = sum;
		
	}

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

	//////////////////////////////////////////////////////////////////////////
	/////// NOTE: Load the multi version fragment program
	/////// Rearrange Aii and R into ONE texture 
	/////// AiiRTex.r = Aii
	/////// AiiRTex.g = R.x
	/////// AiiRTex.b = R.y
	/////// Texture of Aii in 1 channel ----------_iHeight*_iWidth
	/////// Texture of R   in 2 channels----------_iHeight*_iWidth*3
	///////
	/////// Rearrange Aij and C into ONE texture
	/////// AijCTex.r = Aij
	/////// AijCTex.g = C.x
	/////// AijCTex.b = C.y
	/////// Texture of Aij in 1 channel ----------_iIndexHeight*_iIndexWidth
	/////// Texture of C   in 2 channels----------_iIndexHeight*_iIndexWidth*3
	//////////////////////////////////////////////////////////////////////////
	
	float *AiiRData = new float[_iHeight*_iWidth*3];
	memset(AiiRData,   0,sizeof(float)*_iHeight*_iWidth*3);

	float *AijCData = new float[_iIndexHeight*_iIndexWidth*3];
	memset(AijCData,   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;
		AiiRData[i*3] = A[i*_iDim + i];
		for(j=0;j<DimSize;j++)
		{
			/////////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 R & G are used here for dependent texture
					////// ImageRData ---> AijTexID
					AiiRData[i*3+1] = m;         ///R--->U,V-->Aij
					AiiRData[i*3+2] = len[i];    ///G--->record the number nonzero elements
					isFirstNonZero = false;
				}
				
				////// Compact the nonzero elements into AijTexID texture
				AijCData[m*3] = A[i*_iDim + j];
				
				////// ImageRData ---> ImageCData ---> xj in accordance with Aij
				////// ImageCData ---> xj
				AijCData[m*3+1] = (float)(j%_iWidth);    ///R--->U--->Xj
				AijCData[m*3+2] = (float)(j/_iWidth);	 ///G			
				m++;
			}
		}		
	}

	//////////////////////////////////////////////////////////////////////
	//////Generate the corresponding textures
	//////Texture for Aij (Non-diagonal & Non-zero 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, GL_FLOAT_RGB32_NV, _iWidth, _iHeight, 0, GL_RGB, GL_FLOAT, AiiRData );

	//////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, GL_FLOAT_RGB32_NV, _iIndexWidth, _iIndexHeight, 0, GL_RGB, GL_FLOAT, AijCData );

	delete[] AiiRData;
	AiiRData = NULL;
	delete[] AijCData;
	AijCData = NULL;
	delete[] len;
	len = NULL;
	_bInitialized = true;
	
	return;
}