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

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

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

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

char howmny = 'B';

Он рассчитывает каждый собственный вектор для каждого собственного значения. И это работает FINE. Значения верны и для более крупных проблем.

Тем не менее, когда я установил

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