LAPACK(5)——矩阵广义特征值问题和QZ分解

广义特征值问题，即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  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。