Re: [Help-gsl] About create_givens and precision

[Juan Pablo Amorocho D.]

Hi Miguel,

That code really looked very gsl-liked! Have summited a bug report?

-- juan

Am Freitag, 2. September 2011 schrieb Miguel García Torres <
address@hidden>:
> Hi Juan,
> Thanks you. I have tested the code and it is a GSL bug. The function
 "linalg_hesstri_decomp"
> corresponds to GSL function "gsl_linalg_hesstri_decomp". I included it in
my file to debug
> purpose.
> Kins regards,
> MiguelGT
>
> El 2 de septiembre de 2011 20:55, Juan Pablo Amorocho D. <
address@hidden> escribió:
>
> Hi Miguel,
>
> A couple of things. I assume you are trying to do a Hessenberg-Triangular
Reduction. I looked it up in Matrix Computations(MC), 3rd Ed.   and it is
Alg. 7.7.1, page 380.  There is an example there that I ran (see below)
using your code and the rounding doesn't have any effect. In fact, your code
seems to have a bug. Matrices B, U, and V are correct according to the
example of MC which, unfortunately, only provides values up to the 4th
figure after the coma. Now the matrix A is the problem. The right A should
be
>
> [ -2.5849 1.5413 2.4221]
> [-9.7631 0.0874 1.9239 ]
> [0.0000 2.7233 -.7612 ]
>
> so I think you might have a bug in your code.
>
>
>
> #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(int);
> int linalg_hesstri_decomp(gsl_matrix * A, gsl_matrix * B, gsl_matrix * U,
gsl_matrix * V, gsl_vector * work, int do_round);
> void create_givens (const double a, const double b, double *c, double *s);
> void print_matrix(gsl_maatrix *);
> void print_vector(gsl_vector *);
> void apply_threshold(gsl_matrix *, double);
>
> int main (void) {
>   test_hesstri(0);
>   test_hesstri(1);
>
>   return 0;
> }
>
> void test_hesstri(int do_round) {
>   //int n = 4;
>   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);
>
>   linalg_hesstri_decomp(A, B, U, V, work, do_round);
>
>   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);
>   printf("Vector work:\n");
>   print_vector(work);
> }
>
>
> int linalg_hesstri_decomp(gsl_matrix * A, gsl_matrix * B, gsl_matrix * U,
gsl_matrix * V, gsl_vector * work, int do_round) {
>   const double eps = 1e-8;
>   const size_t N = A->size1;
>
>   if ((N != A->size2) || (N != B->size1) || (N != B->size2))
>     {
>       GSL_ERROR ("Hessenberg-triangular reduction requires square
matrices",
>                  GSL_ENOTSQR);
>     }
>   else if (N != work->size)
>     {
>       GSL_ERROR ("length of workspace must match matrix dimension",
>                  GSL_EBADLEN);
>     }
>   else
>     {
>       double cs, sn;          /* rotation parameters */
>       size_t i, j;            /* looping */
>       gsl_vector_view xv, yv; /* temporary views */
>
>       /* B -> Q^T B = R (upper triangular) */
>       gsl_linalg_QR_decomp(B, work);
>       if (do_round) {
>         apply_threshold(B, eps);
>       }
>       /* A -> Q^T A */
>       gsl_linalg_QR_QTmat(B, work, A);
>       if (do_round) {
>         apply_threshold(A, eps);
>       }
>       /* initialize U and V if desired */
>       if (U) {
>           gsl_linalg_QR_unpack(B, work, U, B);
>         }
>       else
>         {
>           /* zero out lower triangle of B */
>           for (j = 0; j < N - 1; ++j)
>             {
>               for (i = j + 1; i < N; ++i)
>                 gsl_matrix_set(B, i, j, 0.0);
>             }
>         }
>
>       if (V)
>         gsl_matrix_set_identity(V);
>
>       if (N < 3)
>         return GSL_SUCCESS; /* nothing more to do */
>
>       /* reduce A and B */
>       for (j = 0; j < N - 2; ++j) {
>           for (i = N - 1; i >= (j + 2); --i)
>             {
>               /* step 1: rotate rows i - 1, i to kill A(i,j) */
>
>               /*
>                * compute G = [ CS SN ] so that G^t [ A(i-1,j) ] = [ * ]
>                *             [-SN CS ]             [ A(i, j)  ]   [ 0 ]
>                */
>
>