[Bug-gsl] Possible bug in Hessenberg-Triangular Reduction.

[Miguel García Torres]

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