又是一個月了，這個月前面大多數的時間都在看網格變形的paper，後面這段時間又要搞Ogre，日子過得很混亂，完全沒有什麼時間來整理學過的東西。在這個月末的時候還是貼篇文章表示還在好好學習。

   網格變形的論文看得理解了一點後就開始嘗試實現，就開始到網上找稀疏矩陣庫。一開始我使用的矩陣庫是Eigen，風格跟matlab很像，但是沒有實現稀疏矩陣運算的功能，雖然它預留了跟其他幾個稀疏矩陣庫的介面，但是我試過的umfpack和superlu都在if(!lu_of_A.succeeded())這一步計算失敗。後來還找嘗試了taucs庫，但是同樣遇到了一個悲劇的問題，編譯成功的庫可以在命令列下使用，但是放到vs2008裡面卻總是無法連結成功。花了幾天的時間，最後還是決定自己封裝下umfpack和Eigen的介面。

  首先要了解稀疏矩陣的儲存格式，稀疏矩陣的標準格式是以列為主序的，使用三個陣列來存放它的資料，int *pAi,int * pAp,和 double *pAx，其中pAi存放存放每個非0元素的行號，pAp則存放每列開頭非0元素在pAx裡面的索引，pAx以列為主序存放非0元素，詳細的資訊可以從http://people.sc.fsu.edu/~jburkardt/pdf/hbsmc.pdf這裡找到，第7頁

  直接貼原始檔

  //eigen_umfpack_bridge.h

#ifndef EIGEN_UMFPACK_BRIDGE_K
#define EIGEN_UMFPACK_BRIDGE_K



#define EIGEN_YES_I_KNOW_SPARSE_MODULE_IS_NOT_STABLE_YET
#include <Eigen/Sparse>
using namespace Eigen;

void ConvEigenToUmfpackType(const SparseMatrix<double>&A,
							double* &pAx,int* &pAi,int* &pAp);

void DecomposeSparseMatrix(const int dim,const double *pAx,const int *pAi,
						   const int *pAp,void* &pNumeric);

void SolveLinEq(const double *pAx,const int *pAi,const int *pAp,
				const double *pB,void* &pNumeric,double *pX);


void FreeNumeric(void *pNumeric );


void WriteEigenMatrixToFile(const SparseMatrix<double>&A,char *szFile);

#endif
 

//eigen_umfpack_bridge.cpp

#include "eigen_umfpack_bridge.h"
#include "umfpack.h"

#include <fstream>
#include <string>
using namespace std;
void ConvEigenToUmfpackType(const SparseMatrix<double>&A,double* &pAx,int* &pAi,int* &pAp)
{
	int size = A.nonZeros();
	int cols = A.cols();
	pAx = new double[size];
	pAi = new int[size];
	pAp = new int[cols+1];
	int len = sizeof(int)*(cols+1);
	memset(pAp,-1,len);
	pAp[cols] = size;
	int ind = 0;

	for (int k=0; k<cols; ++k)
		for (SparseMatrix<double>::InnerIterator it(A,k); it; ++it)
		{
			double val = it.value();
			pAx[ind] = val;

			int r = it.row();   // row index			
			pAi[ind] = r;

			if(pAp[k] == -1)
			{
				pAp[k] = ind;
			}

			ind ++;			
		}

}

void DecomposeSparseMatrix(const int dim,const double *pAx,const int *pAi,
						   const int *pAp,void* &pNumeric)
{
	double *null = (double *) NULL ;

	void *Symbolic ;
	
	int n = dim;
	(void) umfpack_di_symbolic (n, n, pAp, pAi, pAx, &Symbolic, null, null) ;
	(void) umfpack_di_numeric (pAp, pAi, pAx, Symbolic, &pNumeric, null, null) ;
	umfpack_di_free_symbolic (&Symbolic) ;
}

void SolveLinEq(const double *pAx,const int *pAi,const int *pAp,
				const double *pB,void* &pNumeric,double *pX)
{
	double *null = (double *) NULL ;
	(void) umfpack_di_solve (UMFPACK_A, pAp, pAi, pAx, pX, pB, pNumeric, null, null) ;
}
void FreeNumeric(void *pNumeric )
{
	if(pNumeric)umfpack_di_free_numeric (&pNumeric) ;
}

void WriteEigenMatrixToFile(const SparseMatrix<double>&A,char *szFile)
{
	ofstream f(szFile);
	int row = A.rows();
	int col = A.cols();
	f<<row<<" "<<col<<endl;
	string strI="",strJ="",strX="";
	for (int k=0; k<col; ++k)
	{
		for (SparseMatrix<double>::InnerIterator it(A,k); it; ++it)
		{
			double val = it.value();
			
			int r = it.row();   // row index			
			char szI[10],szJ[10],szX[20];
			sprintf(szI,"%d \0",r);
			sprintf(szJ,"%d \0",k);
			sprintf(szX,"%f \0",val);
			strI += szI;
			strJ += szJ;
			strX += szX;
			
		}
	}
	f<<strI<<endl;
	f<<strJ<<endl;
	f<<strX<<endl;
	f.close();
}