最近仿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提供的功能。
vs20008环境配置方式如下：
1、首先要在vs2008中设置所需头文件的路径，在我的电脑上是D:\Program Files\MATLAB\R2009a\extern\include；
2、指定库文件路径，D:\Program Files\MATLAB\R2009a\extern\lib\win32\microsoft；
3、在linker中添加附加依赖库mclmcrrt.lib，Over。
然后在matlab中写好需要的功能：
function x = SparseSolve(r,c,v,m,n,b)%use this function to solve sparse matrix
S = sparse(r,c,v,m,n);
x = S\b;
end
由此生成一个.m文件，下面需要用这个函数生成对应的动态链接库，在matlab中使用如下命令：
mcc -W cpplib:SparseSolver -T link:lib SparseSolve.m
其中SparseSolve是库文件名，-T link:lib 的作用是指定输出文件类型为库文件。
matlab busy一会儿以后，就会产生一系列文件，用到的是其中的三个：SparseSolve.h、SparseSolve.lib和SparseSolve.dll，把这三个文件拷贝到工程文件夹中。
其中SparseSolve.lib是导入库，在vs2008的工程属性中的附加依赖库中添加这个文件。
SparseSolve.h是生成的c/c++函数声明:
......
extern LIB_SparseSolve_CPP_API void MW_CALL_CONV SparseSolve(int nargout
, mwArray& x
, const mwArray& r
, const mwArray& c
, const mwArray& v
, const mwArray& m
, const mwArray& n
, const mwArray& b); .....
在这些函数中，第一个参数是输出参数的个数。其他都是实际的输入输出参数。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 "stdafx.h"
#include "SparseSolve.h"
#include <iostream>
using namespace std;
int _tmain(int argc, _TCHAR* argv[])
{
mclInitializeApplication(NULL,0);//
SparseSolveInitialize();//这两行代码必须有，否则会出现很多问题，本人也是调试之后发现问题之后，纠结了老半天，添加之后才解决的。
int* row = new int[8];
*row = 1;
*(row+1) = 2;
*(row+2) = 3;
*(row+3) = 1;
*(row+4) = 2;
*(row+5) = 3;
*(row+6) = 4;
*(row+7) = 4;
int* col = new int[8];
*col = 1;
*(col+1) = 1;
*(col+2) = 1;
*(col+3) = 2;
*(col+4) = 2;
*(col+5) = 3;
*(col+6) = 3;
*(col+7) = 4;
mxDouble* val = new mxDouble[8];
*val = 1;
*(val+1) = 2;
*(val+2) = 3;
*(val+3) = 3;
*(val+4) = 4;
*(val+5) = 1;
*(val+6) = 2;
*(val+7) = 1;
mxDouble* b = new mxDouble[4];
*b = 1;
*(b+1) = 0;
*(b+2) = 0;
*(b+3) = 1;
mwArray rowArray(8,1,mxDOUBLE_CLASS);
rowArray.SetData(row,8);
mwArray colArray(8,1,mxDOUBLE_CLASS);
colArray.SetData(col,8);
mwArray valArray(8,1,mxDOUBLE_CLASS);
valArray.SetData(val,8);
mwArray bArray(4,1,mxDOUBLE_CLASS);
bArray.SetData(b,4);
mwArray rArray(4,1,mxDOUBLE_CLASS);
int* m = new int[1];
int* n = new int[1];
*m = 4;
*n = 4;
mwArray mArray(1,1,mxDOUBLE_CLASS);
mArray.SetData(m,1);
mwArray nArray(1,1,mxDOUBLE_CLASS);
nArray.SetData(n,1);
cout << "row" <<endl;
cout << rowArray.ToString() << endl;
cout << "col" <<endl;
cout << colArray.ToString() << endl;
cout << "val" <<endl;
cout << valArray.ToString() << endl;
cout << "b" <<endl;
cout << bArray.ToString() << endl;
SparseSolve(1,rArray,rowArray,colArray,valArray,mArray,nArray,bArray);
mxDouble* result = new mxDouble[4];
rArray.GetData(result,4);
cout << "result:" << endl;
for(int i = 0; i < 4; i++)
{
cout << result[i] << endl;
}
delete[] row;
delete[] col;
delete[] val;
delete[] b;
delete[] m;
delete[] n;
SparseSolveTerminate();
mclTerminateApplication();
return 0;
}
//代码没有初始化错误判断和异常处理等，实际使用时，需要注意
最后总结一下，说matlab运算速度慢，效率低的，可能要么是大牛程序员，要么是人云亦云。低水平的coder，像我这种，自己写的数值分析算法肯定不如matlab实现的。所以，对于一般程序员来说，用matlab做一个后台的计算工具是一件很爽的事情。另外matlab本身有成熟的交流平台，file exchange有大量的代码资源可以下载，也能节约不少搜集资源的时间。