bugGNU Scientific Library - Bugs: bug #42220, Division by zero in...

 
 

You are not allowed to post comments on this tracker with your current authentication level.

bug #42220: Division by zero in "hybrid*" when initial guess is root

Submitter:  Patrick Alken <psa>
Submitted:  Mon 28 Apr 2014 03:43:33 PM UTC
   
 
Category:  Runtime error Severity:  3 - Normal
Operating System:  Status:  None
Assigned to:  None Open/Closed:  Open
Release: 

Fri 13 Jun 2014 09:03:22 AM UTC, comment #2: 

Hi!

My last psot was somewhat incomplete and did not contain all information. This is th second try (hope i' doing better this time):

I am using GSL 1.15 on Ububntu 12.04 LTS 64 Bit with gcc 4.6.3.

The example program provided by Curran runs without problems on my machine. If un-comment line 106 in bug_hybrid.c the test fails for gsl_multiroot_fsolver_broyden.

Looking into the problem I found that function broyden_iterate calls GSL_ERROR if lambda=p*v becomes zero. This may happen if the initial guess is the root of the system or even if during iteration the current guess is too close to the root of the system.

My first (quick and dirty) solution was to return GSL_EZERODIV instead of calling GSL_ERROR and checking in the driving loop if the current guess is a root of the system.

My suggestion for solving this problem a bit more satisfying is to pass the desired precision via the vstate argument to the function broyden_iterate. Then the function could check if it has reached a root and stop doing further calculations. I think, this is the way the gsl_ode* functions do it.

As far as I can see, this problem affects all the multiroot solvers.

Kind regards

Oliver

Oliver Ofenloch <ofi>
Thu 12 Jun 2014 01:07:04 PM UTC, comment #1: 

Hi!

In fact lambda = p.v becomes zero if the last x vector matches the root of the system too closely.

The problem can be solved by returning GSL_EZERODIV in function broyden_iterate (file broden.c, line 368, gsl-1.15) instead of calling GSL_ERROR ("approximation to Jacobian has collapsed", GSL_EZERODIV).

The iteration loop then hast to check if the last x vector is a root of the system:

do {
  ++iter;
  status = gsl_multiroot_fsolver_iterate(solver);
  if (status != GSL_SUCCESS) {
    status = gsl_multiroot_test_residual(gsl_multiroot_fsolver_f(solver),
      EPSABS);
    break;
  }
  status = gsl_multiroot_test_residual(gsl_multiroot_fsolver_f(solver),
    EPSABS);
} while ((status == GSL_CONTINUE) && (iter < 1000));


Sorry for the bad formatting, but I'm new here. If you have any questinos, please feel free to contact me.

Best regards

Oliver

Oliver Ofenloch <ofi>
Mon 28 Apr 2014 03:43:33 PM UTC, original submission:  

Hello GSL team,

The GSL implementations of the "Hybrid" family of multidimensional
rootfinders from MINPACK contain a bug where they divide by zero when the
initial guess happens to be the root.  This affects the "fsolver" types
"hybrid" and "hybrids" as well as the "fdfsolver" types "hybridj" and
"hybridsj".  This results in a floating point exception if trapping is
enabled.

I have attached a small program (bug_hybrid.c), based on the manual's
example code, that triggers this bug when GSL_IEEE_MODE=trap-common.  The
current output is:

    GSL_IEEE_MODE="trap-common"
    type = hybrid
    Floating point exception (core dumped)

The expected output is:

    GSL_IEEE_MODE="trap-common"
    type = hybrid
      status = success
      root = 1.000000e+00
    type = hybrids
      status = success
      root = 1.000000e+00
    type = hybridj
      status = success
      root = 1.000000e+00
    type = hybridsj
      status = success
      root = 1.000000e+00


I'll note that when I originally discovered this bug, a colleague of mine
claimed that the problem is not present in the original FORTRAN
implementation in MINPACK, but rather was introduced in the GSL port.

Thank you for looking into this, and if you have any questions, just let me
know.

 - Curran

Patrick Alken <psa>
Group administrator

 

(Note: upload size limit is set to 16384 kB, after insertion of the required escape characters.)

Attached Files
file #31266:  bug_hybrid.c added by psa (3KiB - text/x-csrc)

 

Depends on the following items: None found

Items that depend on this one: None found

 

Carbon-Copy List
  • -email is unavailable- added by ofi (Posted a comment)
  • -email is unavailable- added by psa (Submitted the item)
  •  

    Follows 1 latest change.

    Date Changed by Updated Field Previous Value => Replaced by
    2014-04-28 psa Attached File- Added bug_hybrid.c, #31266

    Back to the top

    Powered by Savane 3.13-758e.
    Corresponding source code