[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

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

reply via email to

[Prev in Thread] Current Thread [Next in Thread]