help-gsl
[Top][All Lists]

## Re: [Help-gsl] About create_givens and precision

 From: Juan Pablo Amorocho D. Subject: Re: [Help-gsl] About create_givens and precision Date: Fri, 2 Sep 2011 20:55:16 +0200

```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",
}
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 ]
*/
create_givens(gsl_matrix_get(A, i - 1, j),
gsl_matrix_get(A, i, j),
&cs,
&sn);
/* invert so drot() works correctly (G -> G^t) */
sn = -sn;
/* compute G^t A(i-1:i, j:n) */
xv = gsl_matrix_subrow(A, i - 1, j, N - j);
yv = gsl_matrix_subrow(A, i, j, N - j);
gsl_blas_drot(&xv.vector, &yv.vector, cs, sn);

/* compute G^t B(i-1:i, i-1:n) */
xv = gsl_matrix_subrow(B, i - 1, i - 1, N - i + 1);
yv = gsl_matrix_subrow(B, i, i - 1, N - i + 1);
gsl_blas_drot(&xv.vector, &yv.vector, cs, sn);

if (U) {
/* accumulate U: U -> U G */
xv = gsl_matrix_column(U, i - 1);
yv = gsl_matrix_column(U, i);
gsl_blas_drot(&xv.vector, &yv.vector, cs, sn);
}

/* step 2: rotate columns i, i - 1 to kill B(i, i - 1) */

create_givens(-gsl_matrix_get(B, i, i),
gsl_matrix_get(B, i, i - 1),
&cs,
&sn);

/* invert so drot() works correctly (G -> G^t) */
sn = -sn;
/* compute B(1:i, i-1:i) G */
xv = gsl_matrix_subcolumn(B, i - 1, 0, i + 1);
yv = gsl_matrix_subcolumn(B, i, 0, i + 1);
gsl_blas_drot(&xv.vector, &yv.vector, cs, sn);

/* apply to A(1:n, i-1:i) */
xv = gsl_matrix_column(A, i - 1);
yv = gsl_matrix_column(A, i);
gsl_blas_drot(&xv.vector, &yv.vector, cs, sn);

if (V)
{
/* accumulate V: V -> V G */
xv = gsl_matrix_column(V, i - 1);
yv = gsl_matrix_column(V, i);
gsl_blas_drot(&xv.vector, &yv.vector, cs, sn);
}
}
}
}
return GSL_SUCCESS;
}

void create_givens (const double a, const double b, double *c, double *s) {
if (b == 0) {
*c = 1;
*s = 0;
} else if (fabs (b) > fabs (a)) {
double t = -a / b;
double s1 = 1.0 / sqrt (1 + t * t);
*s = s1;
*c = s1 * t;
} else {
double t = -b / a;
double c1 = 1.0 / sqrt (1 + t * t);
*c = c1;
*s = c1 * t;
}
}

void apply_threshold(gsl_matrix *A, double eps) {
int i, j;
for (i = 0; i < A->size1; i++) {
for (j = 0; j < A->size2; j++) {
if (fabs(gsl_matrix_get(A, i, j)) <= eps) {
gsl_matrix_set(A, i, j, 0.);
}
}
}
}

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");
}
}

void print_vector(gsl_vector *V) {
int i;

for (i = 0; i < V->size; i++) {
printf("\t%.8f", gsl_vector_get(V, i));
}
printf("\n");
}

void print_vector_complex(gsl_vector_complex *V) {
int i;

for (i = 0; i < V->size; i++) {
gsl_complex z = gsl_vector_complex_get (V, i);
printf("\t(%.8f,%.8fi)", GSL_REAL(z), GSL_IMAG(z));
}
printf("\n");
}

2011/9/2 Miguel García Torres <address@hidden>

> Dear Juan Pablo,
>
> Thanks you for the quick reply. I wasn't sure whether the floating point
> arithmetic was or not the problem. How can I know
> the number of decimals of the machine precision? I have done a small test.
> I have set to 0 all values that are lower or
> equal to 10-8. Then I get the same results. In my case depending on machine
> precision I get very different values
> in the output. Is there any way in GSL to avoid this situation?
>
> Thanks you!
>
> MiguelGT
>
> PS. For comparison purpose, the code is the following (maybe my code is not
> correct) :
>
> #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_matrix *);
> 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;
>   gsl_matrix *A = gsl_matrix_alloc(n, n);
>   gsl_matrix *B = gsl_matrix_alloc(n, n);
>
>   gsl_matrix_set(A, 0, 0, 1.);
>   gsl_matrix_set(A, 0, 1, 2.);
>   gsl_matrix_set(A, 0, 2, 1.);
>   gsl_matrix_set(A, 0, 3, 2.);
>   gsl_matrix_set(A, 1, 0, 3.);
>   gsl_matrix_set(A, 1, 1, 4.);
>   gsl_matrix_set(A, 1, 2, 3.);
>   gsl_matrix_set(A, 1, 3, 4.);
>   gsl_matrix_set(A, 2, 0, 1.);
>   gsl_matrix_set(A, 2, 1, 2.);
>   gsl_matrix_set(A, 2, 2, 1.);
>   gsl_matrix_set(A, 2, 3, 2.);
>   gsl_matrix_set(A, 3, 0, 3.);
>   gsl_matrix_set(A, 3, 1, 4.);
>   gsl_matrix_set(A, 3, 2, 3.);
>   gsl_matrix_set(A, 3, 3, 4.);
>
>   gsl_matrix_set(B, 0, 0, 5.);
>   gsl_matrix_set(B, 0, 1, 6.);
>   gsl_matrix_set(B, 0, 2, 5.);
>   gsl_matrix_set(B, 0, 3, 6.);
>   gsl_matrix_set(B, 1, 0, 7.);
>   gsl_matrix_set(B, 1, 1, 8.);
>   gsl_matrix_set(B, 1, 2, 7.);
>   gsl_matrix_set(B, 1, 3, 8.);
>   gsl_matrix_set(B, 2, 0, 5.);
>   gsl_matrix_set(B, 2, 1, 6.);
>   gsl_matrix_set(B, 2, 2, 5.);
>   gsl_matrix_set(B, 2, 3, 6.);
>   gsl_matrix_set(B, 3, 0, 7.);
>   gsl_matrix_set(B, 3, 1, 8.);
>   gsl_matrix_set(B, 3, 2, 7.);
>   gsl_matrix_set(B, 3, 3, 8.);
>
>   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",
>     }
>   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 ]
>                */
>               create_givens(gsl_matrix_get(A, i - 1, j),
>                             gsl_matrix_get(A, i, j),
>                             &cs,
>                             &sn);
>               /* invert so drot() works correctly (G -> G^t) */
>               sn = -sn;
>               /* compute G^t A(i-1:i, j:n) */
>               xv = gsl_matrix_subrow(A, i - 1, j, N - j);
>               yv = gsl_matrix_subrow(A, i, j, N - j);
>               gsl_blas_drot(&xv.vector, &yv.vector, cs, sn);
>
>               /* compute G^t B(i-1:i, i-1:n) */
>               xv = gsl_matrix_subrow(B, i - 1, i - 1, N - i + 1);
>               yv = gsl_matrix_subrow(B, i, i - 1, N - i + 1);
>               gsl_blas_drot(&xv.vector, &yv.vector, cs, sn);
>
>               if (U) {
>                 /* accumulate U: U -> U G */
>                 xv = gsl_matrix_column(U, i - 1);
>                 yv = gsl_matrix_column(U, i);
>                 gsl_blas_drot(&xv.vector, &yv.vector, cs, sn);
>               }
>
>               /* step 2: rotate columns i, i - 1 to kill B(i, i - 1) */
>
>               create_givens(-gsl_matrix_get(B, i, i),
>                             gsl_matrix_get(B, i, i - 1),
>                             &cs,
>                             &sn);
>
>               /* invert so drot() works correctly (G -> G^t) */
>               sn = -sn;
>               /* compute B(1:i, i-1:i) G */
>               xv = gsl_matrix_subcolumn(B, i - 1, 0, i + 1);
>               yv = gsl_matrix_subcolumn(B, i, 0, i + 1);
>               gsl_blas_drot(&xv.vector, &yv.vector, cs, sn);
>
>               /* apply to A(1:n, i-1:i) */
>               xv = gsl_matrix_column(A, i - 1);
>               yv = gsl_matrix_column(A, i);
>               gsl_blas_drot(&xv.vector, &yv.vector, cs, sn);
>
>               if (V)
>                 {
>                   /* accumulate V: V -> V G */
>                   xv = gsl_matrix_column(V, i - 1);
>                   yv = gsl_matrix_column(V, i);
>                   gsl_blas_drot(&xv.vector, &yv.vector, cs, sn);
>                 }
>             }
>         }
>     }
>   return GSL_SUCCESS;
> }
>
> void create_givens (const double a, const double b, double *c, double *s) {
>   if (b == 0) {
>     *c = 1;
>     *s = 0;
>   } else if (fabs (b) > fabs (a)) {
>     double t = -a / b;
>     double s1 = 1.0 / sqrt (1 + t * t);
>     *s = s1;
>     *c = s1 * t;
>   } else {
>     double t = -b / a;
>     double c1 = 1.0 / sqrt (1 + t * t);
>     *c = c1;
>     *s = c1 * t;
>   }
> }
>
> void apply_threshold(gsl_matrix *A, double eps) {
>   int i, j;
>   for (i = 0; i < A->size1; i++) {
>     for (j = 0; j < A->size2; j++) {
>       if (fabs(gsl_matrix_get(A, i, j)) <= eps) {
>         gsl_matrix_set(A, i, j, 0.);
>       }
>     }
>   }
> }
>
> 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%.20f", gsl_matrix_get(A, i, j));
>     }
>     printf("\n");
>   }
> }
>
> void print_vector(gsl_vector *V) {
>   int i;
>
>   for (i = 0; i < V->size; i++) {
>     printf("\t%.20f", gsl_vector_get(V, i));
>   }
>   printf("\n");
> }
>
> void print_vector_complex(gsl_vector_complex *V) {
>   int i;
>
>   for (i = 0; i < V->size; i++) {
>     gsl_complex z = gsl_vector_complex_get (V, i);
>     printf("\t(%.20f,%.20fi)", GSL_REAL(z), GSL_IMAG(z));
>   }
>   printf("\n");
> }
>
>
>
```

miguel.c
Description: Text Data

reply via email to