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.
[EMAIL PROTECTED]
_______________________________________________
Bug-gsl mailing list
[email protected]
http://lists.gnu.org/mailman/listinfo/bug-gsl
