作为一个遥感图像处理方向研究生，笔者在编程的时候免不得要和矩阵运算打交道，各种矩阵乘法协方差特征值。虽然在自己的研究实验中更喜欢用MATLAB来作为工具，但是在完成工程项目编软件的时候必须要用C++来作为编程语言。在C++中并不能像MATLAB一样，非常方便的进行矩阵运算。当然也可以自己编一些列矩阵运算的函数（笔者以前就干过这事），不过这个运算效率实在无法满足要求。
目前国际上提供了很多矩阵运算线性代数的库，这些库都进行了基于CPU和多核的运算优化，自然比我们自己做的这些函数效率要高的多，而且很多库都有着非常友好的使用方式，这就好比有了又快又舒服汽车自然也就不用再用自己的双腿走到累死。经过了多次折腾，笔者终于大概了解了几个比较常用的矩阵库的安装和使用方式，在这里做一个简要介绍，希望能为其他人提供一个方便。
笔者最早是使用的别人给的MATCOM库，这个库比较古老、不开源而且早就停止维护，唯一的优点就是可以在vc6.0上使用···最近又开始在找更加优质的矩阵运算库，就发现了Armadillo，并且在同学的推荐下用了一下Eigen，发现在vs2010的环境下，这两个库确实比MATCOM要好上N倍。以下就是各个运算库的介绍和安装，本文重点介绍的是Armadillo+OpenBLAS的组合，这个组合是我目前使用的库中最简单效率又最高的。
1、MATCOM
MATCOM矩阵库在国内用的貌似还不少，笔者就是从别的同学那里得到的这个库。根据度娘的介绍，MATCOM其实是mathtools公司提供讲MATLAB代码移植到c++上的一个工具，具体两者的关系笔者搞不太清，但是MATCOM在4.5版本之后据说就被合并到MATLAB中，没有再更新和维护。
MATCOM有两大优点，一是和MATLAB的使用方式非常像，基本函数都叫一个名字，所以可以非常容易地将MATLAB代码移植到C++上；二是可以在VC6上使用，这一点非常重要···因为目前大多数的库都已经不再支持VC6了。不过MATCOM的缺点也非常多，首先MATCOM并不开源，这一点在完成工程项目编软件的时候很麻烦；其次MATCOM在MathTools公司被合并之后就再也没有更新和维护，因此在代码的速度和健康程度上已经是落后于时代了。所以笔者认为除非是必须在VC6下进行编程，不然不建议再使用MATCOM这个矩阵库。
关于MATCOM的具体应用，可以参见http://blog.sina.com.cn/s/blog_60f8483a0100gay7.html这个博客，讲解的非常清楚。在度娘中搜索，也可以找到很多关于MATCOM的文章。不过笔者要提醒一点，在笔者的使用中发现有一些文章和说明中提到的函数和真正使用起来并不一样，这可能是版本差异，所以需要自己动手尝试。
笔者所获得的MATCOM一共只有4个文件，是编译好了的v4501v.lib,v4501v.dll，ago4501.dll和matlib.h文件。这几个文件读者可以在CSDN等网站上搜索下载。
安装方式很简单，只需要在工程中添加上v4501v.lib，在代码头包含#include "matlib.h"，并且把v4501v.dll，ago4501.dll放到应用文件，例如DEBUG或者RELEASE文件夹中，就可以使用。
2、Armadillo
Armadillo是一个开源的矩阵运算库，这个库可以在开源中国下载到
Armadillo的主页也非常友好
这个库其实是对底层的矩阵运算库（例如BLAS，LAPACK）的封装，也就是说，这个库本身的计算能力其实不强，但是这个库可以配合开源的底层库来使用，例如最基本的BLAS+LAPACK或者改进后的OpenBLAS,ACML还有强大的MKL。Armadillo的速度主要就是由底层库决定，其代码使用可以说是非常非常的简单，基本和MATLAB的代码没有什么差别，总之笔者是用的非常爽。
Armadillo的更新很频繁，这正是它还很有生命力的证据。
Armadillo可以从开源中国的连接或者是其主页的连接上下载，其安装说起来就是几步
1-将下载后的文件解压到某一个文件夹，例如D:\Program Files\Armadillo\armadillo-3.910.0；
2-将D:\Program Files\Armadillo\armadillo-3.910.0\include目录包含到在VS工程的文件目录中；
3-如果想要使用BLAS和LAPACK等底层库的加速，需要在将D:\Program Files\Armadillo\armadillo-3.910.0\include\armadillo_bits\config.hpp中的取消一下两行代码的注释
#define ARMA_USE_LAPACK
#define ARMA_USE_BLAS
笔者推荐取消注释，这样才能发挥Armadillo的能力。
4-安装就基本完成了，以上安装可以参考以下两个网址
在度娘上可以搜索到很多Armadillo的安装文章，不过有一些是基于linux的，所以大家可以找针对windows的文章。
安装完毕后需要实践一下才能检验是否知道怎么使用，在压缩文件目录中，有一个examples文件夹，下面有两个例子，安装结束后请运行一下两个例子。
5-因为在Armadillo中使用到了BLAS和LAPACK，在examples\lib_win32中提供了两个已经编译好了的LIB和DLL，在使用到Armadillo的工程中加入两个LIB文件，并把DLL文件拷到执行文件的目录下。在这个文件夹中提供的两个库lapack_win32_MT和blas_win32_MT是提前编译好了的，非常好用，笔者测试后发现比网上自己下载的LAPACK和BLAS要快很多，而且还不需要其他的DLL支持。如果读者不想麻烦，就可以直接用这两个库。
6-请把压缩文件中的Readme.txt好好阅读一遍，这个帮助文件其实就讲的很详细了。还有docs.html，是一份非常详细的使用说明。
下载最新的LAPACK
可能有读者想使用自己下载的LAPACK库，笔者推荐在以下网址下载
读者可以直接在这个网址搜索“Using LAPACKE subroutines in a Visual Studio C/C++ Project”，在这个小节下面有一个编译好的LAPACKE examples，里面有编译好的LIB和DLL，可以拿来直接使用。
不过需要注意的是，网上目前大多数库，包括LAPACK和OpenBLAS，都是使用MinGW进行编译的，所以为了能顺利使用这些库，请读者先去安装一个MinGW吧。请在安装的时候选择C，C++和FORTRAN，这些编译环境都是经常需要的。
读者自己下载LAPACK之后，就可以使用LIB和DLL替换原来Armadillo提供的LIB和DLL，不过因为网上下载的LAPACK是使用MinGW编译的，需要相应的编译环境才能使用，单独使用会提示缺少必要的DLL。有两个办法解决这个问题，一是安装MinGW，并且在系统Path中进行配置，具体请度娘；二是将MinGW安装目录bin文件夹中的相应DLL拷到执行文件目录下，缺什么考什么。笔者推荐第二种方式，因为编好的软件一般都会在其他电脑上运行，不能要求所有的电脑都有MinGW环境。一般需要的DLL是libgcc_s_dw2-1.dll，libgfortran-3.dll和libquadmath-0.dll三个。
使用OpenBLAS
在Armadillo的Readme.txt中，一共推荐了三个库，并说“The OpenBLAS, MKL and ACML libraries are generally the fastest.”其中笔者推荐要用的，就是OpenBLAS。OpenBLAS是一个开源的库，这一点正好能够满足Armadillo联合使用工程项目编程的需要。下面就是OpenBLAS的主页
OpenBLAS有很多优点，计算效率高，开源免费，不过最最重要的一点，这是中国人写的！OpenBLAS是由中科院软件研究所的张先轶开发的一个开源计算库。为什么说这是一个很重要的优点，倒不是说笔者呼吁支持国货，而是如果有问题可以直接用中文问开发者！而且还是一个很热心的开发者！还有什么比这更大的优点。
下面就是作者Xianyi Zhang在开源中国开的博客
笔者在有问题的时候晚上9:29留言，开发者9:35就回答了，第二天更是写了一个“在Windows上Armadillo如何使用OpenBLAS”的博客，笔者当时是感动得泪流满面啊。
这就是开发者提供的如何使用Armadillo+OpenBLAS，里面步骤写的很详细，读者可以按照这个博客来做
唯一要补充一点就是这个库运行也需要MinGW环境，上面关于LAPACK的文字中有提到如何解决。只不过对于OpenBLAS只需要libgfortran-3.dll和libquadmath-0.dll。
在使用OpenBLAS的时候，只需要将LIB文件添加到工程中，再把DLL拷到执行文件目录下就行了。在笔者安装的时候拷到一些英文的安装方法说还需要BLAS或者LAPACK，其实OpenBLAS都已经包含了，除了MinGW环境不需要额外的库。
OpenBLAS已经编译好的文件可以在主页上直接下载到，也可以下载源代码自己编译，关于如何编译作者已经在他的博客中说道了。
Armadillo+OpenBLAS是笔者最推荐的组合。
其他库MKL ACML ATLAS
除了上面提到的BLAS+LAPACK和OpenBLAS，其他著名的库有Intel MKL，AMD ACML以及开源的ATLAS。其中MKL是收费的，ACML和ATLAS是免费的。MKL听说是效率最高的库，但是在Windows下要收费，Linux免费，这是赤裸裸地鄙视Windows程序员的节奏吗。MATLAB就是使用的MKL和ATLAS。当然，有心人总能在网上找到免费的MKL，不过因为只能自己用，所以笔者也就没认真研究。曾经下载到一个MKL的库试过，发现不如OpenBLAS这样容易加到Armadillo，老是出错，也就没再深入过。ACML是免费的，不过貌似需要编译，笔者找到合心的库之后就不想麻烦了，就没弄，如果有兴趣的读者可以自己去网上下载。ATLAS因为是需要编译，笔者也懒了没弄。
3.Eigen
Eigen也是一个开源的矩阵运算库，笔者在同学的推荐下尝试了一下这个库。据笔者同学说，在Windows下Eigen比OpenBLAS效率更高，而且在下面的这个测评中，Eigen比Armadillo要快很多，但是笔者猜这里的Armadillo没有使用OpenBLAS.
Eigen不是以动态库的形式，而是以头文件的方式使用，因此不需要再包含其他的DLL，应该算比较方便的。Eigen的安装方式很简单，就是将下载的压缩包直接解压，例如D:\Program Files\Eigen，然后在VS工程中直接包含路径D:\Program Files\Eigen\Eigen，也就是根目录下的Eigen文件夹就行了。
笔者猜测Eigen可能可以在VC6下使用，如果感兴趣的读者可以自己去尝试一下。
Eigen最大的缺点就是和MATLAB的使用方法差得太多，笔者用起来特别费劲。不过Eigen属于大多数算法都用模板函数，可能有人会觉得用起来更方便。
在这篇博文中，笔者主要进行了针对MATCOM，Armadillo+LAPACK_MT（Armadillo提供的LAPACK和BLAS），Armadillo+LAPACK（网上下载的，需要MinGW环境），Armadillo+OpenBLAS，Eigen以及MATLAB的测评。测评内容大概是一个PCA变换，所有的程序都进行相同的计算内容进行计时。程序都是编译成MT的Release版本进行计算，使用的计时工具是Timer （http://www.songho.ca/misc/timer/timer.html）。
MATCOM的代码
///
#include
using namespace std;
//包含MATCOM头文件
#include "matlib.h"
#include "Timer.h"
//包含GDAL头文件，GDAL是一个用来读取影像的库，搞不懂的就记住这个东西是用来读取一个大矩阵的
#include "gdal_priv.h"
int main()
{
winaxes((HWND) NULL);
initM(MATCOM_VERSION);
Timer time;
cout<<"Hello World!"<<endl;
cout<<"Hello VS2010!"<<endl;
cout<<"Hello Armadillo!"<<endl;
GDALDataset * pGDALDataIn;
GDALAllRegister();
//这里读取了一个影像，大小是1061*1078*4
const char * cPath="D:\\2002qb_wuda_ms.img";
pGDALDataIn=(GDALDataset *)GDALOpen(cPath,GA_ReadOnly);
if (pGDALDataIn==NULL)
{
cout<<"输入文件打开失败！"<<endl;
system("pause");
return 1;
}
int nHeight, nWidth, nBand;
nHeight=pGDALDataIn->GetRasterYSize();
nWidth=pGDALDataIn->GetRasterXSize();
nBand=pGDALDataIn->GetRasterCount();
double * pData;
pData=new double [nHeight*nWidth*nBand];
pGDALDataIn->RasterIO(GF_Read,0,0,nWidth,nHeight,pData,nWidth,nHeight,GDT_Float64,nBand,NULL,0,0,0);
int nImageSize=nHeight*nWidth;
//上面都不用管，总之就是把数据读取到了一个double型的数组里
Mm mImage =zeros(nImageSize,nBand);
double *pmImageData=mxGetPr(&mImage);
//这是用拷内存的方式给MATCOM矩阵赋值，这么做应该是最快的赋值方式了，但是要注意
//矩阵和数组都是double型的才可以
memcpy(pmImageData,pData,nImageSize*nBand*sizeof(double));
mImage=transpose(mImage);
//开始计时
time.start();
Mm mMeanImg, mMean;
mMean=mean(mImage,2);
mMeanImg=mImage-repmat(mMean,1,nImageSize);
Mm mCovImg;
mCovImg=mMeanImg*transpose(mMeanImg);
mCovImg=mCovImg/(nImageSize-1);
Mm mTransW;
Mm mEigVal;
i_o_t iot={0,0};
eig(mCovImg, iot, mTransW, mEigVal);
Mm mInv=inv(mCovImg);
//计时结束
time.stop();
cout<<"The time is "<<time.getElapsedTimeInMilliSec
delete pGDALDataIn;
delete [] pData;
system("pause");
}
///
Armadillo代码
因为Armadillo代码都是一样的，不同的只不过是加入的底层库不同。
///
#include
#include "armadillo"
//命名空间
using namespace arma;
using namespace std;
#include "Timer.h"
#include "gdal_priv.h"
int main()
{
Timer time;
cout<<"Hello World!"<<endl;
cout<<"Hello VS2010!"<<endl;
cout<<"Hello Armadillo!"<<endl;
GDALDataset * pGDALDataIn;
GDALAllRegister();
const char * cPath="D:\\2002qb_wuda_ms.img";
pGDALDataIn=(GDALDataset *)GDALOpen(cPath,GA_ReadOnly);
if (pGDALDataIn==NULL)
{
cout<<"输入文件打开失败！"<<endl;
system("pause");
return 1;
}
int nHeight, nWidth, nBand;
nHeight=pGDALDataIn->GetRasterYSize();
nWidth=pGDALDataIn->GetRasterXSize();
nBand=pGDALDataIn->GetRasterCount();
double * pData;
pData=new double [nHeight*nWidth*nBand];
pGDALDataIn->RasterIO(GF_Read,0,0,nWidth,nHeight,pData,nWidth,nHeight,GDT_Float64,nBand,NULL,0,0,0);
mat Imgmat=zeros(nHeight*nWidth,nBand);
//同样是使用拷内存的方式赋值，好不容易才找到的
double *pmImageData=Imgmat.memptr();
memcpy(pmImageData,pData,nHeight*nWidth*nBand*sizeof(double));
Imgmat=trans(Imgmat);
cout<<Imgmat.n_rows<<" "<<Imgmat.n_cols<<endl;
//开始计时
time.start();
mat mMean, mMeanImg;
mMean=mean(Imgmat,1);
mMeanImg=Imgmat-repmat(mMean,1,nHeight*nWidth);
mat mCov1=mMeanImg*trans(mMeanImg);
mCov1=mCov1/(nHeight*nWidth-1);
vec eigval;
mat eigvec;
eig_sym(eigval,eigvec,mCov1);
mat mInvCov=inv(mCov1);
//计时结束
time.stop();
cout<<eigvec<<endl;
cout<<"The time is "<<time.getElapsedTimeInMilliSec
delete pGDALDataIn;
delete [] pData;
system("pause");
}
///
Eigen代码
///
#include
//包含Eigen头文件
#include
#include
using namespace std;
using namespace Eigen;
#include "Timer.h"
#include "gdal_priv.h"
int main()
{
Timer time;
cout<<"Hello World!"<<endl;
cout<<"Hello VS2010!"<<endl;
cout<<"Hello Armadillo!"<<endl;
GDALDataset * pGDALDataIn;
GDALAllRegister();
const char * cPath="D:\\2002qb_wuda_ms.img";
pGDALDataIn=(GDALDataset *)GDALOpen(cPath,GA_ReadOnly);
if (pGDALDataIn==NULL)
{
cout<<"输入文件打开失败！"<<endl;
system("pause");
return 1;
}
int nHeight, nWidth, nBand;
nHeight=pGDALDataIn->GetRasterYSize();
nWidth=pGDALDataIn->GetRasterXSize();
nBand=pGDALDataIn->GetRasterCount();
double * pData;
pData=new double [nHeight*nWidth*nBand];
pGDALDataIn->RasterIO(GF_Read,0,0,nWidth,nHeight,pData,nWidth,nHeight,GDT_Float64,nBand,NULL,0,0,0);
const int nImageSize=nHeight*nWidth;
//用数组的复制方式
Map mImage(pData,nImageSize,nBand);
cout<<mImage(0,1)<<" "<<mImage(1,0)<<endl;
//定义double型矩阵
MatrixXd mTImage;
mTImage=mImage.transpose();
//开始计时
time.start();
MatrixXd mMeanImg, mMean;
mMean=mTImage.rowwise().mean();
mMeanImg=mTImage-mMean.replicate(1,nImageSize);
MatrixXd mCovImg;
mCovImg=mMeanImg*mMeanImg.transpose();
mCovImg=mCovImg/(nImageSize-1);
EigenSolver es;
MatrixXcd mTranW, mEigVal;
es.compute(mCovImg,false);
//这一步赋值运行老是出错，我也不知道错在哪里，不过不影响
//mTranW=es.eigenvectors();
mEigVal=es.eigenvalues().asDiagonal();
MatrixXd mInv=mCovImg.inverse();
//停止计时
time.stop();
cout<<mTranW<<endl;
cout<<"The time is "<<time.getElapsedTimeInMilliSec
delete pGDALDataIn;
delete [] pData;
system("pause");
}
///
从上面的三段代码就可以看出来，Eigen和MATLAB的使用差得非常远，而Armadillo和MATCOM都和MATLAB非常像。笔者用Eigen用得是很痛苦，可能是用习惯MATLAB的关系。
经过测评，以上这段代码，从计时开始到计时结束，各个库的表现是
MATLAB 61.7101ms
Armadillo+OpenBLAS 73.5832ms
Armadillo+LAPACK_MT 86.5733ms
Eigen 99.3488ms
Armadillo+LAPACK 167.915ms
MATCOM 269.978ms
笔者只大概做了一个测试，从测试结果上看，MATLAB最快，而Armadillo+OpenBLAS次之，Armadillo提供的LAPACK_MT比网上下载的LAPACK要快很多，不知道什么原因，MATCOM最慢。
在网上还有一个测评
这个测评中Armadillo使用了OpenBLAS。其结果也是Armadillo+OpenBLAS最好。
笔者只是大概做了个测评，如果有兴趣，可以自己在做一个更加全面的测评。本次测评使用的环境是Windows7+VS2010。
2万+
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
