[Lapack] dggev fails on non-singular system

[Patrick Alken]

Hello,

  I'd like to report a failure in the generalized nonsymmetric
eigensolver for the system (using lapack-3.1.1):

A =

   1   1   0
   0   0  -1
   1   0   0

B =

  -1   0  -1
   0  -1   0
   0   0  -1

This system has well defined generalized eigenvalues (note
in particular that B is invertible) - maple reports the
eigenvalues as:

-0.5000000007 + 0.8660254030 I, -0.5000000006 - 0.8660254035 I,
0.9999999994

but the DGGEV (DHGEQZ) driver fails to find any eigenvalues.

I have attached a simple test program to illustrate the
problem. Please forgive that it is written in C and
utilizes the GSL library to set up the matrices, as I normally
call lapack routines from C. Here is an illustration of how
to get it working:

> gcc -g -Wall -o testdggev testdggev.c -lgsl -llapack -lf77blas -lcblas 
> -latlas -lg2c
> ./testdggev 
info = 3

Since the system is 3x3, the return value of 'info' indicates no
eigenvalues converged and no eigenvectors were computed. After
stepping through the QZ iterations a bit, it seems like the QZ
iterations just don't split off eigenvalue blocks for this system,
and so possibly the calculation of the exceptional shifts needs
to modified?

Thanks,
Patrick Alken
-------------- next part --------------
#include <stdio.h>
#include <gsl/gsl_matrix.h>

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

void
lapack_dggev(gsl_matrix *A, gsl_matrix *B)
{
  char jobvl, jobvr;
  int N = (int) A->size1;
  double *alphar, *alphai, *beta;
  gsl_matrix *vr = gsl_matrix_alloc(N, N);
  double *work;
  int lwork, info;

  jobvl = 'N';
  jobvr = 'V';
  lwork = 5000;

  alphar = malloc(sizeof(double) * N);
  alphai = malloc(sizeof(double) * N);
  beta = malloc(sizeof(double) * N);
  work = malloc(sizeof(double) * lwork);

  dggev_(&jobvl,
         &jobvr,
         &N,
         A->data,
         (int *) &A->tda,
         B->data,
         (int *) &B->tda,
         alphar,
         alphai,
         beta,
         (double *) 0,
         (int *) &A->tda,
         vr->data, 
         (int *) &vr->tda,
         work,
         &lwork,
         &info);

  printf("info = %d\n", info);
}

int
main()
{
  gsl_matrix *A, *B;
  size_t N;
  double dat_A[] = { 1.0, 1.0, 0.0,
                     0.0, 0.0, -1.0,
                     1.0, 0.0, 0.0 };
  double dat_B[] = { -1.0, 0.0, -1.0,
                     0.0, -1.0, 0.0,
                     0.0, 0.0, -1.0 };
  gsl_matrix_view av, bv;

  N = 3;

  av = gsl_matrix_view_array(dat_A, N, N);
  bv = gsl_matrix_view_array(dat_B, N, N);

  A = gsl_matrix_alloc(N, N);
  B = gsl_matrix_alloc(N, N);

  gsl_matrix_transpose_memcpy(A, &av.matrix);
  gsl_matrix_transpose_memcpy(B, &bv.matrix);

  lapack_dggev(A, B);

  return 0;
}