, 3 min read

Computing Generalized Eigenvalues with LAPACK

Original post is here eklausmeier.goip.de/blog/2025/05-20-computing-generalized-eigenvalues-with-lapack.


Task at hand: compute generalized eigenvalues for

$$ A x = \lambda B x, \qquad A, B \in \mathbb{C}^{n\times n}, \quad x\in\mathbb{C}^n, \quad\lambda\in\mathbb{C} $$

using the C programming language and the Fortran library LAPACK.

Previously I had used cqzhes() and cqzval() from EISPACK.

1. Installing LAPACK. Watch out for any pre-existing BLAS library, e.g., GotoBLAS.

pacman -S blas cblas lapack lapack-doc

Check man page for zgghrd, zhgeqz, and zggev:

man zgghrd
man zhgeqz
man zggev

Checking the so-file with objdump -d and noticing that the Fortran symbols end with _:

objdump -d /usr/lib/liblapack.so

Info on the library liblapack.so to which we link:

ldd /usr/lib/liblapack.so
        linux-vdso.so.1 (0x000077c056586000)
        libblas.so.3 => /usr/lib/libblas.so.3 (0x000077c056481000)
        libgfortran.so.5 => /usr/lib/libgfortran.so.5 (0x000077c055200000)
        libm.so.6 => /usr/lib/libm.so.6 (0x000077c055108000)
        libgcc_s.so.1 => /usr/lib/libgcc_s.so.1 (0x000077c056454000)
        libc.so.6 => /usr/lib/libc.so.6 (0x000077c054f18000)
        /usr/lib64/ld-linux-x86-64.so.2 (0x000077c056588000)

zggev calls zgghrd and zhgeqz.

# zggev ## xerbla ## zgeqrf ## zggbak ## zggbal ## zgghrd ### xerbla ### zlartg ### zlaset ### zrot ## zhgeqz ### xerbla ### zlartg ### zlaset ### zrot ### zscal ## zlacpy ## zlascl ## zlaset ## ztgevc ## zungqr ## zunmqr

The methods in question use the QZ algorithm.

2. Example program. Compile and link with

gcc -Wall qztest.c -o qztest -llapack

Below program uses the following rules when linking to Fortran on Linux:

  1. All arguments passed to Fortran are pointers.
  2. Arrays in Fortran store their values in column order, in contrast to row order in C. Therefore we have to transpose the matrix.
  3. Add underscore as suffix.

The program qztest.c is:

#include <stdio.h>
#include <string.h>
#include <math.h>
#include <complex.h>

extern void zggev_(char *jobvl, char *jobvr, int *n, complex *a, int *lda, complex *b,
    int *ldb, complex *alpha, complex *beta, complex *vl, int *ldvl, complex *vr, int *ldvr,
    complex *work, int *lwork, double *rwork, int *info);



#define N	5
double stoer[] = {	// Stoer/Burlisch (1970):``Introduction to Numerical Analysis", page 374
    12,  1,   0,   0,   0,
    1,   9,   1,   0,   0,
    0,   1,   6,   1,   0,
    0,   0,   1,   3,   1,
    0,   0,   0,   1,   0
};


int main (int argc, char *argv[]) {
    int i, j, lda=N, ldb=N, ldvl=N, ldvr=N, lwork=N*N, n=N, info=0;
    char jobvl[8], jobvr[8];
    complex a[N*N], b[N*N], alpha[N], beta[N], work[N*N], vl[N*N], vr[N*N], lambda;
    double rwork[8*N];

    memset(b,0,sizeof(*b)*N*N);	// set b[] to zero matrix
    for (i=0; i<N; ++i) {
        b[i+i*N] = 1;
        for (j=0; j<N; ++j)
            a[i+j*N] = stoer[i*N+j];	// transpose
    }

    jobvl[0] = 'N';	// do not compute the left generalized eigenvectors
    jobvr[0] = 'N'; // do not compute the right generalized eigenvectors
    zggev_(jobvl, jobvr, &n, a, &lda, b, &ldb, alpha, beta, vl, &ldvl, vr, &ldvr, work, &lwork, rwork, &info);
    printf("zggev_.info=%d\n",info);

    puts("lambda =");
    for (i=0; i<N; ++i) {
        lambda = alpha[i] / beta[i];
        printf("\t(%11.6f %11.6f)\n",creal(lambda),cimag(lambda));
    }

    return 0;
}

It prints

lambda =
        (  12.316876    0.000000)
        (   9.016136    0.000000)
        (   6.000000    0.000000)
        (   2.983864    0.000000)
        (  -0.316876    0.000000)

3. Another example.

$$ A = \pmatrix{ -9 & 0 & -3\cr -10 & 16 & -2\cr -1 & 28 & 1\cr }, \qquad B = \pmatrix{ -1 & 0 & 2\cr 1 & -2 & 3\cr 3 & -2 & 2\cr }. $$

This should print

        (   6.441386    0.000000)
        (  -1.220693    1.495264)
        (  -1.220693   -1.495264)

Crosscheck with Wolfram Alpha:

Eigenvalues[ {{-1,0,2},{1,-2,3},{3,-2,2}}^(-1) * {{-9,0,-3},{-10,16,-2},{-1,28,1}} ]

Added 22-Jul-2025: This post was mentioned in Machine Learning Engineer: 27th May 2025.