LAPACK(5)——矩阵广义特征值问题和QZ分解
2011-07-16 19:12
411 查看
广义特征值问题,即Ax=
Bx,
在Matlab中,使用eig()求解一般特征值问题和广义特征值。[V,D] = eig(A,B,flag), A和B时方阵,flag用来选择算法,'qz'表示选择使用QZ算法。
也可以直接调用qz()来求解,[AA,BB,Q,Z,V] = qz(A,B,flag), flag 表示使用复数或实数计算,默认取值为复数。
在Lapack中,有四个函数都是用来求解广义特征值的,
区别在于前两个分解之后会输出舒尔形式,后两个则输出广义特征向量。而且gegs和gegv都被gges和ggev代替。两个都会用QZ分解求解广义特征值。
LAPACK也给出了QZ分解的函数dhgeqz,但要求输入H,T矩阵,对于一般的方阵,可以使用dgghrd将输入的方阵A,B变换成H,T矩阵。
下面给出这四个函数的原型和测试程序。
结果如下,
代码中调用dhgeqz的时候,没有使用迭代,暂时还没有弄清楚在dhgeqz内部有没有实现迭代,测试中结果是对的。需要注意的是,Matlab的qz函数会给出S,T矩阵,但Lapack的dgges给出的S-T结果并不一致,原因还没有弄明白的。
关于dgges这个函数的一个参数LAPACK_D_SELECT3 selctg,有个参考,http://software.intel.com/sites/products/documentation/hpc/mkl/mklman/lle/lle_cinterface.htm
这里用到了函数指针,http://www.upsdn.net/html/2004-11/40.html。
Bx,
在Matlab中,使用eig()求解一般特征值问题和广义特征值。[V,D] = eig(A,B,flag), A和B时方阵,flag用来选择算法,'qz'表示选择使用QZ算法。
也可以直接调用qz()来求解,[AA,BB,Q,Z,V] = qz(A,B,flag), flag 表示使用复数或实数计算,默认取值为复数。
在Lapack中,有四个函数都是用来求解广义特征值的,
?GEGS Computes the generalized eigenvalues, Schur form, and left and/or right Schur vectors for a pair of non-symmetric matrices. ?GGES Computes the generalized eigenvalues, Schur form, and left and/or right Schur vectors for a pair of non-symmetric matrices. ?GEGV Computes the generalized eigenvalues, and left and/or right generalized eigenvectors for a pair of non-symmetric matrices. ?GGEV Computes the generalized eigenvalues, and left and/or right generalized eigenvectors for a pair of non-symmetric matrices.
区别在于前两个分解之后会输出舒尔形式,后两个则输出广义特征向量。而且gegs和gegv都被gges和ggev代替。两个都会用QZ分解求解广义特征值。
LAPACK也给出了QZ分解的函数dhgeqz,但要求输入H,T矩阵,对于一般的方阵,可以使用dgghrd将输入的方阵A,B变换成H,T矩阵。
下面给出这四个函数的原型和测试程序。
#include <iostream> #include <iomanip> #include <cmath> #include <complex> using namespace std; typedef complex<double> dcomplex_t; //lapacke headers #include "lapacke.h" #include "lapacke_config.h" #include "lapacke_utils.h" extern "C" { lapack_int LAPACKE_dggev( int matrix_order, char jobvl, char jobvr, lapack_int n, double* a, lapack_int lda, double* b, lapack_int ldb, double* alphar, double* alphai, double* beta, double* vl, lapack_int ldvl, double* vr, lapack_int ldvr ); lapack_int LAPACKE_dgges( int matrix_order, char jobvsl, char jobvsr, char sort, LAPACK_D_SELECT3 selctg, lapack_int n, double* a, lapack_int lda, double* b, lapack_int ldb, lapack_int* sdim, double* alphar, double* alphai, double* beta, double* vsl, lapack_int ldvsl, double* vsr, lapack_int ldvsr ); lapack_logical selectg(const double* AR,const double* AI,const double* B){ if(fabs(*AI)<1e-8) return 0; else return 1; } lapack_int LAPACKE_dgghrd( int matrix_order, char compq, char compz, lapack_int n, lapack_int ilo, lapack_int ihi, double* a, lapack_int lda, double* b, lapack_int ldb, double* q, lapack_int ldq, double* z, lapack_int ldz ); lapack_int LAPACKE_dhgeqz( int matrix_order, char job, char compq, char compz, lapack_int n, lapack_int ilo, lapack_int ihi, double* h, lapack_int ldh, double* t, lapack_int ldt, double* alphar, double* alphai, double* beta, double* q, lapack_int ldq, double* z, lapack_int ldz ); } void PrintM(int M,int N,double* A){ int i = 0, j = 0; for(i=0;i<M;i++){ for(j=0;j<N;j++) cout << setw(10) << A[i*N+j] << " "; cout << endl; } } int main(){ int N = 4; double A[16] = {3.9, 12.5, -34.5, -0.5, 4.3, 21.5, -47.5, 7.5, 4.3, 21.5, -43.5, 3.5, 4.4, 26.0, -46.0, 6.0}; double B[16] = {1.0, 2.0, -3.0, 1.0, 1.0, 3.0, -5.0, 4.0, 1.0, 3.0, -4.0, 3.0, 1.0, 3.0, -4.0, 4.0}; double alphar[4]; double alphai[4]; double beta[4]; int sdim[1]; int info = LAPACKE_dgges(LAPACK_ROW_MAJOR,'N', 'N', 'S', selectg, N, A, N, B, N, sdim,alphar,alphai,beta, NULL,N,NULL,N); cout << "dgges result:" << endl; cout << "info = " << info << endl; cout << "sdim = " << sdim[0] << endl; cout << "alpha:" << endl; PrintM(1,4,alphar); PrintM(1,4,alphai); cout << "beta:" << endl; PrintM(1,4,beta); cout << "eigenvalue:" << endl; for(int i=0;i<4;i++) cout << dcomplex_t(alphar[i]/beta[i],alphai[i]/beta[i]) << " "; cout << endl; int k = 0; for(k=0;k<N;k++){ alphar[k] = 0; alphai[k] = 0; beta[k] = 0; } // dggev info = LAPACKE_dggev(LAPACK_ROW_MAJOR,'N','N', N,A,N,B,N, alphar,alphai,beta, NULL,N,NULL,N); cout << "dggev result:" << endl; cout << "info = " << info << endl; cout << "alpha:" << endl; PrintM(1,4,alphar); PrintM(1,4,alphai); cout << "beta:" << endl; PrintM(1,4,beta); cout << "eigenvalue:" << endl; for(int i=0;i<4;i++) cout << dcomplex_t(alphar[i]/beta[i],alphai[i]/beta[i]) << " "; cout << endl; // using QZ iteration method to get eigenvalue cout << "QZ method" << endl; int NA = N; info = LAPACKE_dgghrd(LAPACK_ROW_MAJOR,'N','N', NA,1,NA, A,NA,B,NA,NULL,NA,NULL,NA); if(info!=0){ cout << "error when reduce A/B to H/T" << endl; exit(-1); } for(k=0;k<N;k++){ alphar[k] = 0; alphai[k] = 0; beta[k] = 0; } info = LAPACKE_dhgeqz(LAPACK_ROW_MAJOR,'E','N','N', NA,1,NA,A,NA,B,NA, alphar,alphai,beta, NULL,NA,NULL,NA); if(info!=0){ } cout << "alpha:" << endl; PrintM(1,4,alphar); PrintM(1,4,alphai); cout << "beta:" << endl; PrintM(1,4,beta); cout << "eigenvalue:" << endl; for(int i=0;i<4;i++) cout << dcomplex_t(alphar[i]/beta[i],alphai[i]/beta[i]) << " "; cout << endl; return 0; }
结果如下,
dgges result: info = 0 sdim = 2 alpha: 0.857143 0.857143 0.617213 4 1.14286 -1.14286 0 0 beta: 0.285714 0.285714 0.308607 1 eigenvalue: (3,4) (3,-4) (2,0) (4,0) dggev result: info = 0 alpha: 8.23538 3.54123 0.617213 4 10.9805 -4.72164 0 0 beta: 2.74513 1.18041 0.308607 1 eigenvalue: (3,4) (3,-4) (2,0) (4,0) QZ method alpha: 8.23538 3.54123 0.617213 4 10.9805 -4.72164 0 0 beta: 2.74513 1.18041 0.308607 1 eigenvalue: (3,4) (3,-4) (2,0) (4,0)
代码中调用dhgeqz的时候,没有使用迭代,暂时还没有弄清楚在dhgeqz内部有没有实现迭代,测试中结果是对的。需要注意的是,Matlab的qz函数会给出S,T矩阵,但Lapack的dgges给出的S-T结果并不一致,原因还没有弄明白的。
关于dgges这个函数的一个参数LAPACK_D_SELECT3 selctg,有个参考,http://software.intel.com/sites/products/documentation/hpc/mkl/mklman/lle/lle_cinterface.htm
这里用到了函数指针,http://www.upsdn.net/html/2004-11/40.html。
相关文章推荐
- 大规模稀疏矩阵的广义特征值问题——CLAPACK, IETL库尝试
- 【机器学习】K-Means 聚类是特殊的矩阵分解问题
- 矩阵分解 (特征值/奇异值分解+SVD+解齐次/非齐次线性方程组)
- 矩阵的特征值 => 矩阵的幂 => 广义斐波拉契数列的通项公式
- 矩阵的特征值分解
- 矩阵的特征值 => 矩阵的幂 => 广义斐波拉契数列的通项公式
- 矩阵的特征值分解与奇异值分解的几何意义
- 矩阵的特征值分解与奇异值分解
- 矩阵的特征值分解和奇异值分解
- 矩阵的特征值 => 矩阵的幂 => 广义斐波拉契数列的通项公式
- 关于矩阵分解:特征值分解 svd分解 mf分解 lmf分解 pca 以及个性化推荐 fm ffm als
- (转)QR分解求矩阵的全部特征值
- 【转】矩阵的奇异值分解和特征值分解的异同
- 人工智能里的数学修炼 | 矩阵的花样分解:特征值分解(EVD)、相似对角化、QR分解、Schur分解、奇异值分解(SVD)的概念纠缠与详解
- QR分解求矩阵全部特征值
- 什么是特征向量,特征值,矩阵分解
- 矩阵的特征值分解与奇异值分解
- matlab矩阵及其基本运算—特征值分解和奇异值分解
- QR分解求矩阵全部特征值
- 矩阵的“特征值分解”和“奇异值分解”区别