稀疏矩阵库umfpack和Eigen结合

   又是一个月了，这个月前面大多数的时间都在看网格变形的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 indexpAi[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 indexchar 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();}