[Top][All Lists]

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

Re: [Bug-gsl] GSL zsolve.c bug and the fix

From: Suleyman Guleyupoglu
Subject: Re: [Bug-gsl] GSL zsolve.c bug and the fix
Date: Mon, 1 May 2006 12:40:38 -0400

----- Original Message ----- From: "Brian Gough" <address@hidden>
To: "Suleyman Guleyupoglu" <address@hidden>
Cc: <address@hidden>
Sent: Monday, April 24, 2006 2:49 PM
Subject: Re: [Bug-gsl] GSL zsolve.c bug and the fix

Suleyman Guleyupoglu writes:
> According to gsl_poly_complex_solve API description, "The function
> returns GSL_SUCCESS if all the roots are found and GSL_EFAILED if
> the QR reduction does not converge."  But this is not true. In
> zsolve.c, the default error handler is called instead of the
> function returning GSL_EFAILED when QR reduction does not
> converge. This causes the program to exit as that is the behavior
> of the default error handler.

Thanks for the bug report.  I've amended the documentation to note
that the error handler is called at that point -- it is generally
called whenever there is an error.

I report this as a bug with the code not the documentation though. I would argue that the documentation should remain the same but the code is changed as I suggested earlier. Otherwise, the function's return value is useless. It either returns GSL_SUCCESS or calls the error handler, in which case the program aborts. So for the caller, what is the point of checking the return value? I.e., as long as the function returns, it was successful. So maybe the function should have been declared "void gsl_poly...."

Further support for the need for the function to return GSL_EFAILED instead of calling error handlers comes from the "Error Handler" documentation. The first paragraph says: "The default behavior of the GSL error handler is to print a short message and call abort(). When this default is in use programs will stop with a core-dump whenever a library routine reports an error. This is intended as a fail-safe default for programs which do not check the return status of library routines (we don't encourage you to write programs this way)."

Even though it is not encouraged, I could use error handlers to deal with this problem but it leads to spaghetti code. You need to have the handler figure out where the error occurred and handle it (even though it does not have all the information at hand) and then call back the approprite section of the code that originated the problem. This is much harder than simply checking the function return value and dealing with the problem there.

Here is an example code inspired by the example code in the documentation. This code causes the default error handler to be called. However, if the gsl library is modified as I suggested in the beginning of this thread, it is very easy for the client code (library user) to check the return value and adjust parameters to get answers. This code illustrates that addition of a small number (1.E-15) to one of the polynomial coefiicients can push gsl_poly_complex_solve to generate an answer instead of failing. Sometimes more than one coefficient needs to be fudged to generate an answer though.

The code is followed by the program output (run on CentOS 4.3) below and it is also attached:

#include <stdio.h>
#include <gsl/gsl_poly.h>
#include <gsl/gsl_errno.h>

#define polyOrderPlus1  6

/* This is a test program to show failure of gsl_poly_complex_solve method.
* According to the documentation, the method should return failure code
* but it does call default handler and exit the application instead. This
* code relies on the GSL_EFAILED return to fudge the coefficients to find
* roots.
* -- Suleyman G

int main (void) {
  int i, status, arrayIndex, failCount = 0;
  const int polyOrder = polyOrderPlus1 - 1;

 double a[polyOrderPlus1] = {

  double z[2*polyOrder];

gsl_poly_complex_workspace * w = gsl_poly_complex_workspace_alloc (polyOrderPlus1);

while (status=gsl_poly_complex_solve(a, polyOrderPlus1, w, z) != GSL_SUCCESS) {
       printf("gsl_poly_complex_solve FAILED for coefficients...\n");
       for (i = 0; i < polyOrderPlus1; i++) {
           printf ("  a[%d]=%+.15lf\n", i, a[i]);

       if (++failCount > 10) {
       } else {
           arrayIndex = (failCount-1) % polyOrderPlus1;
           a[arrayIndex] += 1e-15;

   gsl_poly_complex_workspace_free (w);

if (status == GSL_SUCCESS) {
 printf("gsl_poly_complex_solve FOUND the roots...\n  Coefficients:\n");
    for (i = 0; i < polyOrderPlus1; i++) {
  printf ("    a[%d] = %+.15lf\n", i, a[i]);
    printf("  Roots:\n");
    for (i = 0; i < polyOrder; i++) {
  printf ("    z[%d] = %+.8f %+.8fi\n", i, z[2*i], z[2*i+1]);

  return 0;

gsl_poly_complex_solve FAILED for coefficients...
gsl_poly_complex_solve FOUND the roots...
   a[0] = -0.058690412428658
   a[1] = +0.028123386142471
   a[2] = +0.008523354273975
   a[3] = -0.002639648481342
   a[4] = -0.000378023608125
   a[5] = +0.000099773554248
   z[0] = +1.73113792 +0.00000000i
   z[1] = -3.63775203 +1.16095773i
   z[2] = -3.63775203 -1.16095773i
   z[3] = +4.66659092 +1.23569845i
   z[4] = +4.66659092 -1.23569845i

Attachment: zsolveTest.c
Description: Text document

Attachment: out.txt
Description: Text document

reply via email to

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