How do I get specified Eigenvectors from the gener

Question

I am grad student trying to rewrite my MATLAB prototype code into C++ code using Eigen and LAPACK. Generalised eigenvalue solver (A*x=lamba*B*x) takes some part in this program. Because Eigen's generalised eigen solver could result in wrong eigenvalue (in my experience), I decided to go with LAPACK (using Accelerate framework in OS X)

The code below is the C++ translation of tgevc example http://www.nag.com/numeric/fl/manual/pdf/F08/f08ykf.pdf that I wrote. Everything looks fine unless I specify the eigenvalues for which right eigenvectors are calculated.

when I set

char howmny = 'B';

It calculates every eigenvectors for every eigenvalues. And it works FINE. The values are correct also for larger problems.

However, when I set

char howmny = 'S';

The result just doesn't make sense. It looks like nothing.

The reason for doing this was for possible speed up in the first place. However, profiling shows that it has less than 4% of the entire process (for larger problems). So it won't change much. Anyway I am still wondering whether I understand something Wrong.

Here is my sample code.

#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;
}