Re: generalized eigensystems

[Brian Gough <bjg at network-theory dot co dot uk>]

At Mon, 12 Mar 2007 23:03:49 -0600,
Patrick Alken wrote:
> But you can send along any output from the test program to me and
> I'll have a look at it.

Here is an interesting test case, it aborts with a "shouldn't get
here" error in extended precision but works in double precision (via
GSL_IEEE_MODE).

Might be worth also looking at the previous case you had that went
away when changing the optimisation - did you try it in
double-precision?

-- 
Brian Gough

$ ./a.out
gsl: gen.c:373: ERROR: gen_schur_standardize2 failed on complex block
Default GSL error handler invoked.
Aborted
$ GSL_IEEE_MODE=double-precision ./a.out
GSL_IEEE_MODE="double-precision,trap-common"
solution:
-8.817336437946092624e-16  0.000000000000000000e+00  6.919269373264336664e+00
 1.343730449186862330e+01  0.000000000000000000e+00  1.349256613371829339e+01
 6.035613520034543528e-08  0.000000000000000000e+00  8.100863692583372355e+00
-9.140256551631864568e-08  0.000000000000000000e+00  1.226784413824448983e+01

#include <stdio.h>
#include <gsl/gsl_eigen.h>
#include <gsl/gsl_ieee_utils.h>

int main ()
{
  double AA[16]   = {   8,  -5,  -5,  -7,
                      -10, -10, -10, -10,
                      -10, -10, -10, -10,
                      -10, -10, -10, -10 };
  double BB[16]   = {  7, -10,   -9,   -1,
                       4,   3,   -7,    9,
                      -4,   1,   -9,    9,
                      -7,   5,   -5,   -1} ;

  gsl_matrix_view A = gsl_matrix_view_array(AA, 4, 4);
  gsl_matrix_view B = gsl_matrix_view_array(BB, 4, 4);

  gsl_vector_complex * alpha = gsl_vector_complex_alloc(4);
  gsl_vector * beta = gsl_vector_alloc(4);

  gsl_eigen_gen_workspace * w = gsl_eigen_gen_alloc(4);
  gsl_eigen_gen_params (0, 0, 0, w);

  gsl_ieee_env_setup();
  gsl_eigen_gen (&A.matrix, &B.matrix, alpha, beta, w);

  printf("solution:\n");

  {
    gsl_vector_view a_real = gsl_vector_complex_real (alpha);
    gsl_vector_view a_imag = gsl_vector_complex_imag (alpha);
    int k;

    for (k = 0; k < 4; k++) { 
      printf("% .18e % .18e % .18e\n", 
             gsl_vector_get(&a_real.vector,k), 
             gsl_vector_get(&a_imag.vector,k),
             gsl_vector_get(beta,k));
    }
  }
}