Re: [Bug-gsl] singular Jacobian not detectedby Newton

[Brian Gough]

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++)
     {