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

[Miguel García Torres]

I would like to say thanks to all of you to have a look at this issue!

Regards,

MiguelGT


2011/9/23 Patrick Alken <address@hidden>

> I just had a chance to look at this and agree with Brian's earlier
> response. Your test program satisfies the relations
>
>
> A = U H V^T
> B = U R V^T
>
> so the routine is working correctly for that case. This decomposition may
> not be unique which is why you may see a different result in that book you
> referenced.
>
> Patrick
>
>
> On 09/02/2011 10:19 PM, Miguel García Torres wrote:
>
>> 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<https://lists.gnu.org/mailman/listinfo/bug-gsl>
>>
>
>