sparsematrixvectormultiply1.cpp

PDE simulator on GPU.

// SparseMatrixVectorMultiply1.cpp: implementation of the SparseMatrixVectorMultiply class.
//
//////////////////////////////////////////////////////////////////////

#include "stdafx.h"
#include "shadow.h"
#include "SparseMatrixOnGPU1.h"
#include "SparseMatrixVectorMultiply1.h"

#include "FragmentProgram.h"
#include "CommonHW.h"
#include "pbuffer.h"

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

//////////////////////////////////////////////////////////////////////
// Modified Bolz Method in Siggraph2003
//////////////////////////////////////////////////////////////////////

//////////////////////////////////////////////////////////////////////
/// Implement y = A*x
///           Vector = Matrix * Vector 
/// GL_NV_fragment_program
/// 
///           DO NOT WORK
///
//////////////////////////////////////////////////////////////////////
void BolzSparseMatVectAnotherWay(const CBolzMatrix& A, const CVectorOnGPU& x, CVectorOnGPU& result)
{
	int i,nPass = A.GetPassNum();
	////////////////////////////////////////////
	///If the size is compatible, then error occurs
	if(A.GetDim() != x.GetVectorLength())
	{
		AfxMessageBox("Error to Operate on Matrix-Vector Multiply for different size!");
		return;
	}

	////////////////////////////////////////////
	///Ax---->Matrix-Vector Product
	unsigned int ShaderID;

	unsigned int _iTextureTarget = x.GetTextureTarget();
	int Width  = x.GetWidth();
	int Height = x.GetHeight(); 

	if(_iTextureTarget == GL_TEXTURE_2D)
	{
		AfxMessageBox("Not support GL_TEXTURE_2D in current system(BolzSparseMatVectAnotherWay) for the reason of float precision!");
		return;
		//LoadFragmentProgramFromFile("..\\Bolzmatvector.fp", ShaderID);
	}
	else
	{
		ShaderID = GlobeShaderID[CL_BOLZ_SPARSE_MAT_VECT].ShaderID;
		//LoadFragmentProgramFromFile("..\\BolzmatvectorRECT.fp", ShaderID);
	}

	//////////////////////////////////////////////
	///// Clear the framebuffer firstly to avoid the previous computation
	///// This step is important to avoid computation errors
	glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);	
	glDisable(GL_DEPTH_TEST);

	///// Bind the Fragment Program to current drawing
	glEnable(GL_FRAGMENT_PROGRAM_NV);

	unsigned char strTexPara[] = "TexPARA";
	int len = strlen((const char*)strTexPara);
	float TexelPara[4];
	TexelPara[0] = 0;
	TexelPara[1] = A.GetWidth();
	TexelPara[2] = 0;
	TexelPara[3] = 0;
	
	glBindProgramNV(GL_FRAGMENT_PROGRAM_NV, ShaderID);

	//////////////////////////////////////////////////////////////////////////
	//# Vector X  ----TEX0 
	//# Diagonal Element   ----TEX1 -- AiiTexID  // record the diagonal elements
	//# NonZero NonDiagonal----TEX2 -- AijTexID; // record the nondiagonal & nonzero elements
	//# For reference in X ----TEX3 -- CTexID;   // record the according position for Vector X which multiply with Aij
	//# First NonZero Postion--TEX4 -- RTexID;   // record the first non zero & non diagonal element position J
	//# Previous  ====TEX5 //////record the previous result
	//////////////////////////////////////////////
	/// Tex0 --x
	glActiveTextureARB(GL_TEXTURE0_ARB);
	glBindTexture(_iTextureTarget, x.GetTextureID());	

	//////////////////////////////////////////////
	/// Tex1 ---A.Diagonal Aii Texture
	glActiveTextureARB(GL_TEXTURE1_ARB);
	glBindTexture(_iTextureTarget, A.GetAiiTextureID());	

	//////////////////////////////////////////////
	/// Tex2 ---A.Non Diagonal Aij Texture
	glActiveTextureARB(GL_TEXTURE2_ARB);
	glBindTexture(_iTextureTarget, A.GetAijTextureID());

	//////////////////////////////////////////////
	/// Tex3 ---A.C Texture
	glActiveTextureARB(GL_TEXTURE3_ARB);
	glBindTexture(_iTextureTarget, A.GetCTextureID());

	//////////////////////////////////////////////
	/// Tex4 ---A.R Texture
	glActiveTextureARB(GL_TEXTURE4_ARB);
	glBindTexture(_iTextureTarget, A.GetRTextureID());

	//////////////////////////////////////////////
	/// Tex5 ---result
	glActiveTextureARB(GL_TEXTURE5_ARB);
	glBindTexture(_iTextureTarget, result.GetTextureID());
	result.SetZero();
	glActiveTextureARB(GL_TEXTURE0_ARB);

	///Following: 
	//nPass = 1;
	for(i=0; i<nPass; i++)
	{	
		/// Transfer the Pass No. into fragment program
		TexelPara[0] = (float)i;
		glProgramNamedParameter4fvNV(ShaderID,  len, strTexPara, TexelPara);
		
		glBegin( GL_QUADS );
			glVertex2i( 0, 0 );
			glVertex2i( Width, 0 );
			glVertex2i( Width, Height );
			glVertex2i( 0, Height );
		glEnd();

		// Now we need to copy the resulting pixels into the output
		// save the result texture into result
		glActiveTextureARB(GL_TEXTURE5_ARB);
		glBindTexture(_iTextureTarget, result.GetTextureID());
		glCopyTexSubImage2D(_iTextureTarget, 0, 0, 0, 0, 0, Width, Height);
	}

	glDisable(GL_FRAGMENT_PROGRAM_NV);


	return;
}

//////////////////////////////////////////////////////////////////////////////////////////
void BolzSparseMatVect(const CBolzMatrix& A, const CVectorOnGPU& x, CVectorOnGPU& result)
{
	int i,nPass = A.GetPassNum();
	////////////////////////////////////////////
	///If the size is compatible, then error occurs
	if(A.GetDim() != x.GetVectorLength())
	{
		AfxMessageBox("Error to Operate on Matrix-Vector Multiply for different size!");
		return;
	}

	////////////////////////////////////////////
	///Ax---->Matrix-Vector Product
	unsigned int ShaderID, LastShaderID;

	unsigned int _iTextureTarget = x.GetTextureTarget();
	int Width  = x.GetWidth();
	int Height = x.GetHeight(); 

	if(_iTextureTarget == GL_TEXTURE_2D)
	{
		AfxMessageBox("Not support GL_TEXTURE_2D in current system(BolzSparseMatVect) for the reason of float precision!");
		return;
	}

	ShaderID = GlobeShaderID[CL_BOLZ_SPARSE_MAT_VECT].ShaderID;
	LastShaderID = GlobeShaderID[CL_BOLZ_SPARSE_MAT_VECT_LAST].ShaderID;
	//LoadFragmentProgramFromFile("..\\BolzmatvectorRECT.fp", ShaderID);
	//LoadFragmentProgramFromFile("..\\BolzmatvectorRECTLast.fp", LastShaderID);

	//////////////////////////////////////////////
	///// Clear the framebuffer firstly to avoid the previous computation
	///// This step is important to avoid computation errors
	//glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);	
	//glDisable(GL_DEPTH_TEST);

	unsigned char strTexPara[] = "TexPARA";
	int len = strlen((const char*)strTexPara);
	float TexelPara[4]={0,0,0,0};
	TexelPara[1] = A.GetIndexWidth();
	TexelPara[2] = 1.0/TexelPara[1];

	//////////////////////////////////////////////////////////////////////////
	//# Vector X  ----TEX0 
	//# Diagonal Element   ----TEX1 -- AiiTexID  // record the diagonal elements
	//# NonZero NonDiagonal----TEX2 -- AijTexID; // record the nondiagonal & nonzero elements
	//# For reference in X ----TEX3 -- CTexID;   // record the according position for Vector X which multiply with Aij
	//# First NonZero Postion--TEX4 -- RTexID;   // record the first non zero & non diagonal element position J
	//# Previous  ====TEX5 //////record the previous result
	////##################################################################
	for(i=0;i<2;i++)
	{
		pbuffer[i]->BeginCapture();
		glClear(GL_COLOR_BUFFER_BIT);	

		///// Bind the Fragment Program to current drawing
		glEnable(GL_FRAGMENT_PROGRAM_NV);
		glBindProgramNV(GL_FRAGMENT_PROGRAM_NV, ShaderID);
		//////////////////////////////////////////////
		/// Tex0 --x
		glActiveTextureARB(GL_TEXTURE0_ARB);
		glBindTexture(_iTextureTarget, x.GetTextureID());	

		//////////////////////////////////////////////
		/// Tex1 ---A.Diagonal Aii Texture
		glActiveTextureARB(GL_TEXTURE1_ARB);
		glBindTexture(_iTextureTarget, A.GetAiiTextureID());	

		//////////////////////////////////////////////
		/// Tex2 ---A.Non Diagonal Aij Texture
		glActiveTextureARB(GL_TEXTURE2_ARB);
		glBindTexture(_iTextureTarget, A.GetAijTextureID());

		//////////////////////////////////////////////
		/// Tex3 ---A.C Texture
		glActiveTextureARB(GL_TEXTURE3_ARB);
		glBindTexture(_iTextureTarget, A.GetCTextureID());

		//////////////////////////////////////////////
		/// Tex4 ---A.R Texture
		glActiveTextureARB(GL_TEXTURE4_ARB);
		glBindTexture(_iTextureTarget, A.GetRTextureID());

		//////////////////////////////////////////////
		/// Tex5 ---result
		glActiveTextureARB(GL_TEXTURE5_ARB);
		glBindTexture(_iTextureTarget, result.GetTextureID());
		result.SetZero();
		pbuffer[i]->EndCapture();
	}

	////##################################################################		
	int WhichPBuffer;
	for(i=0; i<nPass; i++)
	{
		WhichPBuffer = i%2;
		pbuffer[WhichPBuffer]->BeginCapture();
		/// Transfer the Pass No. into fragment program
		TexelPara[0] = (float)i;
		glProgramNamedParameter4fvNV(ShaderID,  len, strTexPara, TexelPara);

		glActiveTextureARB(GL_TEXTURE5_ARB);
		glBindTexture(_iTextureTarget, result.GetTextureID());

		glBegin( GL_QUADS );
			glVertex2i( 0, 0 );
			glVertex2i( Width, 0 );
			glVertex2i( Width, Height );
			glVertex2i( 0, Height );
		glEnd();

		// Now we need to copy the resulting pixels into the output
		// save the result texture into result
		glBindTexture(_iTextureTarget, result.GetTextureID());
		//glCopyTexSubImage2D(_iTextureTarget, 0, 0, 0, 0, 0, Width, Height);
		pbuffer[WhichPBuffer]->EndCapture();
	}

	////##################################################################		
	WhichPBuffer = i%2;
	pbuffer[WhichPBuffer]->BeginCapture();
	//// The last Pass to sum all components
	glBindProgramNV(GL_FRAGMENT_PROGRAM_NV, LastShaderID);
	glBegin( GL_QUADS );
		glVertex2i( 0, 0 );
		glVertex2i( Width, 0 );
		glVertex2i( Width, Height );
		glVertex2i( 0, Height );
	glEnd();

	// Now we need to copy the resulting pixels into the output
	// save the result texture into result
	glBindTexture(_iTextureTarget, result.GetTextureID());
	//glCopyTexSubImage2D(_iTextureTarget, 0, 0, 0, 0, 0, Width, Height);
	pbuffer[WhichPBuffer]->EndCapture();

/*
	for(i=1;i>=0;i--)
	{
		pbuffer[i]->BeginCapture();
		glDisable(GL_FRAGMENT_PROGRAM_NV);
		pbuffer[i]->EndCapture();
	}
*/

	return;
}

//////////////////////////////////////////////////////////////////////////////////////////
void MultiBolzSparseMatVect(const CBolzMatrix& A, const CVectorOnGPU& x, CVectorOnGPU& result)
{
	int i,nPass = A.GetPassNum();
	////////////////////////////////////////////
	///If the size is compatible, then error occurs
	if(A.GetDim() != x.GetVectorLength())
	{
		AfxMessageBox("Error to Operate on Matrix-Vector Multiply for different size!");
		return;
	}

	////////////////////////////////////////////
	///Ax---->Matrix-Vector Product
	unsigned int ShaderID, LastShaderID;

	unsigned int _iTextureTarget = x.GetTextureTarget();
	int Width  = x.GetWidth();
	int Height = x.GetHeight(); 

	if(_iTextureTarget == GL_TEXTURE_2D)
	{
		AfxMessageBox("Not support GL_TEXTURE_2D in current system(BolzSparseMatVect) for the reason of float precision!");
		return;
	}

	ShaderID = GlobeShaderID[CL_MULTI_BOLZ_SPARSE_MAT_VECT].ShaderID;
	LastShaderID = GlobeShaderID[CL_MULTI_BOLZ_SPARSE_MAT_VECT_LAST].ShaderID;
	//LoadFragmentProgramFromFile("..\\MultiBolzmatvectorRECT.fp", ShaderID);
	//LoadFragmentProgramFromFile("..\\MultiBolzmatvectorRECTLast.fp", LastShaderID);

	//////////////////////////////////////////////
	///// Clear the framebuffer firstly to avoid the previous computation
	///// This step is important to avoid computation errors
	//glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);	
	//glDisable(GL_DEPTH_TEST);

	unsigned char strTexPara[] = "TexPARA";
	int len = strlen((const char*)strTexPara);
	float TexelPara[4]={0,0,0,0};
	TexelPara[1] = A.GetIndexWidth();
	TexelPara[2] = 1.0/TexelPara[1];

	//////////////////////////////////////////////////////////////////////////
	//# Vector X  ----TEX0 
	//# Diagonal Element   ----TEX1 -- AiiRTexID  // record the diagonal elements
	//# First NonZero Postion--TEX1 -- RTexID;   // record the first non zero & non diagonal element position J
	//# NonZero NonDiagonal----TEX2 -- AijCTexID; // record the nondiagonal & nonzero elements
	//# For reference in X ----TEX2 -- CTexID;   // record the according position for Vector X which multiply with Aij
	//# Previous  ====TEX3 //////record the previous result
	////##################################################################
	for(i=0;i<2;i++)
	{
		pbuffer[i]->BeginCapture();
		glClear(GL_COLOR_BUFFER_BIT);	
		//////////////////////////////////////////////
		/// Tex0 --x
		glActiveTextureARB(GL_TEXTURE0_ARB);
		glBindTexture(_iTextureTarget, x.GetTextureID());	

		//////////////////////////////////////////////
		/// Tex1 ---A.Diagonal Aii Texture & R Texture
		glActiveTextureARB(GL_TEXTURE1_ARB);
		glBindTexture(_iTextureTarget, A.GetAiiTextureID());	

		//////////////////////////////////////////////
		/// Tex2 ---A.Non Diagonal Aij Texture & C Texture
		glActiveTextureARB(GL_TEXTURE2_ARB);
		glBindTexture(_iTextureTarget, A.GetAijTextureID());

		//////////////////////////////////////////////
		/// Tex3 ---result
		glActiveTextureARB(GL_TEXTURE3_ARB);
		glBindTexture(_iTextureTarget, result.GetTextureID());
		result.SetZero();
		glEnable(GL_FRAGMENT_PROGRAM_NV);
		glBindProgramNV(GL_FRAGMENT_PROGRAM_NV, ShaderID);

		pbuffer[i]->EndCapture();
	}

	////##################################################################
	int WhichPBuffer;
	for(i=0; i<nPass; i++)
	{	
		WhichPBuffer = i%2;
		pbuffer[WhichPBuffer]->BeginCapture();
		/// Transfer the Pass No. into fragment program
		TexelPara[0] = (float)i;
		glProgramNamedParameter4fvNV(ShaderID,  len, strTexPara, TexelPara);
		
		glBegin( GL_QUADS );
			glVertex2i( 0, 0 );
			glVertex2i( Width, 0 );
			glVertex2i( Width, Height );
			glVertex2i( 0, Height );
		glEnd();

		// Now we need to copy the resulting pixels into the output
		// save the result texture into result
		glBindTexture(_iTextureTarget, result.GetTextureID());
		//glCopyTexSubImage2D(_iTextureTarget, 0, 0, 0, 0, 0, Width, Height);
		pbuffer[WhichPBuffer]->EndCapture();
	}

	////##################################################################		
	WhichPBuffer = i%2;
	pbuffer[WhichPBuffer]->BeginCapture();
	//// The last Pass to sum all components
	glBindProgramNV(GL_FRAGMENT_PROGRAM_NV, LastShaderID);
	glBegin( GL_QUADS );
		glVertex2i( 0, 0 );
		glVertex2i( Width, 0 );
		glVertex2i( Width, Height );
		glVertex2i( 0, Height );
	glEnd();

	// Now we need to copy the resulting pixels into the output
	// save the result texture into result
	glBindTexture(_iTextureTarget, result.GetTextureID());
	//glCopyTexSubImage2D(_iTextureTarget, 0, 0, 0, 0, 0, Width, Height);
	pbuffer[WhichPBuffer]->EndCapture();	

/*
	for(i=1;i>=0;i--)
	{
		pbuffer[i]->BeginCapture();
		glDisable(GL_FRAGMENT_PROGRAM_NV);
		pbuffer[i]->EndCapture();
	}
*/	
	return;
}