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

[Suleyman Guleyupoglu]

----- 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] = {
  -0.058690412428659,
  0.028123386142471,
  0.008523354273975,
 -0.002639648481342,
 -0.000378023608125,
     0.000099773554248
  };

  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) {
           break;
       } 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...
 a[0]=-0.058690412428659
 a[1]=+0.028123386142471
 a[2]=+0.008523354273975
 a[3]=-0.002639648481342
 a[4]=-0.000378023608125
 a[5]=+0.000099773554248
gsl_poly_complex_solve FOUND the roots...
 Coefficients:
   a[0] = -0.058690412428658
   a[1] = +0.028123386142471
   a[2] = +0.008523354273975
   a[3] = -0.002639648481342
   a[4] = -0.000378023608125
   a[5] = +0.000099773554248
 Roots:
   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

zsolveTest.c
Description: Text document

out.txt
Description: Text document