[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] = {
-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
- Re: [Bug-gsl] GSL zsolve.c bug and the fix,
Suleyman Guleyupoglu <=