[Top][All Lists]
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[Bug-gsl] Possible bug in cholesky.c: gsl_linalg_cholesky_decomp()
From: |
Stephen Milborrow |
Subject: |
[Bug-gsl] Possible bug in cholesky.c: gsl_linalg_cholesky_decomp() |
Date: |
Tue, 2 Aug 2005 17:44:40 -0700 |
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
- [Bug-gsl] Possible bug in cholesky.c: gsl_linalg_cholesky_decomp(),
Stephen Milborrow <=