在VisualStudio2005中求解实对称矩阵特征值算法包调查

这几天为了帮MM解决一个技术问题，在网上查了很多关于如何使用C/C++算法包计算大型实对称矩阵特征值的资料，这里小结一下。

开发平台：win32, Visual Studio 2005
待解决问题：在C/C++代码中求解大约2000*2000的实对称稀疏矩阵的特征值和特征向量。
Matlab中令人吃惊的效率：使用稀疏矩阵的存储方式，调用eigs()函数可以在几秒内解决该问题。(多么令人畏惧的效率啊！对比以下解决方案可知。)

C/C++解决方案：
(1) 使用OpenCV的cvEigenVV()或者cvSVD()求解，结果大约需要花费4分钟。
(2) 使用Matcom45转换Matlab代码的eig()函数为C/C++代码(注意，eig是矩阵在稠密存储方式下的特征值求解，而eigs是求稀疏矩阵特征值，两者效率差别很大)，因为eigs()函数无法用Matcom45转换(可怜的MM，在老板的催促下花了好多精力尝试这个……)。结果求解时间比OpenCV费时还多……
(3) 开始上网找。首先是MM找到的一篇介绍SLEPc的survey (A Survey of Software for Sparse Eigenvalue Problems)，里面提到了不少解决稀疏矩阵特征值问题的软件。
其中提到的部分算法包如下：
• Fortran77 libraries: ARPACK, BLZPACK, PDACG, QMRPACK, NA18 and SRRIT. Note that these libraries can usually be called from C/C++ code, if appropriate calling conventions are used.
• C libraries: PRIMME, JDBSYM and MPB. These can be called from C++ code without problems. primme also provides a Fortran77 interface.
• Fortran90 libraries: SPAM and TRLAN. These can be called from C/C++ and Fortran77.
• C++ libraries: ANASAZI and IETL.
可惜大多数算法包的开发环境貌似都是Unix环境，在Windows平台下使用起来不是很方便，有些甚至是噩梦。所谓顺藤摸瓜，就从这里列到的几个C/C++算法包开始查起，调查调查每个算法包是否好用。
(4) SLEPc(the Scalable Library for Eigenvalue Problem Computations)：使用C语言编写的算法包。需要安装PETSc。因为是Unix/Linux平台下开发的，看起来在windows下不是很好编译，所以没有深入尝试。
(5) TRILINOS/ANASAZI：使用C++开发的基于面向对象框架的算法库。貌似在windows平台不太容易使用，没有进行深入尝试。
(6) PRIMME(PReconditioned Iterative MultiMethod Eigensolver)：使用C语言编写的算法包，专门针对对称阵或者Hermitian矩阵的特征值求解，速度应该不错。这里(State-of-the-art numerical solution of large Hermitian eigenvalue problems)还有一篇作者自己的广告PPT，声称PRIMME是一个state-of-the-art的特征值求解器。在Win32平台下用VS2005直接将其源代码编译，只需要修改几个小地方，可以成功生成静态链接库。使用该算法包时还需要BLAS和LAPACK(Win32平台下可以使用CLAPACK-VisualStudio版的静态链接库)。代码使用方法可以参见上面的PPT里的例子，貌似比较简单，但是需要搞清楚其矩阵输入的格式。结论：没时间摸索其矩阵输入的格式，就此放弃，但是个人感觉该库还是值得一试，毕竟专门针对对称阵，算法做了有针对性的优化，还提供了很多种可选的解法。
(7) IETL(The Iterative Eigensolver Template Library)：使用C++写的一个模板库，使用了boost库。使用其中的部分算法(如Lanczos算法)时需要BLAS和LAPACK库。该库的主页上说，LAPACK提供了对稠密矩阵求解特征值的最好的算法，但是对于稀疏矩阵来讲，如果只需要求解其部分特征值而不是全部的话，使用迭代算法会更快(这也说明了为什么matlab中eigs和eig算法效率差距为何如此巨大)。该库在C++中比较易用，环境也比较容易配置，基本只需要安装boost库即可。结论：测试其例子程序提供的使用power算法求解特征值，2000阶矩阵耗时3分钟左右。
(8) ARPACK：该算法包貌似受到的评价很不错，是使用Fortran语言写的解决大型稀疏矩阵特征值问题的算法包。需要用到外部库BLAS和LAPACK。这里有一篇介绍在windows平台下使用ARPACK的方法，只是我没有试。这里有一些为windows平台编译好的库，我试了一下，因为貌似依赖外部库比较乱，没有成功。结论：该算法包应该不错，只是使用起来比较麻烦。
(9) SuperLU：在ARPACK中提到了该算法包，貌似评价不错。查了一下，该算法包是用于求解大型稀疏线性方程组的，是由C语言编写的。只是不知道有没有直接求解特征值的功能，又因为其在windows平台下不太容易编译，所以没有进行深入摸索。
(10)Eigen：另一个用C++写的模板库，用于解决常见的线性代数问题，当然由其名字可以看出来肯定包括了求解特征值的功能。从其网站上提供的tutorial中可以看到，该库很方便使用，提供的矩阵表示方法非常简洁方便，与IETL有些类似。结论：使用其提供的EigenSolver和SVD试了一下，速度比较慢，求解了一个500阶的矩阵特征值花费了20多分钟(GOD!)，网上看到一篇拿此库和LAPACK同时求解对称阵特征值进行对比的文章，看这里，好像说Eigen是about a factor of 3-4 slower than LAPACK，看来是要慢不少。结论：好用，但是太慢。
(11)LAPACK/CLAPACK：大名鼎鼎的主角登场了，前面几乎所有的算法包都提到了LAPACK，可见其影响力之大。LAPACK是用Fortran编写的算法库，顾名思义，Linear Algebra PACKage，是为了解决通用的线性代数问题的。另外必须要提的算法包是BLAS(Basic Linear Algebra Subprograms)，其实LAPACK底层是使用了BLAS库的。不少计算机厂商都提供了针对不同处理器进行了优化的BLAS/LAPACK算法包，例如Intel的MKL(Math Kernel Library，很不幸是收费的)，AMD的ACML等。在Matlab的bin目录里可以发现MKL和ACML动态链接库的踪影，所以由此推断，Matlab底层应该也是使用了BLAS/LAPACK库的。闲话少说，如何在windows下使用LAPACK解决问题呢？这里介绍了在windows下使用LAPACK的方法和很多相关资源的下载。CLAPACK是使用f2c工具将LAPACK 的Fortran代码转换成C语言代码的C语言算法包。这里提供了在VS2005下的CLAPACK工程及已经编译好的静态链接库，这里有一篇中文写的如何在VC中调用CLAPACK。LAPACK的函数命名规则(naming scheme)比较奇怪，这里可以找到相关介绍，这里有一篇中文介绍。我使用CLAPACK测试了几个函数，其中一个是DSYEVD，命名含义如下：D表示双精度，SY表示对称阵，EV表示特征值问题，D表示求解算法中使用了divide-and-conquer。所以DSYEVD就是求解对称阵特征值的函数。调用该函数的C语言代码参考了Intel MKL库里提供的一个例程，在这里可以找到。结论：我修改后的代码如下，测试2000*2000的矩阵求解特征值，耗时不超过1分钟(果然比较强悍，注意，这里还是求解的稠密矩阵！)。

/*
使用CLAPACK中dsyevd_函数求解对称矩阵特征值的例程。
CLAPACK库：http://www.netlib.org/clapack/CLAPACK-3.1.1-VisualStudio.zip
参考：http://software.intel.com/sites/products/documentation/hpc/mkl/lapack/mkl_lapack_examples/dsyevd_ex.c.htm
修改：MulinB
日期：2010-04-16
*/
#include <stdio.h>
#include "f2c.h"
#include "clapack.h"

/* Auxiliary routines prototypes */
extern void print_matrix( char* desc, int m, int n, double* a, int lda );

/* Parameters */
#define N 2000
#define LDA N

/* Main program */
int main() {
/* Locals */
int xx, yy; //for loop
int n = N, lda = LDA, info, lwork, liwork;
char JOBZ = 'V';
char UPLO = 'U';
int iwkopt;
int* iwork;
double wkopt;
double* work;
double* w;
double* a;

a = (double*)malloc(sizeof(double)*LDA*N);
for (xx=0; xx<LDA; xx++)
{
for (yy=0; yy<N; yy++)
{
a[xx*LDA + yy] = xx+yy;
}
}

w =  (double*)malloc(sizeof(double)*N);

/* Executable statements */
printf( " DSYEVD Example Program Results: (eigen of sym matrix: %d * %d)/n", N, N );

lwork = -1;
liwork = -1;
/* Query and allocate the optimal workspace */
dsyevd_( &JOBZ, &UPLO, &n, a, &lda, w, &wkopt, &lwork, &iwkopt, &liwork, &info );

lwork = (int)wkopt;
work = (double*)malloc( lwork*sizeof(double) );
liwork = iwkopt;
iwork = (int*)malloc( liwork*sizeof(int) );
/* Solve eigenproblem */
dsyevd_( &JOBZ, &UPLO, &n, a, &lda, w, work, &lwork, iwork, &liwork, &info );

/* Check for convergence */
if( info > 0 ) {
printf( "The algorithm failed to compute eigenvalues./n" );
exit( 1 );
}
else
{
printf("INFO=%d (succeed!)/n", info);
}

/* Print eigenvalues */
//print_matrix( "Eigenvalues", 1, n, w, 1 );

///* Print eigenvectors */
//print_matrix( "Eigenvectors (stored columnwise)", n, n, a, lda );

/* Free workspace */
free( (void*)iwork );
free( (void*)work );

return 0;
} /* End of DSYEVD Example */

/* Auxiliary routine: printing a matrix */
void print_matrix( char* desc, int m, int n, double* a, int lda ) {
int i, j;
printf( "/n %s/n", desc );
for( i = 0; i < m; i++ ) {
for( j = 0; j < n; j++ ) printf( " %6.2f", a[i+j*lda] );
printf( "/n" );
}
}