At Wed, 24 Mar 2010 12:38:53 -0600,
Daniel Neilson wrote:
> I have an application in which I'm using gsl_linalg_cholesky_decomp()
> to decompose a symmetric positive-definite matrix, A. However, the
> decomp fails in a very small percentage of the matrices that I feed to
> this function due to round-off error and all of that fun stuff.
>
> Does anyone happen to know of a good method that will still give me
> the L L^T decomposition for A when one of these errors is detected? For
> example, is there a method whereby I can add a small-magnitude diagonal
> matrix to A, Cholesky decomp that, and modify the result to obtain the
> decomp as though I'd done it with A instead of the modified A?
In linalg/cholesky.c there is a function which handles all the square roots
double
quiet_sqrt (double x)
/* avoids runtime error, for checking matrix for positive definiteness
*/
{
return (x >= 0) ? sqrt(x) : GSL_NAN;
}
You make a copy of cholesky.c in your application and modify this to
return 0 if x is negative but "close enough" to zero.
--
Brian Gough
GNU Scientific Library -
http://www.gnu.org/software/gsl/
_______________________________________________
Help-gsl mailing list
[email protected]
http://lists.gnu.org/mailman/listinfo/help-gsl
