Re: Trouble with Cholesky

[Frank Boehme <f dot boehme at cs dot ucc dot ie>]

Frank Boehme wrote:

Ok, I made a mistake:

>   gsl_vector_set (v, 0, 3.0);
>   gsl_vector_set (v, 0, 3.0);
                      ^^^

But anyway, I believe there is a bug in the Cholesky algorithms. Here is
the correct demo. Householder works, Cholesky doesn't:

======cut=here=====

#include <stdio.h>
#include <gsl/gsl_linalg.h>
int main() {
  gsl_matrix * m = gsl_matrix_alloc (2, 2);
  gsl_vector * v = gsl_vector_alloc (2);
/* this matrix is symmetrical and positiv definit: */ 
  gsl_matrix_set (m, 0, 0, 3.0);
  gsl_matrix_set (m, 0, 1, 0.0);
  gsl_matrix_set (m, 1, 0, 0.0);
  gsl_matrix_set (m, 1, 1, 3.0);
  gsl_vector_set (v, 0, 3.0);
  gsl_vector_set (v, 1, 3.0);
/* uncomment one of the following two lines to see the difference */
  gsl_linalg_cholesky_svx (m, v);
/*  gsl_linalg_HH_svx (m, v); */
  printf("%g\n", gsl_vector_get (v, 0)); 
  printf("%g\n", gsl_vector_get (v, 1)); 
}

======end===========

-- 
Dr Frank Boehme                      | Email: f.boehme@cs.ucc.ie
National University of Ireland, Cork | phone: +353-21-4903163
Dept of Computer Science             |   fax: +353-21-4903113
Cork, Ireland                        |