Incorrect eigenvectors with dggev

[abdullahhyder]

Hi. I sent this message in the LAPACK Google group. I am repeating my question here because I am unsure if the Google group is still in use. Thank you.

I have been trying to use dggev to get the eigenvectors of a generalized eigenvalue problem. The eigenvalues returned are correct, but with everything I have tried, the eigenvectors returned are incorrect, and the transpose of the matrix of eigenvectors is also incorrect.

Below is a simple example of how I am calling/using the dggev function in C++:

#include <iostream>
#include <cmath>
// #include <Accelerate/Accelerate.h>
using namespace std;

extern "C" {
    void dggev_(char *jobvl, char *jobvr, int *n, double *a, int *lda,
    double *b, int *ldb, double *alphar, double *alphai, double *beta, 
    double *vl, int *ldvl, double *vr, int *ldvr, double *work, int *lwork,
    int *info);
}

int main() {
    int N = 3;

    char jobvl = 'N'; // do not compute the left generalized eigenvectors
    char jobvr = 'V'; // compute the right generalized eigenvectors
    int orderN = N;
    int lda = N;
    int ldb = N;
    double alphar[N];
    double alphai[N];
    double beta[N];
    double vl[N*N];
    int ldvl = N;
    double vr[N*N];
    int ldvr = N;
    double work[N];
    int lwork = N*16;
    int info = 0;

    double* A_mat = new double[N*N];
    double* B_mat = new double[N*N];

    // Set to identity matrix
    for (int i = 0; i<N; i++) {
        for (int j = 0; j<N; j++) {
            if (i == j) {
                A_mat[i*N+j] = 1;
                B_mat[i*N+j] = 1;
            } else {
                A_mat[i*N+j] = 0;
                B_mat[i*N+j] = 0;
            }
        }
    }

    // initializing alphar, alphai, beta, etc.
    for (int i=0; i<N; i++) {
        alphar[i] = -1;
        alphai[i] = -1;
        beta[i] = -1;
        work[i] = -1;

        for (int j=0; j<N; j++) {
            vl[i*N+j] = -1;
            vr[i*N+j] = -1;
        }
    }

    // Print matrix A
    std::cout << "Matrix A: " << endl;
    for (int i = 0; i<N; i++) {
        for (int j = 0; j<N; j++) {
            std::cout << A_mat[i*N+j] << " ";
        }
        std::cout << endl;
    }

    // Print matrix B
    std::cout << "Matrix B: " << endl;
    for (int i = 0; i<N; i++) {
        for (int j = 0; j<N; j++) {
            std::cout << A_mat[i*N+j] << " ";
        }
        std::cout << endl;
    }

    dggev_(&jobvl, &jobvr, &orderN, (double *) A_mat, &lda,
        (double *) B_mat, &ldb, (double *) alphar, (double *) alphai,
        (double *) beta, (double *) vl, &ldvl, (double *) vr, &ldvr,
        (double *) work, &lwork, &info);

    // Print eigenvalues
    std::cout << "Eigenvalues: " << endl;
    for (int i = 0; i<N; i++) {
        if (beta[i]==0) {
            std::cout << "NaN ";
        } else {
            std::cout << alphar[i]/beta[i] << " ";
        }
    }
    std::cout << endl;

    // Print right eigenvector matrix
    std::cout << "Eigenvectors: " << endl;
    for (int i = 0; i<N; i++) {
        for (int j = 0; j<N; j++) {
            std::cout << vr[i*N+j] << " ";
        }
        std::cout << endl;
    }

}

The output from this code is:

>> g++ -c example.cpp -llapack -lblas -lgfortran -lm
>> g++ -o a.out example.o -llapack -lblas -lgfortran -lm
>> ./a.out 
Matrix A: 
1 0 0 
0 1 0 
0 0 1 
Matrix B: 
1 0 0 
0 1 0 
0 0 1 
Eigenvalues: 
1 1 1 
Eigenvectors: 
0 1 0 
0 0 0 
1 0 0

As you can see, the returned eigenvectors of these identity matrices include a vector of zeros instead of the <0, 1, 0> vector.

I have tried using "#include <Accelerate/Accelerate.h>" instead of "extern" to import the dggev method, but the result is the same. Any help would be greatly appreciated!

ScaLAPACK version: 2.1.0 (installed using macports) g++ version: g++ (MacPorts gcc12 12.3.0_0+stdlib_flag) 12.3.0

[weslleyspereira]

I quickly checked the code and

    double work[N];
    int lwork = N*16;

looks wrong to me. Did you check info in the output of dggev_()? You call dggev_() with lwork=-1 so that you have the optimal workspace size, or even use LAPACKE and then you don't need to worry about the workspace.

[weslleyspereira]

I have just verified and the code works if we replace double work[N]; by double work[N*16];. It is advisable to perform a workspace query to obtain the workspace size instead of setting it like you did.

Please, reopen if this actually does not solve your problem.