From: Brian Gough <bjg@network-theory.co.uk>
To: Patrick Alken <patrick.alken@colorado.edu>
Cc: gsl-discuss@sourceware.org
Subject: Re: generalized eigensystems
Date: Wed, 21 Mar 2007 10:10:00 -0000 [thread overview]
Message-ID: <87ejni7vqw.wl%bjg@network-theory.co.uk> (raw)
In-Reply-To: <20070313050349.GA22421@hippogriff.physics.drexel.edu>
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));
}
}
}
next prev parent reply other threads:[~2007-03-21 10:10 UTC|newest]
Thread overview: 15+ messages / expand[flat|nested] mbox.gz Atom feed top
2007-03-13 5:03 Patrick Alken
2007-03-16 15:49 ` Brian Gough
2007-03-20 19:17 ` Patrick Alken
2007-03-20 22:35 ` Leopold Palomo Avellaneda
2007-03-21 0:01 ` Patrick Alken
2007-03-21 7:32 ` Leopold Palomo Avellaneda
2007-03-21 14:47 ` Patrick Alken
2007-03-21 15:02 ` Leopold Palomo-Avellaneda
2007-03-21 10:10 ` Brian Gough [this message]
2007-03-21 18:02 ` Patrick Alken
2007-03-21 19:38 ` Brian Gough
2007-03-21 19:18 ` Patrick Alken
2007-03-30 21:34 ` Patrick Alken
2007-04-03 14:32 ` Brian Gough
2007-04-23 22:15 ` Patrick Alken
Reply instructions:
You may reply publicly to this message via plain-text email
using any one of the following methods:
* Save the following mbox file, import it into your mail client,
and reply-to-all from there: mbox
Avoid top-posting and favor interleaved quoting:
https://en.wikipedia.org/wiki/Posting_style#Interleaved_style
* Reply using the --to, --cc, and --in-reply-to
switches of git-send-email(1):
git send-email \
--in-reply-to=87ejni7vqw.wl%bjg@network-theory.co.uk \
--to=bjg@network-theory.co.uk \
--cc=gsl-discuss@sourceware.org \
--cc=patrick.alken@colorado.edu \
/path/to/YOUR_REPLY
https://kernel.org/pub/software/scm/git/docs/git-send-email.html
* If your mail client supports setting the In-Reply-To header
via mailto: links, try the mailto: link
Be sure your reply has a Subject: header at the top and a blank line
before the message body.
This is a public inbox, see mirroring instructions
for how to clone and mirror all data and code used for this inbox;
as well as URLs for read-only IMAP folder(s) and NNTP newsgroup(s).