[Ipopt] bug in finite difference Jacobian

[Enrico Bertolazzi]

I wish to signal an error in finite difference Jacobian computation.
In line 2388 of file src/Interfaces/lpTNLPAdapter.cpp there is this piece of code:

        // Compute the finite difference Jacobian
        for (Index ivar = 0; ivar<n_full_x_; ivar++) {
          if (findiff_x_l_[ivar] < findiff_x_u_[ivar]) {
            const Number xorig = full_x_pert[ivar];
            Number this_perturbation =
              findiff_perturbation_*Max(1., fabs(full_x_[ivar]));
            full_x_pert[ivar] += this_perturbation;
            if (full_x_pert[ivar] > findiff_x_u_[ivar]) {
              full_x_pert[ivar] = xorig - this_perturbation;
            }
            retval = tnlp_->eval_g(n_full_x_, full_x_pert, true, n_full_g_,
                                   full_g_pert);
            if (!retval) break;
            for (Index i=findiff_jac_ia_[ivar]; i<findiff_jac_ia_[ivar+1]; i++) {
              const Index& icon = findiff_jac_ja_[i];
              const Index& ipos = findiff_jac_postriplet_[i];
              jac_g_[ipos] =
                (full_g_pert[icon]-full_g_[icon])/this_perturbation;
            }
            full_x_pert[ivar] = xorig;
          }
        }

the code (al line 2395) when the test full_x_pert[ivar] > findiff_x_u_[ivar]
is true change the increment (computed at line 2392)

full_x_pert[ivar] += this_perturbation ; (= xorig + this_perturbation)

to (see line 2396)

full_x_pert[ivar] = xorig - this_perturbation;

but the finite difference approximation (at line 2404) is 

jac_g_[ipos] = (full_g_pert[icon]-full_g_[icon])/this_perturbation;

thus the sign of the ratio is wrong.
In practice when finite difference are computed on a point touching
the limit the computed value has sign changed.

A simple workaround id to change line 2396

              full_x_pert[ivar] = xorig - this_perturbation;

with

              this_perturbation = -this_perturbation ; full_x_pert[ivar] = xorig + this_perturbation;

can some developer fix the error?


Enrico Bertolazzi
__________________________________________________________________________
Dipartimento di Ingegneria Industriale             
Università degli Studi di Trento
http://www.ing.unitn.it/~bertolaz