From mboxrd@z Thu Jan  1 00:00:00 1970
Return-Path: <gsl-discuss-return-5065-listarch-gsl-discuss=sources.redhat.com@sourceware.org>
Received: (qmail 1091 invoked by alias); 21 Mar 2007 10:10:25 -0000
Received: (qmail 1079 invoked by uid 22791); 21 Mar 2007 10:10:23 -0000
X-Spam-Check-By: sourceware.org
Received: from network-theory.com (HELO mail.network-theory.co.uk) (66.199.228.187)
    by sourceware.org (qpsmtpd/0.31) with ESMTP; Wed, 21 Mar 2007 10:10:19 +0000
Date: Wed, 21 Mar 2007 10:10:00 -0000
Message-ID: <87ejni7vqw.wl%bjg@network-theory.co.uk>
From: Brian Gough <bjg@network-theory.co.uk>
To: Patrick Alken <patrick.alken@colorado.edu>
Cc: gsl-discuss@sourceware.org
Subject: Re: generalized eigensystems
In-Reply-To: <20070313050349.GA22421@hippogriff.physics.drexel.edu>
References: <20070313050349.GA22421@hippogriff.physics.drexel.edu>
User-Agent: Wanderlust/2.14.0 (Africa) Emacs/21.3 Mule/5.0 (SAKAKI)
MIME-Version: 1.0 (generated by SEMI 1.14.6 - "Maruoka")
Content-Type: text/plain; charset=US-ASCII
X-Message-Mac: d4c8db55499c246067813c83b43bdde6
Mailing-List: contact gsl-discuss-help@sourceware.org; run by ezmlm
Precedence: bulk
List-Id: <gsl-discuss.sourceware.org>
List-Subscribe: <mailto:gsl-discuss-subscribe@sourceware.org>
List-Archive: <http://sourceware.org/ml/gsl-discuss/>
List-Post: <mailto:gsl-discuss@sourceware.org>
List-Help: <mailto:gsl-discuss-help@sourceware.org>, <http://sourceware.org/ml/#faqs>
Sender: gsl-discuss-owner@sourceware.org
X-SW-Source: 2007-q1/txt/msg00061.txt.bz2

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