在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" ); }}