[Bug-gsl] Possible bug in cholesky.c: gsl_linalg_cholesky_decomp()

[Stephen Milborrow]

gsl_linalg_cholesky_decomp() appears to have a bug which causes it take
square roots of negative numbers when passed a non positive-definite matrix.
I believe this is contrary to the intent of the library which is that such
errors should be handled by GSL_ERROR.

The following lines in gsl_linalg_cholesky_decomp are the culprit:

   double A_00 = gsl_matrix_get (A, 0, 0);

   double L_00 = sqrt(A_00);

   if (A_00 <= 0)
     {
       status = GSL_EDOM ;
     }

The code does the check against zero AFTER taking the square root and thus
can try to take the root of a negative A_00 (which is a run time error in
some environments). The correct code is, I think:

   double A_00 = gsl_matrix_get (A, 0, 0);
   double L_00;

   if (A_00 <= 0)
     {
       status = GSL_EDOM ;
     }

   L_00 = sqrt(A_00);

There are similar switched-around checks later on in the same routine, grep
for sqrt to see them. They should be fixed too.

An easy way to reproduce the problem is:

   #include <stdio.h>
   #include "gsl/gsl_linalg.h"

   void main(void)
   {
   double a_data[] = {
       -1, 0,
         0, 1 };

   gsl_matrix_view m = gsl_matrix_view_array (a_data, 2, 2);

   gsl_linalg_cholesky_decomp (&m.matrix);  /* causes sqrt domain err */
   }


This bug report is for the released version:

  gsl-1.6.tar.gz 1/4/2005 ftp.gnu.org/gnu/gsl

but I looked at the source code in the CVS archive and that has the same
problem.

You'll only see this problem if you try to do a cholesky on a non pos-def
matrix, as far as I can tell.  Therefore it's arguable if it even matters I
suppose. I only saw it because I tweaked the routine to use it as a test for
positive definiteness.

Many thanks for the excellent library, it has been very useful for me.

Regards,
Stephen Milborrow.
address@hidden