在VisualStudio2005中求解实对称矩阵特征值算法包调查
来源:互联网 发布:东莞cnc编程教学 编辑:程序博客网 时间:2024/06/05 07:07
这几天为了帮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" );
- }
- }
- 在VisualStudio2005中求解实对称矩阵特征值算法包调查
- 在VisualStudio2005中求解实对称矩阵特征值算法包调查
- 实对称矩阵特征值求解算法:Jacobi行循环法
- mkl中使用dsyevd求解实对称矩阵的所有特征向量和特征值
- Householder+QR法求解实对称矩阵全部特征值和特征向量(VC++)
- 基于 intel MKL 的对称矩阵特征值求解器
- 计算实对称矩阵特征值和特征向量
- 利用QR算法求解矩阵的特征值和特征向量
- 雅可比旋转求解对称二维矩阵的特征值和特征向量
- mkl中dgeev求解矩阵的特征向量和特征值
- 矩阵特征值的求解例子
- 矩阵特征值和特征向量求解——特征值分解
- 正定对称矩阵快速求解
- 幂法求解矩阵特征值及特征向量
- 神经网络求解 矩阵特征值和特征向量
- 三阶矩阵的特征值一般求解
- 在VisualStudio2005中调试客户端javascript
- 对称矩阵算法2
- sql2005的增删改和函数
- js中会用jstl
- CSS Sprite 网页优化
- PopupWindow使用
- DWZ-JUI中碰到的问题解决方法详解
- 在VisualStudio2005中求解实对称矩阵特征值算法包调查
- Android的网络应用 - 简单的C/S聊天室
- vim编辑器的使用
- UML几种图的绘制
- strtol 函数
- HDU 1395 2^x mod n = 1
- jquery 页面静态化 asp.net
- 【解题报告】 HDU 1272 小希的迷宫 并查集 判连通+判环
- 如何判断SD卡的剩余空间小于某个值