[Bug-gsl] lmder

[Butel Rémi]

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         *
************************************************