Hi, the information about the possible bug is
[*] The version number of GSL: 1.15
[*] The hardware and operating system: Intel(R) Core(TM)2 Duo CPU P9600
@ 2.66GHz 8GB RAM, Linux (Ubuntu karmic)
[*] The compiler used, including version number and compilation options: gcc
4.4.3 options:-lgsl -lgslcblas
[*] A description of the bug behavior
I have tried to port the Hessenberg-Triangular Reduction from gsl to Java.
However I found some wrong outputs. I sent an email to gsl-help and told
me the it could be a bug in the code:
1) In GSL-help, Juan Pablo told me that Looking up the book Matrix
Computation, 3rd Ed., Alg. 7.7.1 and page 380, he found is an example
of the Hessenberg-Triangular Reduction. The values of matrices match
all but A. The right output A should be
[ -2.5849 1.5413 2.4221]
[-9.7631 0.0874 1.9239 ]
[0.0000 2.7233 -.7612 ]
The output I get with GSL is
[-2.58492131056598850591 -1.10096376512636062728
-2.20192753025272169864]
[9.76310308345568600430 -0.06246341106135822052 1.93636574290210594640]
[-0.00000000000000044409 -1.94524474299697525126 1.18406201747641981470]
2) When I tried ti run the Hessenberg-Triangular Reduction for the following
matrices:
A =
[1 2 1 2]
[3 4 3 4]
[1 2 1 2]
[3 4 3 4]
B =
[5 6 5 6]
[7 8 7 8]
[5 6 5 6]
[7 8 7 8]
I get different output depending on the machine precision. When it has
to calculate the
givens rotation (create_givens) in gsl_linalg_hesstri_decomp function,
values are
very small (~10-15) and so I get vry different outputs for different
machines. Should
this thing be controlled?
[*] A short program which exercises the bug
include<stdio.h>
#include<stdlib.h>
#include<gsl/gsl_math.h>
#include<gsl/gsl_vector.h>
#include<gsl/gsl_matrix.h>
#include<gsl/gsl_blas.h>
#include<gsl/gsl_eigen.h>
#include<gsl/gsl_linalg.h>
void test_hesstri();
void print_matrix(gsl_matrix *);
int main (void) {
test_hesstri();
return 0;
}
void test_hesstri() {
int n = 3;
gsl_matrix *A = gsl_matrix_alloc(n, n);
gsl_matrix *B = gsl_matrix_alloc(n, n);
gsl_matrix_set(A, 0, 0, 10);
gsl_matrix_set(A, 0, 1, 1);
gsl_matrix_set(A, 0, 2, 2);
gsl_matrix_set(A, 1, 0, 1);
gsl_matrix_set(A, 1, 1, 2);
gsl_matrix_set(A, 1, 2, -1);
gsl_matrix_set(A, 2, 0, 1);
gsl_matrix_set(A, 2, 1, 1);
gsl_matrix_set(A, 2, 2, 2);
gsl_matrix_set(B, 0, 0, 1);
gsl_matrix_set(B, 0, 1, 2);
gsl_matrix_set(B, 0, 2, 3);
gsl_matrix_set(B, 1, 0, 4);
gsl_matrix_set(B, 1, 1, 5);
gsl_matrix_set(B, 1, 2, 6);
gsl_matrix_set(B, 2, 0, 7);
gsl_matrix_set(B, 2, 1, 8);
gsl_matrix_set(B, 2, 2, 9);
gsl_matrix *U = gsl_matrix_alloc(n, n);
gsl_matrix *V = gsl_matrix_alloc(n, n);
gsl_vector *work = gsl_vector_alloc(n);
gsl_linalg_hesstri_decomp(A, B, U, V, work);
printf(":::::::::::::::::::::::::::::::::::::::\n");
printf("[D]Matriz A:\n");
print_matrix(A);
printf("[D]Matriz B:\n");
print_matrix(B);
printf("[D]Matriz U:\n");
print_matrix(U);
printf("[D]Matriz V:\n");
print_matrix(V);
}
void print_matrix(gsl_matrix *A) {
int i, j;
for (i = 0; i< A->size1; i++) {
for (j = 0; j< A->size2; j++) {
printf("\t%.8f", gsl_matrix_get(A, i, j));
}
printf("\n");
}
}
Thanks you,
MiguelGT
_______________________________________________
Bug-gsl mailing list
address@hidden
https://lists.gnu.org/mailman/listinfo/bug-gsl