[Top][All Lists]
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Re: [Bug-gsl] singular Jacobian not detectedby Newton
From: |
Brian Gough |
Subject: |
Re: [Bug-gsl] singular Jacobian not detectedby Newton |
Date: |
Mon, 24 Jan 2011 12:42:30 +0000 |
User-agent: |
Wanderlust/2.15.6 (Almost Unreal) Emacs/23.2 Mule/6.0 (HANACHIRUSATO) |
At Fri, 21 Jan 2011 15:44:53 +0100,
Michel Kern wrote:
> I'm solving small dense (and sometimes nasty) non-linear systems coming
> from chemical equilibrium. Sometimes the solvers fail, but that's life.
> Depending on the system, different methods perform better.
> We found a problem with Newton's method: sometimes we encounter a
> singular Jacobian. This is correctly detected by the LU solve routine
> gsl_linalg_LU_solve, but the return code is ignored by newton_iterate
> (line 95 in newton.c).
>
> It would be useful if this routine checked the error status of the
> functions it calls, so the user can do something about it in the calling
> program.
> It may require an additional error number, like E_SINGJAC or something
> similar. Is that feasible ?
Thanks for the bug report and suggestion. For simplicity I will
propagate the error code from the LU_solve function (GSL_EDOM in this
case).
--
Brian Gough
=== modified file 'multiroots/newton.c'
--- multiroots/newton.c 2007-09-10 18:09:35 +0000
+++ multiroots/newton.c 2011-01-24 12:34:04 +0000
@@ -106,7 +106,12 @@
gsl_linalg_LU_decomp (state->lu, state->permutation, &signum);
- gsl_linalg_LU_solve (state->lu, state->permutation, f, dx);
+ {
+ int status = gsl_linalg_LU_solve (state->lu, state->permutation, f, dx);
+
+ if (status)
+ return status;
+ }
for (i = 0; i < n; i++)
{