[Top][All Lists]
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[Bug-gsl] lmder
From: |
Butel Rémi |
Subject: |
[Bug-gsl] lmder |
Date: |
Mon, 19 Jan 2004 18:02:28 +0100 |
I discovered with interest the existence of GSL,
and started to use it for optimisation problems I have to solve.
The version used is 1.4
I am sorry, I found 4 bugs which doesnt seem to be in the GSL-database.
I'm unfortunatly unable to send you the contexts in which they appears,
but I can send the corrections done to the sources to work around them.
These are the following :
========================================================================
In lmiterate.c, line 42, add :
>> if (state->fnorm == 0.0) {
>> return GSL_EZERODIV;
>> };
The function I try to minimize is awfully complicated, and generates in some
cases
(in some cases only) a null fnorm.
========================================================================
At lmiterate.c, line 55, change into :
>> status = lmpar (r, perm, qtf, diag, state->delta, &(state->par), newton,
>> gradient, sdiag, dx, w);
>>
>> if (status) {
>> return status;
>> };
The status must be checked (as prescribed in section 3.3 of the excellent
documentation of GSL).
========================================================================
At vector_source.c, line 36, I change the return statement by the following
piece of code,
assuming a ''BASE result;'' a few lines before.
>> #if defined(BASE_DOUBLE)
>> (void)signal(SIGFPE, SIG_IGN);
>> result = *(BASE *) (v->data + MULTIPLICITY * i * v->stride);
>> (void)signal(SIGFPE, SIG_IGN);
>> if (result > 1.0E+150 || result < -1.0E+150)
>> {
>> fprintf(stderr, "ERROR : my_gsl_vector_get : v=%x, i=%u, result=%g\n",
>> v, i, result);
>> }
>> return result;
>> #else
>> return *(BASE *) (v->data + MULTIPLICITY * i * v->stride);
>> #endif
This was done to work around a SIGFPE exception; gdb reported that this
exception happens in the return
statement. I'm not sure that my solution is efficient, as gsl_vector_get is a
heavily used function (or macro).
So the bug was eliminated at an upper level, in
lmpar.c at line 244 :
>> {
>> double vmin = gsl_vector_min(newton);
>> double vmax = gsl_vector_max(newton);
>> if (vmin < -1.0E+150 || vmax > 1.0E+150) {
>> printf ("newton = ");
>> gsl_vector_fprintf (stdout, newton, "%g");
>> printf ("\n");
>> return GSL_EOVRFLW;
>> };
>> }
Again, the function optimized is very complex, and seems to generate
instabilities.
As my program is designed to work around the cases (rare) of non convergence
from the GSL library,
I don not worry about an occasionnal failure, if it lets the user able to take
action for repairing,
that is not the case when SIGFPE is raised (or GSL_ERROR, or abort());
the "try { } catch" construct of Java is the correct one for handling exception
!
The way to implement this construct may use the long jump of C ?
I hope you will find usefull these remarks.
And thanks for GSL !
Sincerely yours,
R.Butel
Software Engineer
************************************************
* Rémi BUTEL *
* *
* D.G.O. / E.P.O.C. / C.N.R.S. *
* Avenue des Universités *
* 33405 TALENCE *
* FRANCE *
* *
* tél: (33) 5.40.00.88.31 *
* fax: (33) 5.56.84.08.48 *
* e-mail: address@hidden *
************************************************
- [Bug-gsl] lmder,
Butel Rémi <=