最近仿08年TOG上一篇骨架提取的文章Skeleton extraction by mesh contraction，其中涉及到線性方程組的最小二乘解問題，即Ax = b。
最開始使用了Armadillo庫進行求解，程序寫完后發現矩陣A的規模與頂點數的平方成正比，不使用稀疏矩陣的話只能計算很小的模型，但Armadillo沒有提供稀疏矩陣模塊。聽說Eigen庫有稀疏矩陣模塊，又查了下Eigen庫，但是發現Eigen庫的稀疏矩陣求解線性方程組的功能只能用於A為方陣的情況。最后考慮用matlab結合vs2010的方式實現。
Armadillo和Eigen都是很綠色的線性代數庫，都是以泛型編程的方式實現。Armadillo非常簡練，文檔也小巧精悍，上手很快，底層依賴lapack和blas庫，環境配置方式寫在了之前的一篇文檔里http://www.cnblogs.com/youthlion/archive/2012/05/15/2501465.html。Eigen庫更加重量級一些，功能更加全面，文檔詳細，不依賴於任何其他底層庫。環境配置也很簡單，不多說了。下面主要記錄一下vs2010和matlab混用的方法。
其實matlab提供了多種工具，既可以在matlab中調用其他語言寫好的模塊，也可以在其他語言中調用matlab生成的模塊。因為不熟悉matlab，所以我選擇的方式是用matlab生成動態鏈接庫，在c/c++中調用。這是Matlab Compiler提供的功能。
vs2010環境配置方式如下：
1、首先要在vs2010中設置所需頭文件的路徑，在我的電腦上是D:\Program Files\MATLAB\R2011a\extern\include；
2、指定庫文件路徑，D:\Program Files\MATLAB\R2011a\extern\lib\win32\microsoft；
3、在linker中添加附加依賴庫mclmcrrt.lib，Over。
然后在matlab中寫好需要的功能，第一次用matlab，代碼蠢笨，大家見笑了。
首先，要在matlab中實現稀疏矩陣的創建（文件名createsparse.m），其中m，n分別是行數和列數：
function A = createsparse(m, n)
A = sparse(m,n);
end
接下來要實現對稀疏矩陣中的非零元素賦值（文件名eleassign.m），其中i，j分別是行索引和列索引，value是元素的值：
function A = eleassign( A, i, j, value )
A(i,j)=value;
end
矩陣元素賦值完成后就是線性方程組的求解。這個函數是直接在mathworks官網上下載的，代碼比較長，這里只給出聲明（文件名lsmr.m）：
function [x]...
= lsmr(A, b, lambda, atol, btol, conlim, itnlim)
實際上原始的輸入和輸出參數要更多一些，由於我只用到一個輸出，即方程組的解向量，所以把其他輸出刪掉了。
現在有了這三個函數，分別對應三個.m文件，下面需要用這三個函數生成對應的動態鏈接庫，在matlab中使用如下命令：
mcc -W cpplib:SparseSolver -T link:lib createsparse.m eleassign.m lsmr.m
其中SparseSolver是庫文件名，-T link:lib 的作用是指定輸出文件類型為庫文件。
matlab busy一會兒以后，就會產生一系列文件，用到的是其中的三個：SparseSolver.h、SparseSolver.lib和SparseSolver.dll，把這三個文件拷貝到工程文件夾中。
其中SparseSolver.lib是導入庫，在vs2010的工程屬性中的附加依賴庫中添加這個文件。
SparseSolver.h是生成的c/c++函數聲明，比如上面實現的幾個函數，在頭文件中分別為：
extern LIB_SparseSolver_CPP_API void MW_CALL_CONV createsparse(int nargout,
mwArray& A,
const mwArray& m,
const mwArray& n);
extern LIB_SparseSolver_CPP_API void MW_CALL_CONV eleassign(int nargout,
mwArray& A,
const mwArray& A_in1,
const mwArray& i,
const mwArray& j,
const mwArray& value);
extern LIB_SparseSolver_CPP_API void MW_CALL_CONV lsmr(int nargout,
mwArray& x,
const mwArray& A,
const mwArray& b,
const mwArray& lambda,
const mwArray& atol,
const mwArray& btol,
const mwArray& conlim,
const mwArray& itnlim);
在這些函數中，第一個參數是輸出參數的個數。其他都是實際的輸入輸出參數。C++中與matlab互相傳數據都是用這種mwArray的數據類型。比如要C++中需要向matlab中傳遞一個向量，可能要這樣：
...
double vec[] = {1,2,3};
mwArray _vec(3,1,mxDOUBLE_CLASS,mxREAL);
_vec.SetData(vec,3);
...
mwArray _vec(3,1,...) 中，3是行數，1是列數。
有兩件事要注意，一是matlab中矩陣的存儲是列優先存儲，比如矩陣
在matlab中存儲的順序為1，1，1，1，-1，1。而在C++中我們習慣的存儲方式是{1,1,1,-1,1,1}，所以在用SetData拷貝數據的時候必須先轉換行列順序。
二是向matlab傳數據只能是以mwArray傳，即使只是一個普通數值，也必以mwArray的方式傳過去，比如這樣：
...
double somevalue[] = {1};
mwArray _somevalue(1,1,mxDOUBLE_CLASS,mxREAL);
_somevalue.SetData(somevalue,1);
...
當然mwArray也支持隨機訪問，比如：
...
mwArray _somevalue(1,1,mxDOUBLE_CLASS,mxREAL);
_somevalue(1,1) = 5.0f;
...
確實不太方便，不過比起尋找各種庫所花費的力氣，值了。
下面是一些功能測試代碼，很亂：
#include <iostream>
#include "SparseSolver.h"
using namespace std;
int main()
{
//mclmcrInitialize();
if (!mclInitializeApplication(NULL,0))
{
std::cerr << "could not initialize the application properly"
<< std::endl;
return -1;
}
if( !SparseSolverInitialize() )
{
std::cerr << "could not initialize the library properly"
<< std::endl;
return -1;
}
else
{
try
{
//parameters initialization
double lamda[] = {0.0f};
double atol[] = {0.0f};
double btol[] = {0.0f};
double conlim[] = {0.0f};
double itnlim[] = {50.0f};
mwArray _lamda(1,1,mxDOUBLE_CLASS,mxREAL);
//_lamda.SetData(lamda,1);
_lamda(1,1) = lamda[0];
mwArray _atol(1,1,mxDOUBLE_CLASS,mxREAL);
//_atol.SetData(atol,1);
_atol(1,1) = atol[0];
mwArray _btol(1,1,mxDOUBLE_CLASS,mxREAL);
//_btol.SetData(btol,1);
_btol(1,1) = btol[0];
mwArray _conlim(1,1,mxDOUBLE_CLASS,mxREAL);
//_conlim.SetData(conlim,1);
_conlim(1,1) = conlim[0];
mwArray _itnlim(1,1,mxDOUBLE_CLASS,mxREAL);
//_itnlim.SetData(itnlim,1);
_itnlim(1,1) = itnlim[0];
//data
double matA[] = {1.0,1.0,1.0,-1.0,1.0,1.0};
double vecb[] = {2.0,1.0,3.0};
mwArray _vecb(3,1,mxDOUBLE_CLASS,mxREAL);
_vecb.SetData(vecb,3);
//create matrix A
mwArray _matA;
double rows[] = {3.0f};
double cols[] = {2.0f};
mwArray _rows(1,1,mxDOUBLE_CLASS,mxREAL);
mwArray _cols(1,1,mxDOUBLE_CLASS,mxREAL);
_rows.SetData(rows,1);
_cols.SetData(cols,1);
createsparse(1,_matA,_rows,_cols);
cout<<_matA<<endl;
for(int i = 0; i < 3; i++)
{
for(int j = 0; j<2; j++)
{
double temp[1];
temp[0] = matA[i*2+j];
mwArray _temp(1,1,mxDOUBLE_CLASS,mxREAL);
_temp.SetData(temp,1);
double indexr[1],indexc[1];
indexr[0] = i+1;
indexc[0] = j+1;
// mwArray _indexr(1,1,mxDOUBLE_CLASS,mxREAL);
// mwArray _indexc(1,1,mxDOUBLE_CLASS,mxREAL);
// _indexr.SetData(indexr,1);
// _indexc.SetData(indexc,1);
// eleassign(1,_matA,_matA,_indexr,_indexc,_temp);
_matA(i+1,j+1) = matA[i*2+j];
}
}
cout<<_matA<<endl;
int cnt = 0;
for(int i = 0; i<3;i++)
{
for(int j= 0;j<2;j++)
_matA(i+1,j+1) = cnt++;
}
cout<<_matA<<endl;
for(int i=0;i<3;i++)
{
for(int j=0;j<2;j++)
{
double temp = _matA(i+1,j+1);
if(temp == 0)
{
cout<<"zero pos:"<<endl;
cout<<i+1<<" "<<j+1<<endl;
}
}
}
mwArray _matB;
double rowsb[] = {4.0f};
double colsb[] = {4.0f};
mwArray _rowsb(1,1,mxDOUBLE_CLASS,mxREAL);
mwArray _colsb(1,1,mxDOUBLE_CLASS,mxREAL);
_rowsb.SetData(rowsb,1);
_colsb.SetData(colsb,1);
createsparse(1,_matB,_rowsb,_colsb);
double rbegin[] = {1.0f};
double rend[] = {3.0f};
double lbegin[] = {1.0f};
double lend[] = {2.0f};
mwArray _rbegin(1,1,mxDOUBLE_CLASS,mxREAL);
mwArray _rend(1,1,mxDOUBLE_CLASS,mxREAL);
mwArray _lbegin(1,1,mxDOUBLE_CLASS,mxREAL);
mwArray _lend(1,1,mxDOUBLE_CLASS,mxREAL);
_rbegin.SetData(rbegin,1);
_rend.SetData(rend,1);
_lbegin.SetData(lbegin,1);
_lend.SetData(lend,1);
cout<<"B"<<endl;
cout<<_matB<<endl;
submat(1,_matB,_matB,_rbegin,_rend,_lbegin,_lend,_matA);
cout<<"A:"<<endl;
cout<<_matA<<endl;
cout<<"B"<<endl;
cout<<_matB<<endl;
double temp[1];
double ok = 2;
temp[0] = ok;
mwArray _temp(1,1,mxDOUBLE_CLASS,mxREAL);
_temp.SetData(temp,1);
matmult(1,_matB,_matB,_temp);
cout<<"B after mult"<<endl;
cout<<_matB<<endl;
cout<<_matB(4,4)<<endl;
_matB(4,4) = 5.0;
cout<<_matB(4,4)<<endl;
mwArray _x;
lsmr(1,_x,_matA,_vecb,_lamda,_atol,_btol,_conlim,_itnlim);
cout<<_x<<endl;
}
catch (const mwException& e)
{
cerr << e.what() << endl;
system("pause");
return -2;
}
catch (...)
{
cerr << "Unexpected error thrown" << endl;
system("pause");
return -3;
}
// Call the application and library termination routine
SparseSolverTerminate();
}
mclTerminateApplication();
system("pause");
return 0;
}
最后總結一下，說matlab運算速度慢，效率低的，可能要么是大牛程序員，要么是人雲亦雲。低水平的coder，像我這種，自己寫的數值分析算法肯定不如matlab實現的。所以，對於一般程序員來說，用matlab做一個后台的計算工具是一件很爽的事情。另外matlab本身有成熟的交流平台，file exchange有大量的代碼資源可以下載，也能節約不少搜集資源的時間。
就寫到這里，關於matlab計算效率的問題，等把骨架提取的代碼改完后再補充吧。