Как получить заданные собственные векторы из обобщенной факторизации Шура пары матриц с использованием LAPACK?

Я аспирант, пытающийся переписать мой код прототипа MATLAB в код C++, используя Eigen и LAPACK. Обобщенный решатель собственных значений (A*x=lamba*B*x) принимает некоторое участие в этой программе. Поскольку обобщенный собственный решатель Eigen может привести к неправильному собственному значению (по моему опыту), я решил использовать LAPACK (используя платформу Accelerate в OS X)

Приведенный ниже код представляет собой перевод примера tgevc на C++ http://www.nag.com/numeric/fl/manual/pdf/F08/f08ykf.pdf, который я написал. Все выглядит нормально, если я не укажу собственные значения, для которых вычисляются правильные собственные векторы.

когда я установил

char howmny = 'B';

Он вычисляет все собственные векторы для каждого собственного значения. И это работает ОТЛИЧНО. Значения верны и для более крупных задач.

Однако, когда я установил

char howmny = 'S';

Результат просто не имеет смысла. Похоже ничего.

Причиной для этого было возможное ускорение в первую очередь. Однако профилирование показывает, что на него приходится менее 4% всего процесса (для более крупных проблем). Так что это не сильно изменится. Во всяком случае, мне все еще интересно, понимаю ли я что-то не так.

Вот мой пример кода.

#include <iostream>
#include <cmath>
#include <Eigen/Core>
#include <Accelerate/Accelerate.h>

using namespace std;
using namespace Eigen;

int main(int argc, const char * argv[])
{
    // parameters
    int nmax = 5, lda = nmax, ldb = nmax, ldq = nmax, ldz = nmax, lwork=ldq*6, N=5, M;

    // Local scalars
    int ihi, ilo, info, irows, icols, jwork;
    char compq, compz, job, job1, job2, side;

    // Local arrays
    MatrixXd A(lda,nmax), B(ldb, nmax), Q(ldq,ldq), Z(ldz, ldz);
    VectorXd alphai(nmax), alphar(nmax), beta(nmax), lscale(nmax), rscale(nmax), tau(nmax), work(lwork);

    cout << " F08XEF Example Program Results " << endl;

    for (int i = 0; i < N ; i++) {
        A(i,0) = (i+1)*1.0;
    }
    for (int j = 1; j < N ; j++) {
        A.col(j) = A.col(0).array()*A.col(j-1).array();
    }

    B.block(0,0,N,N) = A.block(0, 0, N, N).transpose();
    cout << "Matrix A is" << endl;
    cout << A.block(0, 0, N, N) << endl;
    cout << "Matrix B is" << endl;
    cout << B.block(0,0,N,N) << endl;

    // Balance matrix pair (A, B)
    job = 'B';

    dggbal_(&job, &N, A.data(), &lda, B.data(), &ldb, &ilo, &ihi, lscale.data(), rscale.data(), work.data(), &info);
    cout << "Matrix A after balacing" << endl;
    cout << A.block(0, 0, N, N) << endl;
    cout << "Matrix B after balacing" << endl;
    cout << B.block(0,0,N,N) << endl;

    // Reduce B to triangular form using QR
    irows = ihi + 1 - ilo;
    icols = N + 1 - ilo;
    dgeqrf_(&irows, &icols, B.block(ilo-1,ilo-1,irows,icols).data(), &ldb, tau.data(), work.data(), &lwork, &info);

    // Apply the orthogonal transformation t matrix A
    job1 = 'L';
    job2 = 'T';
    dormqr_(&job1, &job2, &irows, &icols, &irows, B.block(ilo-1,ilo-1,irows,icols).data(), &ldb, tau.data(), A.block(ilo-1,ilo-1,irows,irows).data(), &lda, work.data(), &lwork, &info);
    Z.block(0, 0, N, N) = MatrixXd::Identity(N, N);
    // Compute the generalized Hassenberg form of (A, B)
    compq = 'V';
    compz = 'V';
    dgghrd_(&compq,&compz,&N,&ilo,&ihi,A.block(ilo-1,ilo-1,irows,irows).data(),&lda, B.block(ilo-1,ilo-1,irows,irows).data(),&ldb,Q.data(),&ldq,Z.data(),&ldz,&info);
    cout << "Matrix A in Hessenber form" << endl;
    cout << A.block(0, 0, N, N) << endl;
    cout << "Matrix B is triangular" << endl;
    cout << B.block(0, 0, N, N) << endl;

    // Routine dhgeqz
    // Worspace query : jwork = -1
    jwork = -1;
    job = 'S';
    dhgeqz_(&job, &compq, &compz, &N, &ilo, &ihi, A.data(), &lda, B.data(), &ldb, alphar.data(), alphai.data(), beta.data(), Q.data(), &ldq, Z.data(), &ldz, work.data(), &jwork, &info);
    cout << "Minimal required lwork = " << round(work(0)) << endl;
    cout << "Actual value of lwork = " << lwork << endl;
    // Compute the generalized Schur form if the workspace lwork is adequate
    if (round(work(0)) <= lwork)
    {
        dhgeqz_(&job, &compq, &compz, &N, &ilo, &ihi, A.data(), &lda, B.data(), &ldb, alphar.data(), alphai.data(), beta.data(), Q.data(), &ldq, Z.data(), &ldz, work.data(), &lwork, &info);
    }

    // Print the generalized eigenvalue

    for (int i = 0; i < N ; i++) {
        if (beta(i) != 0.0)
        {
            cout << alphar(i)/beta(i) << "," << alphai(i)/beta(i)<< endl;
        }
    }


    // Compute right generalized eigenvectors of the balaced matrix

    // !This is where the problem shows up!
    char howmny = 'S';
    side = 'R';
    __CLPK_logical select[5]; // int selec[5] works same.

    select[0]=false;
    select[1]=false;
    select[2]=false;
    select[3]=false;
    select[4]=false;


    dtgevc_(&side, &howmny, select, &N, A.data(), &lda, B.data(), &ldb, Q.data(), &ldq, Z.data(), &ldz, &N, &M, work.data(), &info);

    job = 'B';
    dggbak_(&job, &side, &N, &ilo, &ihi, lscale.data(), rscale.data(), &N, Z.data(), &ldz, &info);

    cout << "Right eigenvectors" << endl;
    cout << Z << endl;

    return 0;
}