Hi Juan,
Yesterday I sent an email to bug-gsl and but I couldn't check that they had
received it
because I only could see messages until August (
http://lists.gnu.org/archive/html/bug-gsl/).
The subject of the email was "Possible bug in Hessenberg-Triangular
Reduction". Is there
any other way to check it?
Regards,
MiguelGT
El 3 de septiembre de 2011 09:36, Juan Pablo Amorocho D. <jamor...@gmail.com
> escribió:
> 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 <
> mgarcia...@gmail.com>:
>
> > 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. <
> jamor...@gmail.com> 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 ]
> > */
> >
> >
>
_______________________________________________
Help-gsl mailing list
Help-gsl@gnu.org
https://lists.gnu.org/mailman/listinfo/help-gsl
