Changeset 424


Ignore:
Timestamp:
Aug 4, 2005 2:55:26 AM (14 years ago)
Author:
andreasw
Message:

Synchronized with branches/dev r423

Location:
trunk
Files:
3 deleted
114 edited
8 copied

Legend:

Unmodified
Added
Removed
  • trunk/AUTHORS

    • Property svn:keywords set to Author Date Id Revision
  • trunk/Algorithm/IpAlgBuilder.cpp

    r336 r424  
    1515#include "IpFilterLineSearch.hpp"
    1616#include "IpMonotoneMuUpdate.hpp"
    17 #include "IpNonmonotoneMuUpdate.hpp"
     17#include "IpAdaptiveMuUpdate.hpp"
    1818#include "IpLoqoMuOracle.hpp"
    1919#include "IpProbingMuOracle.hpp"
     
    6969      "Update strategy for barrier parameter.",
    7070      "monotone",
    71       "monotone", "use a monotone (Fiacco-McCormick) strategy",
    72       "nonmonotone", "use a nonmonotone update strategy",
    73       "Determines whether a nonmonotone barrier "
    74       "parameter strategy is to be used.");
     71      "monotone", "use the monotone (Fiacco-McCormick) strategy",
     72      "adaptive", "use the adaptive update strategy",
     73      "Determines which barrier parameter strategy is to be used.");
    7574    roptions->AddStringOption3(
    7675      "mu_oracle",
    77       "Oracle for a new barrier parameters in the nonmonotone strategy",
     76      "Oracle for a new barrier parameters in the adaptive strategy",
    7877      "probing",
    7978      "probing", "Mehrotra's probing heuristic",
     
    8180      "quality_function", "minimize a quality function",
    8281      "Determines how a new barrier parameter is computed in each "
    83       "\"free-mode\" iteration of the nonmonotone barrier parameter "
    84       "strategy. (Only considered if \"nonmonotone\" is selected for "
     82      "\"free-mode\" iteration of the adaptive barrier parameter "
     83      "strategy. (Only considered if \"adaptive\" is selected for "
    8584      "option \"mu_strategy\".");
    8685    roptions->AddStringOption4(
     
    9392      "average_compl", "base on current average complementarity",
    9493      "Determines how the first value of the barrier parameter should be "
    95       "computed when switching to the \"monotone mode\" in the nonmonotone "
    96       "strategy. (Only considered if \"nonmonotone\" is selected for option "
     94      "computed when switching to the \"monotone mode\" in the adaptive "
     95      "strategy. (Only considered if \"adaptive\" is selected for option "
    9796      "\"mu_strategy\".");
    9897    roptions->AddStringOption2(
     
    186185      new PDFullSpaceSolver(*resto_AugSolver);
    187186
     187    // Convergence check in the restoration phase
     188    SmartPtr<RestoFilterConvergenceCheck> resto_convCheck =
     189      new RestoFilterConvergenceCheck();
     190
    188191    // Line search method for the restoration phase
    189192    SmartPtr<RestoRestorationPhase> resto_resto =
    190193      new RestoRestorationPhase();
    191194    SmartPtr<FilterLineSearch> resto_LineSearch =
    192       new FilterLineSearch(GetRawPtr(resto_resto), GetRawPtr(resto_PDSolver));
     195      new FilterLineSearch(GetRawPtr(resto_resto), GetRawPtr(resto_PDSolver),
     196                           GetRawPtr(resto_convCheck));
    193197
    194198    // Create the mu update that will be used by the restoration phase
     
    200204    std::string resto_smuoracle;
    201205    std::string resto_sfixmuoracle;
    202     if (resto_smuupdate=="nonmonotone" ) {
     206    if (resto_smuupdate=="adaptive" ) {
    203207      options.GetValue("mu_oracle", resto_smuoracle, "resto."+prefix);
    204208      options.GetValue("fixed_mu_oracle", resto_sfixmuoracle, "resto."+prefix);
     
    208212      resto_MuUpdate = new MonotoneMuUpdate(GetRawPtr(resto_LineSearch));
    209213    }
    210     else if (resto_smuupdate=="nonmonotone") {
     214    else if (resto_smuupdate=="adaptive") {
    211215      SmartPtr<MuOracle> resto_MuOracle;
    212216      if (resto_smuoracle=="loqo") {
     
    233237      }
    234238      resto_MuUpdate =
    235         new NonmonotoneMuUpdate(GetRawPtr(resto_LineSearch),
    236                                 resto_MuOracle, resto_FixMuOracle);
    237     }
    238 
    239     // Convergence check in the restoration phase
    240     SmartPtr<RestoFilterConvergenceCheck> resto_convCheck =
    241       new RestoFilterConvergenceCheck();
     239        new AdaptiveMuUpdate(GetRawPtr(resto_LineSearch),
     240                             resto_MuOracle, resto_FixMuOracle);
     241    }
    242242
    243243    // Initialization of the iterates for the restoration phase
     
    268268    // Create the line search to be used by the main algorithm
    269269    SmartPtr<FilterLineSearch> lineSearch =
    270       new FilterLineSearch(GetRawPtr(resto_phase), GetRawPtr(PDSolver));
     270      new FilterLineSearch(GetRawPtr(resto_phase), GetRawPtr(PDSolver),
     271                           convCheck);
    271272
    272273    // The following cross reference is not good: We have to store a
     
    282283    std::string smuoracle;
    283284    std::string sfixmuoracle;
    284     if (smuupdate=="nonmonotone" ) {
     285    if (smuupdate=="adaptive" ) {
    285286      options.GetValue("mu_oracle", smuoracle, prefix);
    286287      options.GetValue("fixed_mu_oracle", sfixmuoracle, prefix);
     
    290291      MuUpdate = new MonotoneMuUpdate(GetRawPtr(lineSearch));
    291292    }
    292     else if (smuupdate=="nonmonotone") {
     293    else if (smuupdate=="adaptive") {
    293294      SmartPtr<MuOracle> muOracle;
    294295      if (smuoracle=="loqo") {
     
    314315        FixMuOracle = NULL;
    315316      }
    316       MuUpdate = new NonmonotoneMuUpdate(GetRawPtr(lineSearch),
    317                                          muOracle, FixMuOracle);
     317      MuUpdate = new AdaptiveMuUpdate(GetRawPtr(lineSearch),
     318                                      muOracle, FixMuOracle);
    318319    }
    319320
  • trunk/Algorithm/IpAlgStrategy.hpp

    r2 r424  
    4646    virtual ~AlgorithmStrategyObject()
    4747    {}
    48     //@}
    49 
    50     /** @name Exceptions */
    51     //@{
    52     /** Exception FAILED_INITIALIZATION for problem during
    53      *  initialization of a strategy object.  This is thrown by a
    54      *  strategy object, if a problem arises during initialization,
    55      *  such as a value out of a feasible range.
    56      */
    57     DECLARE_STD_EXCEPTION(FAILED_INITIALIZATION);
    5848    //@}
    5949
  • trunk/Algorithm/IpConvCheck.hpp

    r2 r424  
    3737      CONTINUE,
    3838      CONVERGED,
     39      CONVERGED_TO_ACCEPTABLE_POINT,
    3940      MAXITER_EXCEEDED,
    4041      FAILED
     
    4748    /** Pure virtual method for performing the convergence test */
    4849    virtual ConvergenceStatus CheckConvergence()=0;
     50
     51    /** Method for testing if the current iterate is considered to
     52     *  satisfy the "accptable level" of accuracy.  The idea is that
     53     *  if the desired convergence tolerance cannot be achieved, the
     54     *  algorithm might stop after a number of acceptable points have
     55     *  been encountered. */
     56    virtual bool CurrentIsAcceptable()=0;
    4957
    5058  private:
  • trunk/Algorithm/IpDefaultIterateInitializer.cpp

    r336 r424  
    4646      "than this value, it is discarded and all constraint multipliers are "
    4747      "set to zero.  This options is also used in the classes "
    48       "\"RestoIterateInitializer\" and \"MinC_1NrmRestorationPhase\".  "
    49       "In the latter it determines when the "
    50       "least-square estimate after returning from the restoration phase is "
    51       "to be discareded.");
     48      "\"RestoIterateInitializer\".  By default, "
     49      "\"resto.constr_mult_init_max\" (the one "
     50      "used in RestoIterateInitializer) is set to zero.");
    5251    reg_options->AddLowerBoundedNumberOption(
    5352      "bound_mult_init_val",
     
    9089
    9190    IpData().InitializeDataStructures(IpNLP(), true, false, false,
    92                                       false, false, false, false);
     91                                      false, false);
    9392
    9493    // get a container of the current point. We will modify parts of this
     
    105104    // Push the x iterates sufficiently inside the bounds
    106105    // Calculate any required shift in x0 and s0
    107     const double dbl_min = std::numeric_limits<double>::min();
    108     const double tiny_double = 100.0*dbl_min;
    109     SmartPtr<const Vector> x = iterates->x();
    110     SmartPtr<const Vector> x_L = IpNLP().x_L();
    111     SmartPtr<const Vector> x_U = IpNLP().x_U();
    112     SmartPtr<const Matrix> Px_L = IpNLP().Px_L();
    113     SmartPtr<const Matrix> Px_U = IpNLP().Px_U();
    114 
    115     SmartPtr<Vector> tmp = x->MakeNew();
    116     SmartPtr<Vector> tmp_l = x_L->MakeNew();
    117     SmartPtr<Vector> tmp_u = x_U->MakeNew();
    118     SmartPtr<Vector> tiny_l = x_L->MakeNew();
    119     tiny_l->Set(tiny_double);
    120 
    121     // Calculate p_l
    122     SmartPtr<Vector> q_l = x_L->MakeNew();
    123     SmartPtr<Vector> p_l = x_L->MakeNew();
    124 
    125     DBG_PRINT_MATRIX(2,"Px_L", *Px_L);
    126     DBG_PRINT_VECTOR(2, "x_L", *x_L);
    127     DBG_PRINT_VECTOR(2, "tmp", *tmp);
    128     Px_L->MultVector(1.0, *x_L, 0.0, *tmp);
    129     Px_U->TransMultVector(1.0, *tmp, 0.0, *tmp_u);
    130     tmp_u->AddOneVector(1., *x_U, -1.);
    131     Px_U->MultVector(1.0, *tmp_u, 0.0, *tmp);
    132     Px_L->TransMultVector(1.0, *tmp, 0.0, *q_l);
    133     q_l->AddOneVector(-1.0, *tiny_l, bound_frac_);
    134 
    135     tmp_l->Set(1.0);
    136     p_l->Copy(*x_L);
    137     p_l->ElementWiseSgn();
    138     p_l->ElementWiseMultiply(*x_L);
    139     p_l->ElementWiseMax(*tmp_l);
    140     p_l->AddOneVector(-1.0, *tiny_l, bound_push_);
    141 
    142     q_l->ElementWiseReciprocal();
    143     p_l->ElementWiseReciprocal();
    144 
    145     p_l->ElementWiseMax(*q_l);
    146     p_l->ElementWiseReciprocal();
    147     p_l->Axpy(1.0, *tiny_l);
    148 
    149     // Calculate p_u
    150     SmartPtr<Vector> q_u = x_U->MakeNew();
    151     SmartPtr<Vector> p_u = x_U->MakeNew();
    152     SmartPtr<Vector> tiny_u = x_U->MakeNew();
    153     tiny_u->Set(tiny_double);
    154 
    155     Px_U->MultVector(1.0, *x_U, 0.0, *tmp);
    156     Px_L->TransMultVector(1.0, *tmp, 0.0, *tmp_l);
    157     tmp_l->Axpy(-1.0, *x_L);
    158     Px_L->MultVector(1.0, *tmp_l, 0.0, *tmp);
    159     Px_U->TransMultVector(1.0, *tmp, 0.0, *q_u);
    160     q_u->AddOneVector(-1.0, *tiny_u, bound_frac_);
    161     DBG_PRINT_VECTOR(2,"q_u",*q_u);
    162 
    163     tmp_u->Set(1.0);
    164     p_u->Copy(*x_U);
    165     p_u->ElementWiseSgn();
    166     p_u->ElementWiseMultiply(*x_U);
    167     p_u->ElementWiseMax(*tmp_u);
    168     p_u->AddOneVector(-1.0, *tiny_u, bound_push_);
    169     DBG_PRINT_VECTOR(2,"p_u",*p_u);
    170 
    171     q_u->ElementWiseReciprocal();
    172     p_u->ElementWiseReciprocal();
    173 
    174     p_u->ElementWiseMax(*q_u);
    175     p_u->ElementWiseReciprocal();
    176     p_u->Axpy(1.0, *tiny_u);
    177     DBG_PRINT_VECTOR(2,"actual_p_u",*p_u);
    178 
    179 
    180     // Calculate the new x
    181     SmartPtr<Vector> delta_x = x->MakeNew();
    182 
    183     SmartPtr<Vector> zero_l = x_L->MakeNew();
    184     zero_l->Set(0.0);
    185     SmartPtr<Vector> zero_u = x_U->MakeNew();
    186     zero_u->Set(0.0);
    187 
    188     Px_L->TransMultVector(-1.0, *x, 0.0, *tmp_l);
    189     tmp_l->AddTwoVectors(1.0, *x_L, 1.0, *p_l, 1.);
    190     tmp_l->ElementWiseMax(*zero_l);
    191     Number nrm_l = tmp_l->Amax();
    192     if (nrm_l>0.) {
    193       Px_L->MultVector(1.0, *tmp_l, 0.0, *delta_x);
    194     }
    195     else {
    196       delta_x->Set(0.);
    197     }
    198 
    199     Px_U->TransMultVector(1.0, *x, 0.0, *tmp_u);
    200     tmp_u->AddTwoVectors(-1.0, *x_U, 1.0, *p_u, 1.);
    201     tmp_u->ElementWiseMax(*zero_u);
    202     Number nrm_u = tmp_u->Amax();
    203     if (nrm_u>0.) {
    204       Px_U->MultVector(-1.0, *tmp_u, 1.0, *delta_x);
    205     }
    206     //delta_x->Axpy(-1.0, *tmp);
    207 
    208     SmartPtr<Vector> new_x = delta_x;
    209     new_x->Axpy(1.0, *x);
    210 
    211     iterates->Set_x_NonConst(*new_x);
     106
     107    SmartPtr<const Vector> new_x;
     108    push_variables(Jnlst(), bound_push_, bound_frac_,
     109                   "x", *iterates->x(), new_x, *IpNLP().x_L(),
     110                   *IpNLP().x_U(), *IpNLP().Px_L(), *IpNLP().Px_U());
     111
     112    iterates->Set_x(*new_x);
    212113    IpData().set_trial(iterates);
    213 
    214     if (nrm_l > 0 || nrm_u > 0) {
    215       Jnlst().Printf(J_DETAILED, J_INITIALIZATION, "Moved initial values of x sufficiently inside the bounds.\n");
    216       Jnlst().PrintVector(J_VECTOR, J_INITIALIZATION, "original x", *IpData().curr()->x());
    217       Jnlst().PrintVector(J_VECTOR, J_INITIALIZATION, "new x", *IpData().trial()->x());
    218     }
    219114
    220115    // Calculate the shift in s...
     
    222117    DBG_PRINT_VECTOR(2, "s", *s);
    223118
    224     SmartPtr<const Vector> d_L = IpNLP().d_L();
    225     SmartPtr<const Vector> d_U = IpNLP().d_U();
    226     SmartPtr<const Matrix> Pd_L = IpNLP().Pd_L();
    227     SmartPtr<const Matrix> Pd_U = IpNLP().Pd_U();
    228     DBG_PRINT_VECTOR(2, "d_L", *d_L);
    229     DBG_PRINT_VECTOR(2, "d_U", *d_U);
    230 
    231     // bump the starting point if necessary
    232     tmp = s->MakeNew();
    233     tmp_l = d_L->MakeNew();
    234     tmp_u = d_U->MakeNew();
    235     tiny_l = d_L->MakeNew();
    236     tiny_l->Set(tiny_double);
    237 
    238     // Calculate p_l
    239     q_l = d_L->MakeNew();
    240     p_l = d_L->MakeNew();
    241 
    242     Pd_L->MultVector(1.0, *d_L, 0.0, *tmp);
    243     Pd_U->TransMultVector(1.0, *tmp, 0.0, *tmp_u);
    244     tmp_u->AddOneVector(1., *d_U, -1.);
    245     Pd_U->MultVector(1.0, *tmp_u, 0.0, *tmp);
    246     Pd_L->TransMultVector(1.0, *tmp, 0.0, *q_l);
    247     q_l->AddOneVector(-1.0, *tiny_l, bound_frac_);
    248 
    249     tmp_l->Set(1.0);
    250     p_l->Copy(*d_L);
    251     p_l->ElementWiseSgn();
    252     p_l->ElementWiseMultiply(*d_L);
    253     p_l->ElementWiseMax(*tmp_l);
    254     p_l->AddOneVector(-1.0, *tiny_l, bound_push_);
    255 
    256     q_l->ElementWiseReciprocal();
    257     p_l->ElementWiseReciprocal();
    258 
    259     p_l->ElementWiseMax(*q_l);
    260     p_l->ElementWiseReciprocal();
    261     p_l->Axpy(1.0, *tiny_l);
    262     DBG_PRINT_VECTOR(2, "p_l", *p_l);
    263 
    264     // Calculate p_u
    265     q_u = d_U->MakeNew();
    266     p_u = d_U->MakeNew();
    267     tiny_u = d_U->MakeNew();
    268     tiny_u->Set(tiny_double);
    269 
    270     Pd_U->MultVector(1.0, *d_U, 0.0, *tmp);
    271     Pd_L->TransMultVector(1.0, *tmp, 0.0, *tmp_l);
    272     tmp_l->Axpy(-1.0, *d_L);
    273     Pd_L->MultVector(1.0, *tmp_l, 0.0, *tmp);
    274     Pd_U->TransMultVector(1.0, *tmp, 0.0, *q_u);
    275     q_u->AddOneVector(-1.0, *tiny_u, bound_frac_);
    276 
    277     tmp_u->Set(1.0);
    278     p_u->Copy(*d_U);
    279     p_u->ElementWiseSgn();
    280     p_u->ElementWiseMultiply(*d_U);
    281     p_u->ElementWiseMax(*tmp_u);
    282     p_u->AddOneVector(-1.0, *tiny_u, bound_push_);
    283 
    284     q_u->ElementWiseReciprocal();
    285     p_u->ElementWiseReciprocal();
    286 
    287     p_u->ElementWiseMax(*q_u);
    288     p_u->ElementWiseReciprocal();
    289     p_u->Axpy(1.0, *tiny_u);
    290     DBG_PRINT_VECTOR(2, "p_u", *p_u);
    291 
    292     zero_l = d_L->MakeNew();
    293     zero_l->Set(0.0);
    294     zero_u = d_U->MakeNew();
    295     zero_u->Set(0.0);
    296 
    297     Pd_L->TransMultVector(-1.0, *s, 0.0, *tmp_l);
    298     tmp_l->AddTwoVectors(1.0, *d_L, 1.0, *p_l, 1.);
    299     tmp_l->ElementWiseMax(*zero_l);
    300     nrm_l = tmp_l->Amax();
    301     SmartPtr<Vector> delta_s = s->MakeNew();
    302     if (nrm_l>0.) {
    303       Pd_L->MultVector(1.0, *tmp_l, 0.0, *delta_s);
    304     }
    305     else {
    306       delta_s->Set(0.);
    307     }
    308 
    309     Pd_U->TransMultVector(1.0, *s, 0.0, *tmp_u);
    310     tmp_u->AddTwoVectors(-1.0, *d_U, 1.0, *p_u, 1.);
    311     tmp_u->ElementWiseMax(*zero_u);
    312     nrm_u = tmp_u->Amax();
    313     Pd_U->MultVector(-1.0, *tmp_u, 1.0, *delta_s);
    314 
    315     SmartPtr<Vector> new_s = delta_s;
    316     new_s->Axpy(1.0, *s);
     119    SmartPtr<const Vector> new_s;
     120    push_variables(Jnlst(), bound_push_, bound_frac_,
     121                   "s", *s, new_s, *IpNLP().d_L(),
     122                   *IpNLP().d_U(), *IpNLP().Pd_L(), *IpNLP().Pd_U());
    317123
    318124    iterates = IpData().trial()->MakeNewContainer();
    319     iterates->Set_s_NonConst(*new_s);
    320     if (nrm_l > 0 || nrm_u > 0) {
    321       Jnlst().Printf(J_DETAILED, J_INITIALIZATION,
    322                      "Moved initial values of s sufficiently inside the bounds.\n");
    323       Jnlst().PrintVector(J_VECTOR, J_INITIALIZATION,
    324                           "original s", *s);
    325       Jnlst().PrintVector(J_VECTOR, J_INITIALIZATION,
    326                           "new s", *iterates->s());
    327     }
     125    iterates->Set_s(*new_s);
    328126
    329127    /////////////////////////////////////////////////////////////////////
     
    347145    /////////////////////////////////////////////////////////////////////
    348146
    349     iterates = IpData().trial()->MakeNewContainer();
     147    least_square_mults(Jnlst(), IpNLP(), IpData(), IpCq(),
     148                       eq_mult_calculator_, constr_mult_init_max_);
     149
     150    // upgrade the trial to the current point
     151    IpData().AcceptTrialPoint();
     152
     153    return true;
     154  }
     155
     156  void DefaultIterateInitializer::push_variables(
     157    const Journalist& jnlst,
     158    Number bound_push,
     159    Number bound_frac,
     160    std::string name,
     161    const Vector& orig_x,
     162    SmartPtr<const Vector>& new_x,
     163    const Vector& x_L,
     164    const Vector& x_U,
     165    const Matrix& Px_L,
     166    const Matrix& Px_U)
     167  {
     168    DBG_START_FUN("DefaultIterateInitializer::push_variables",
     169                  dbg_verbosity);
     170    // Calculate any required shift in x0 and s0
     171    const double dbl_min = std::numeric_limits<double>::min();
     172    const double tiny_double = 100.0*dbl_min;
     173
     174    SmartPtr<Vector> tmp = orig_x.MakeNew();
     175    SmartPtr<Vector> tmp_l = x_L.MakeNew();
     176    SmartPtr<Vector> tmp_u = x_U.MakeNew();
     177    SmartPtr<Vector> tiny_l = x_L.MakeNew();
     178    tiny_l->Set(tiny_double);
     179
     180    // Calculate p_l
     181    SmartPtr<Vector> q_l = x_L.MakeNew();
     182    SmartPtr<Vector> p_l = x_L.MakeNew();
     183
     184    DBG_PRINT_VECTOR(2,"orig_x", orig_x);
     185    DBG_PRINT_MATRIX(2,"Px_L", Px_L);
     186    DBG_PRINT_VECTOR(2, "x_L", x_L);
     187    DBG_PRINT_MATRIX(2,"Px_U", Px_U);
     188    DBG_PRINT_VECTOR(2, "x_U", x_U);
     189
     190    Px_L.MultVector(1.0, x_L, 0.0, *tmp);
     191    Px_U.TransMultVector(1.0, *tmp, 0.0, *tmp_u);
     192    tmp_u->AddOneVector(1., x_U, -1.);
     193    Px_U.MultVector(1.0, *tmp_u, 0.0, *tmp);
     194    Px_L.TransMultVector(1.0, *tmp, 0.0, *q_l);
     195    q_l->AddOneVector(-1.0, *tiny_l, bound_frac);
     196    DBG_PRINT_VECTOR(2, "q_l", *q_l);
     197    // At this point, q_l is
     198    // bound_frac * Px_L^T Px_U(x_U - Px_U^T Px_L x_L)  -  tiny_double
     199    // i.e., it is bound_frac*(x_U - x_L) for those components that have
     200    //          two bounds
     201    //       and - tiny_double for those that have only one or no bounds
     202
     203    tmp_l->Set(bound_push);
     204    p_l->AddOneVector(bound_push, x_L, 0.);
     205    p_l->ElementWiseAbs();
     206    p_l->ElementWiseMax(*tmp_l);
     207    // now p_l is bound_push * max(|x_L|,1)
     208
     209    q_l->ElementWiseReciprocal();
     210    p_l->ElementWiseReciprocal();
     211
     212    p_l->ElementWiseMax(*q_l);
     213    p_l->ElementWiseReciprocal();
     214    //    p_l->Axpy(1.0, *tiny_l);  we shouldn't need this
     215
     216    // At this point, p_l is
     217    //  min(bound_push * max(|x_L|,1), bound_frac*(x_U-x_L)  for components
     218    //                                                       with two bounds
     219    //  bound_push * max(|x_L|,1)                            otherwise
     220    // This is the margin we want to the lower bound
     221    DBG_PRINT_VECTOR(1, "p_l", *p_l);
     222
     223    // Calculate p_u
     224    SmartPtr<Vector> q_u = x_U.MakeNew();
     225    SmartPtr<Vector> p_u = x_U.MakeNew();
     226    SmartPtr<Vector> tiny_u = x_U.MakeNew();
     227    tiny_u->Set(tiny_double);
     228
     229    Px_U.MultVector(1.0, x_U, 0.0, *tmp);
     230    Px_L.TransMultVector(1.0, *tmp, 0.0, *tmp_l);
     231    tmp_l->Axpy(-1.0, x_L);
     232    Px_L.MultVector(1.0, *tmp_l, 0.0, *tmp);
     233    Px_U.TransMultVector(1.0, *tmp, 0.0, *q_u);
     234    q_u->AddOneVector(-1.0, *tiny_u, bound_frac);
     235    DBG_PRINT_VECTOR(2,"q_u",*q_u);
     236    // q_u is now the same as q_l above, but of the same dimension as x_L
     237
     238    tmp_u->Set(bound_push);
     239    p_u->Copy(x_U);
     240    p_u->AddOneVector(bound_push, x_U, 0.);
     241    p_u->ElementWiseAbs();
     242    p_u->ElementWiseMax(*tmp_u);
     243    DBG_PRINT_VECTOR(2,"p_u",*p_u);
     244
     245    q_u->ElementWiseReciprocal();
     246    p_u->ElementWiseReciprocal();
     247
     248    p_u->ElementWiseMax(*q_u);
     249    p_u->ElementWiseReciprocal();
     250    p_u->Axpy(1.0, *tiny_u);
     251    // At this point, p_l is
     252    //  min(bound_push * max(|x_U|,1), bound_frac*(x_U-x_L)  for components
     253    //                                                       with two bounds
     254    //  bound_push * max(|x_U|,1)                            otherwise
     255    // This is the margin we want to the upper bound
     256    DBG_PRINT_VECTOR(2,"actual_p_u",*p_u);
     257
     258    // Calculate the new x
     259    SmartPtr<Vector> delta_x = orig_x.MakeNew();
     260
     261    SmartPtr<Vector> zero_l = x_L.MakeNew();
     262    zero_l->Set(0.0);
     263    SmartPtr<Vector> zero_u = x_U.MakeNew();
     264    zero_u->Set(0.0);
     265
     266    Px_L.TransMultVector(-1.0, orig_x, 0.0, *tmp_l);
     267    tmp_l->AddTwoVectors(1.0, x_L, 1.0, *p_l, 1.);
     268    tmp_l->ElementWiseMax(*zero_l);
     269    // tmp_l is now max(x_L + p_l - x, 0), i.e., the amount by how
     270    // much need to correct the variable
     271    Number nrm_l = tmp_l->Amax();
     272    if (nrm_l>0.) {
     273      Px_L.MultVector(1.0, *tmp_l, 0.0, *delta_x);
     274    }
     275    else {
     276      delta_x->Set(0.);
     277    }
     278
     279    Px_U.TransMultVector(1.0, orig_x, 0.0, *tmp_u);
     280    tmp_u->AddTwoVectors(-1.0, x_U, 1.0, *p_u, 1.);
     281    tmp_u->ElementWiseMax(*zero_u);
     282    // tmp_u is now max(x - (x_U-p_u), 0), i.e., the amount by how
     283    // much need to correct the variable
     284    Number nrm_u = tmp_u->Amax();
     285    if (nrm_u>0.) {
     286      Px_U.MultVector(-1.0, *tmp_u, 1.0, *delta_x);
     287    }
     288
     289    if (nrm_l > 0 || nrm_u > 0) {
     290      delta_x->Axpy(1.0, orig_x);
     291      new_x = ConstPtr(delta_x);
     292      jnlst.Printf(J_DETAILED, J_INITIALIZATION, "Moved initial values of %s sufficiently inside the bounds.\n", name.c_str());
     293      jnlst.PrintVector(J_VECTOR, J_INITIALIZATION, "original vars", orig_x);
     294      jnlst.PrintVector(J_VECTOR, J_INITIALIZATION, "new vars", *new_x);
     295    }
     296    else {
     297      new_x = &orig_x;
     298      jnlst.Printf(J_DETAILED, J_INITIALIZATION, "Initial values of %s sufficiently inside the bounds.\n", name.c_str());
     299    }
     300  }
     301
     302  void DefaultIterateInitializer::least_square_mults(
     303    const Journalist& jnlst,
     304    IpoptNLP& ip_nlp,
     305    IpoptData& ip_data,
     306    IpoptCalculatedQuantities& ip_cq,
     307    const SmartPtr<EqMultiplierCalculator>& eq_mult_calculator,
     308    Number constr_mult_init_max)
     309  {
     310    DBG_START_FUN("DefaultIterateInitializer::least_square_mults",
     311                  dbg_verbosity);
     312
     313    SmartPtr<IteratesVector> iterates = ip_data.trial()->MakeNewContainer();
    350314    iterates->create_new_y_c();
    351315    iterates->create_new_y_d();
    352     if (IsValid(eq_mult_calculator_) && constr_mult_init_max_>0.) {
     316
     317    if (IsValid(eq_mult_calculator) && constr_mult_init_max>0. &&
     318        iterates->y_c_NonConst()->Dim()+iterates->y_d_NonConst()->Dim()>0) {
    353319      // First move all the trial data into the current fields, since
    354320      // those values are needed to compute the initial values for
    355321      // the multipliers
    356       IpData().CopyTrialToCurrent();
     322      ip_data.CopyTrialToCurrent();
    357323      SmartPtr<Vector> y_c = iterates->y_c_NonConst();
    358324      SmartPtr<Vector> y_d = iterates->y_d_NonConst();
    359       bool retval = eq_mult_calculator_->CalculateMultipliers(*y_c, *y_d);
     325      bool retval = eq_mult_calculator->CalculateMultipliers(*y_c, *y_d);
    360326      if (!retval) {
    361327        y_c->Set(0.0);
     
    363329      }
    364330      else {
    365         Jnlst().Printf(J_DETAILED, J_INITIALIZATION,
    366                        "Least square estimates max(y_c) = %e, max(y_d) = %e\n",
    367                        y_c->Amax(), y_d->Amax());
     331        /*
     332        {
     333          ip_data.set_trial(iterates);
     334          printf("grad_x = %e grad_s = %e y_c = %e y_d = %e\n",
     335                 ip_cq.trial_grad_lag_x()->Amax(),
     336                 ip_cq.trial_grad_lag_s()->Amax(),
     337                 y_c->Amax(),
     338                 y_d->Amax());
     339          iterates = ip_data.trial()->MakeNewContainer();
     340        }
     341        */
     342        jnlst.Printf(J_DETAILED, J_INITIALIZATION,
     343                     "Least square estimates max(y_c) = %e, max(y_d) = %e\n",
     344                     y_c->Amax(), y_d->Amax());
    368345        Number yinitnrm = Max(y_c->Amax(), y_d->Amax());
    369         if (yinitnrm > constr_mult_init_max_) {
     346        if (yinitnrm > constr_mult_init_max) {
    370347          y_c->Set(0.0);
    371348          y_d->Set(0.0);
    372349        }
     350        else {
     351          ip_data.Append_info_string("y");
     352        }
    373353      }
    374354    }
     
    377357      iterates->y_d_NonConst()->Set(0.0);
    378358    }
    379     IpData().set_trial(iterates);
    380 
    381     //Qu: why do you print curr here? they have not been updated yet?
    382     DBG_PRINT_VECTOR(2, "y_c", *IpData().curr()->y_c());
    383     DBG_PRINT_VECTOR(2, "y_d", *IpData().curr()->y_d());
    384 
    385     DBG_PRINT_VECTOR(2, "z_L", *IpData().curr()->z_L());
    386     DBG_PRINT_VECTOR(2, "z_U", *IpData().curr()->z_U());
    387     DBG_PRINT_VECTOR(2, "v_L", *IpData().curr()->v_L());
    388     DBG_PRINT_VECTOR(2, "v_U", *IpData().curr()->v_U());
    389 
    390     // upgrade the trial to the current point
    391     IpData().AcceptTrialPoint();
    392 
    393     return true;
     359    ip_data.set_trial(iterates);
     360
     361    DBG_PRINT_VECTOR(2, "y_c", *ip_data.trial()->y_c());
     362    DBG_PRINT_VECTOR(2, "y_d", *ip_data.trial()->y_d());
    394363  }
    395364
  • trunk/Algorithm/IpDefaultIterateInitializer.hpp

    r336 r424  
    4848    virtual bool SetInitialIterates();
    4949
     50    /** Auxilliary function for moving the initial point.  This is
     51     *  declared static so that it can also be used from
     52     *  WarmStartIterateInitializer. */
     53    static void push_variables(const Journalist& jnlst,
     54                               Number bound_push,
     55                               Number bound_frac,
     56                               std::string name,
     57                               const Vector& orig_x,
     58                               SmartPtr<const Vector>& new_x,
     59                               const Vector& x_L,
     60                               const Vector& x_U,
     61                               const Matrix& Px_L,
     62                               const Matrix& Px_U);
     63
     64    /** Auxilliary function for computing least_square multipliers.
     65     *  The multipliers are computed based on the values in the trial
     66     *  fields (current is overwritten).  On return, the multipliers
     67     *  are in the trial fields as well.  The value of
     68     *  constr_mult_init_max determines if the computed least square
     69     *  estimate should be used, or if the initial multipliers are set
     70     *  to zero. */
     71    static void least_square_mults(const Journalist& jnlst,
     72                                   IpoptNLP& ip_nlp,
     73                                   IpoptData& ip_data,
     74                                   IpoptCalculatedQuantities& ip_cq,
     75                                   const SmartPtr<EqMultiplierCalculator>& eq_mult_calculator,
     76                                   Number constr_mult_init_max);
     77
     78
    5079    /** Methods for IpoptType */
    5180    //@{
  • trunk/Algorithm/IpFilter.cpp

    r157 r424  
    4040  {
    4141    DBG_START_METH("FilterLineSearch::Filter::Acceptable", dbg_verbosity);
    42     DBG_ASSERT(vals.size()==dim_);
     42    DBG_ASSERT((Index)vals.size()==dim_);
    4343    bool acceptable = true;
    4444    std::list<FilterEntry*>::iterator iter;
     
    5656  {
    5757    DBG_START_METH("FilterLineSearch::Filter::AddEntry", dbg_verbosity);
    58     DBG_ASSERT(vals.size()==dim_);
     58    DBG_ASSERT((Index)vals.size()==dim_);
    5959    std::list<FilterEntry*>::iterator iter;
    6060    for (iter = filter_list_.begin(); iter != filter_list_.end();
  • trunk/Algorithm/IpFilter.hpp

    r324 r424  
    3737    {
    3838      Index ncoor = (Index)vals_.size();
    39       DBG_ASSERT(vals.size() == ncoor);
     39      DBG_ASSERT((Index)vals.size() == ncoor);
    4040
    4141      // ToDo decide where we can get Compare_le from
     
    5757    {
    5858      Index ncoor = (Index)vals_.size();
    59       DBG_ASSERT(vals.size() == ncoor);
     59      DBG_ASSERT((Index)vals.size() == ncoor);
    6060
    6161      bool retval = true;
  • trunk/Algorithm/IpFilterLineSearch.cpp

    r389 r424  
    1010#include "IpJournalist.hpp"
    1111#include "IpRestoPhase.hpp"
     12#include "IpAlgTypes.hpp"
    1213
    1314#ifdef OLD_C_HEADERS
     
    2728
    2829  FilterLineSearch::FilterLineSearch(const SmartPtr<RestorationPhase>& resto_phase,
    29                                      const SmartPtr<PDSystemSolver>& pd_solver)
     30                                     const SmartPtr<PDSystemSolver>& pd_solver,
     31                                     const SmartPtr<ConvergenceCheck>& conv_check)
    3032      :
    3133      LineSearch(),
     34      theta_max_(-1.0),
     35      theta_min_(-1.0),
     36      filter_(2),
    3237      resto_phase_(resto_phase),
    3338      pd_solver_(pd_solver),
    34       theta_min_(-1.0),
    35       theta_max_(-1.0),
    36       filter_(2)
     39      conv_check_(conv_check)
    3740  {
    3841    DBG_START_FUN("FilterLineSearch::FilterLineSearch",
     
    6871      "eta_phi",
    6972      "Relaxation factor in the Armijo condition.",
    70       0.0, true, 0.5, true, 1e-4,
     73      0.0, true, 0.5, true, 1e-8,
    7174      "(See Eqn. (20) in implementation paper)");
    7275    roptions->AddLowerBoundedNumberOption(
     
    201204      "infeasible.  In the filter line search procedure, the restoration "
    202205      "phase is called more qucikly than usually, and more reduction in "
    203       "the constraint violation is enforced.");
     206      "the constraint violation is enforced. If the problem is square, this "
     207      "is enabled automatically.");
    204208    roptions->AddLowerBoundedNumberOption(
    205209      "expect_infeasible_problem_ctol",
     
    208212      "If the constraint violation becomes small than this threshold, "
    209213      "the \"expect_infeasible_problem\" heuristics in the filter line "
    210       "search will are disabled.");
     214      "search will are disabled. If the problem is square, this is set to 0.");
    211215    roptions->AddLowerBoundedNumberOption(
    212216      "soft_resto_pderror_reduction_factor",
     
    221225      "phase is called.  Choosing \"0\" here disables the soft "
    222226      "restoration phase.");
     227    roptions->AddStringOption2(
     228      "start_with_resto",
     229      "Tells algorithm to switch to restoration phase in first iteration.",
     230      "no",
     231      "no", "don't force start in restoration phase",
     232      "yes", "force start in restoration phase",
     233      "Setting this option to yes forces the algorithm to switch to the "
     234      "restoration phase in the first iteration.  If the initial point "
     235      "is feasible, the algorithm will abort with a failure.");
    223236    roptions->AddLowerBoundedNumberOption(
    224237      "tiny_step_tol",
     
    243256      "Determines the number of trial iterations before the watchdog "
    244257      "procedure is aborted and the algorithm returns to the stored point.");
    245     roptions->AddLowerBoundedNumberOption(
    246       "acceptable_tol",
    247       "Threshold for NLP error to consider iterate as acceptable.",
    248       0.0, false, 1e-6,
    249       "Determines tolerance for which an iterate is considered as accepable "
    250       "solution.  If the algorithm would trigger the restoration phase at "
    251       "such a point, it instead terminates, returning this acceptable "
    252       "point.  Further, if the algorithm encounters \"acceptable_iter_max\" "
    253       "successive points satisfying this NLP tolerance, the algorithm "
    254       "terminates.");
    255     roptions->AddLowerBoundedIntegerOption(
    256       "acceptable_iter_max",
    257       "Number of acceptable iterates to trigger termination.",
    258       0, 15,
    259       "if the algorithm encounters so many successive acceptable iterates "
    260       "(see \"acceptable_tol\"), it terminates, assuming that the problem "
    261       "has been solved to best possible accuracy given round-off.");
    262258  }
    263259
     
    294290    alpha_for_y_ = AlphaForYEnum(enum_int);
    295291    options.GetNumericValue("corrector_compl_avrg_red_fact", corrector_compl_avrg_red_fact_, prefix);
     292    options.GetNumericValue("expect_infeasible_problem_ctol", expect_infeasible_problem_ctol_, prefix);
    296293    options.GetBoolValue("expect_infeasible_problem", expect_infeasible_problem_, prefix);
    297     options.GetNumericValue("expect_infeasible_problem_ctol", expect_infeasible_problem_ctol_, prefix);
     294
     295    options.GetBoolValue("start_with_resto", start_with_resto_, prefix);
    298296
    299297    bool retvalue = true;
     
    308306    options.GetIntegerValue("watchdog_trial_iter_max", watchdog_trial_iter_max_, prefix);
    309307    options.GetIntegerValue("watchdog_shortened_iter_trigger", watchdog_shortened_iter_trigger_, prefix);
    310     options.GetNumericValue("acceptable_tol", acceptable_tol_, prefix);
    311     options.GetIntegerValue("acceptable_iter_max", acceptable_iter_max_, prefix);
    312308
    313309    // ToDo decide if also the PDSystemSolver should be initialized here...
     
    320316
    321317    count_successive_shortened_steps_ = 0;
    322     count_acceptable_iter_ = 0;
    323318
    324319    acceptable_iterate_ = NULL;
     
    335330                   IpData().iter_count());
    336331
     332    // If the problem is square, we want to enable the
     333    // expect_infeasible_problem option automatically so that the
     334    // restoration phase is entered soon
     335    if (IpCq().IsSquareProblem()) {
     336      expect_infeasible_problem_ = true;
     337      expect_infeasible_problem_ctol_ = 0.;
     338    }
     339
    337340    // Store current iterate if the optimality error is on acceptable
    338341    // level to restored if things fail later
    339     if (IpCq().curr_nlp_error()<=acceptable_tol_) {
     342    if (CurrentIsAcceptable()) {
    340343      Jnlst().Printf(J_DETAILED, J_LINE_SEARCH,
    341344                     "Storing current iterate as backup acceptable point.\n");
    342       IpData().Append_info_string("A");
    343345      StoreAcceptablePoint();
    344346    }
     
    369371      // off of a something else than the norm of the search direction
    370372      goto_resto = true;
     373    }
     374
     375    if (start_with_resto_) {
     376      // If the use requested to start with the restoration phase,
     377      // skip the lin e search and do exactly that.  Reset the flag so
     378      // that this happens only once.
     379      goto_resto = true;
     380      start_with_resto_= false;
    371381    }
    372382
     
    522532            AugmentFilter();
    523533          }
    524           if (IpCq().curr_nlp_error()<=acceptable_tol_) {
     534          if (CurrentIsAcceptable()) {
    525535            THROW_EXCEPTION(ACCEPTABLE_POINT_REACHED,
    526536                            "Restoration phase called at acceptable point.");
     
    528538
    529539          if (!IsValid(resto_phase_)) {
    530             //ToDo
    531540            THROW_EXCEPTION(IpoptException, "No Restoration Phase given to this Filter Line Search Object!");
    532541          }
     542          // ToDo make the 1e-2 below a parameter?
    533543          if (IpCq().curr_constraint_violation()<=
    534               1e-2*Min(IpData().tol(),IpData().primal_inf_tol())) {
     544              1e-2*IpData().tol()) {
    535545            bool found_acceptable = RestoreAcceptablePoint();
    536546            if (found_acceptable) {
    537547              THROW_EXCEPTION(ACCEPTABLE_POINT_REACHED,
    538                               "Restoration phase called at almost feasible point, but acceptable point could be restore.\n");
     548                              "Restoration phase called at almost feasible point, but acceptable point could be restored.\n");
    539549            }
    540550            else {
     551              // ToDo does that happen too often?
    541552              THROW_EXCEPTION(RESTORATION_FAILED,
    542553                              "Restoration phase called, but point is almost feasible.");
     
    565576          }
    566577          count_successive_shortened_steps_ = 0;
    567           count_acceptable_iter_ = 0;
    568578          if (expect_infeasible_problem_) {
    569579            expect_infeasible_problem_ = false;
     
    584594      PerformDualStep(alpha_primal, alpha_dual_max, actual_delta);
    585595
    586       if (true || n_steps==0) { // The original heuristic only
     596      if (n_steps==0) {
    587597        // accepted this if a full step was
    588598        // taken
    589599        count_successive_shortened_steps_ = 0;
    590         if (acceptable_iter_max_>0) {
    591           if (IpCq().curr_nlp_error()<=acceptable_tol_) {
    592             count_acceptable_iter_++;
    593             if (count_acceptable_iter_>=acceptable_iter_max_) {
    594               IpData().AcceptTrialPoint();
    595               THROW_EXCEPTION(ACCEPTABLE_POINT_REACHED,
    596                               "Algorithm seems stuck at acceptable level.");
    597             }
    598           }
    599           else {
    600             count_acceptable_iter_=0;
    601           }
    602         }
    603600        watchdog_shortened_iter_ = 0;
    604601      }
     
    606603        count_successive_shortened_steps_++;
    607604        watchdog_shortened_iter_++;
    608         count_acceptable_iter_ = 0;
    609605      }
    610606
     
    778774    DBG_START_METH("FilterLineSearch::IsFtype",
    779775                   dbg_verbosity);
    780 
     776    DBG_ASSERT(reference_theta_>0. || reference_gradBarrTDelta_ < 0.0);
    781777    return (reference_gradBarrTDelta_ < 0.0 &&
    782778            alpha_primal_test*pow(-reference_gradBarrTDelta_,s_phi_) >
     
    807803      // evaluation of the functions (in that case, we want to further
    808804      // reduce the step size
    809       Number trial_barr = IpCq().trial_barrier_obj();
    810       Number trial_theta = IpCq().trial_constraint_violation();
     805      /* Number trial_barr = */ IpCq().trial_barrier_obj();
     806      /* Number trial_theta = */
     807      IpCq().trial_constraint_violation();
    811808      return true;
    812809    }
     
    875872  bool FilterLineSearch::ArmijoHolds(Number alpha_primal_test)
    876873  {
     874    /*
     875    Jnlst().Printf(J_DETAILED, J_LINE_SEARCH,
     876                   "ArmijoHolds test with trial_barr = %25.16e reference_barr = %25.16e\n        alpha_primal_test = %25.16e reference_gradBarrTDelta = %25.16e\n", IpCq().trial_barrier_obj(), reference_barr_,alpha_primal_test,reference_gradBarrTDelta_);
     877    */
    877878    return Compare_le(IpCq().trial_barrier_obj()-reference_barr_,
    878879                      eta_phi_*alpha_primal_test*reference_gradBarrTDelta_,
     
    11001101
    11011102    Number alpha_y;
    1102     switch (corrector_type_) {
     1103    switch (alpha_for_y_) {
    11031104      case PRIMAL_ALPHA_FOR_Y:
    11041105      alpha_y = alpha_primal;
     
    11131114      alpha_y = Max(alpha_dual, alpha_primal);
    11141115      break;
     1116      case FULL_STEP_FOR_Y:
     1117      alpha_y = 1;
     1118      break;
    11151119      case MIN_DUAL_INFEAS_ALPHA_FOR_Y:
    11161120      case SAFE_MIN_DUAL_INFEAS_ALPHA_FOR_Y:
     
    11511155    // Set the eq multipliers from the step now that alpha_y
    11521156    // has been calculated.
     1157    DBG_PRINT((1, "alpha_y = %e\n", alpha_y));
     1158    DBG_PRINT_VECTOR(2, "delta_y_c", *delta->y_c());
     1159    DBG_PRINT_VECTOR(2, "delta_y_d", *delta->y_d());
    11531160    IpData().SetTrialEqMultipliersFromStep(alpha_y, *delta->y_c(), *delta->y_d());
    11541161
     
    11711178
    11721179    Number theta_soc_old = 0.;
    1173     Number theta_curr = IpCq().curr_constraint_violation();
    11741180    Number theta_trial = IpCq().trial_constraint_violation();
    11751181    Number alpha_primal_soc = alpha_primal;
     
    12221228        Jnlst().Printf(J_WARNING, J_MAIN, "Warning: SOC step rejected due to evaluation error\n");
    12231229        accept = false;
     1230        // There is no point in continuing SOC procedure
     1231        break;
    12241232      }
    12251233
     
    15861594  }
    15871595
     1596  bool FilterLineSearch::CurrentIsAcceptable()
     1597  {
     1598    return (IsValid(conv_check_) &&
     1599            conv_check_->CurrentIsAcceptable());
     1600  }
     1601
    15881602  void FilterLineSearch::StoreAcceptablePoint()
    15891603  {
  • trunk/Algorithm/IpFilterLineSearch.hpp

    r336 r424  
    1616#include "IpPDSystemSolver.hpp"
    1717#include "IpIpoptType.hpp"
     18#include "IpConvCheck.hpp"
    1819
    1920namespace Ipopt
     
    3233    /** Constructor.  The PDSystemSolver object only needs to be
    3334     *  provided (i.e. not NULL) if second order correction is to be
    34      *  used. */
     35     *  used.  The ConvergenceCheck object is used to determine
     36     *  whether the current iterate is acceptable (for example, the
     37     *  restoration phase is not started if the acceptability level
     38     *  has been reached).  If conv_check is NULL, we assume that the
     39     *  current iterate is not acceptable. */
    3540    FilterLineSearch(const SmartPtr<RestorationPhase>& resto_phase,
    36                      const SmartPtr<PDSystemSolver>& pd_solver
     41                     const SmartPtr<PDSystemSolver>& pd_solver,
     42                     const SmartPtr<ConvergenceCheck>& conv_check
    3743                    );
    3844
     
    242248     *  found. Returns true if such as point is available. */
    243249    bool RestoreAcceptablePoint();
     250
     251    /** Method for determining if the current iterate is acceptable.
     252     *  This is a wrapper for same method from ConvergenceCheck, but
     253     *  returns false, if no ConvergenceCheck object is provided. */
     254    bool CurrentIsAcceptable();
    244255
    245256    /** @name Parameters for the filter algorithm.  Names as in the paper */
     
    280291      MIN_ALPHA_FOR_Y,
    281292      MAX_ALPHA_FOR_Y,
     293      FULL_STEP_FOR_Y,
    282294      MIN_DUAL_INFEAS_ALPHA_FOR_Y,
    283295      SAFE_MIN_DUAL_INFEAS_ALPHA_FOR_Y
     
    334346    Index watchdog_shortened_iter_trigger_;
    335347
    336     /** Acceptable tolerance for the problem to terminate earlier if
    337      *  algorithm seems stuck or cycling */
    338     Number acceptable_tol_;
    339     /** Maximum number of iterations with acceptable level of accuracy
    340      *  and full steps, after which the algorithm terminates.  If 0,
    341      *  this heuristic is disabled. */
    342     Index acceptable_iter_max_;
     348    /** Indicates whether the algorithm should start directly with the
     349     *  restoratin phase */
     350    bool start_with_resto_;
    343351    //@}
    344352
     
    403411    Index count_successive_shortened_steps_;
    404412
    405     /** Counter for the number of successive iterations in which the
    406      *  nlp error was below the acceptable tolerance and a full step
    407      *  was accepted. */
    408     Index count_acceptable_iter_;
    409 
    410413    /** Flag indicating if a tiny step was detected in previous
    411414     *  iteration */
    412415    bool tiny_step_last_iteration_;
    413416
     417    /** @name Strategy objective that are used */
     418    //@{
    414419    SmartPtr<RestorationPhase> resto_phase_;
    415420    SmartPtr<PDSystemSolver> pd_solver_;
     421    SmartPtr<ConvergenceCheck> conv_check_;
     422    //@}
    416423  };
    417424
  • trunk/Algorithm/IpIpoptAlg.cpp

    r352 r424  
    99#include "IpIpoptAlg.hpp"
    1010#include "IpJournalist.hpp"
     11#include "IpRestoPhase.hpp"
    1112
    1213namespace Ipopt
     
    6667    bool retvalue = IpData().Initialize(Jnlst(),
    6768                                        options, prefix);
    68     ASSERT_EXCEPTION(retvalue, AlgorithmStrategyObject::FAILED_INITIALIZATION,
     69    ASSERT_EXCEPTION(retvalue, FAILED_INITIALIZATION,
    6970                     "the IpIpoptData object failed to initialize.");
    7071
     
    7273    retvalue = IpCq().Initialize(Jnlst(),
    7374                                 options, prefix);
    74     ASSERT_EXCEPTION(retvalue, AlgorithmStrategyObject::FAILED_INITIALIZATION,
     75    ASSERT_EXCEPTION(retvalue, FAILED_INITIALIZATION,
    7576                     "the IpIpoptCalculatedQuantities object failed to initialize.");
    7677
     
    7879    retvalue = IpNLP().Initialize(Jnlst(),
    7980                                  options, prefix);
    80     ASSERT_EXCEPTION(retvalue, AlgorithmStrategyObject::FAILED_INITIALIZATION,
     81    ASSERT_EXCEPTION(retvalue, FAILED_INITIALIZATION,
    8182                     "the IpIpoptNLP object failed to initialize.");
    8283
     
    124125  }
    125126
    126   IpoptAlgorithm::SolverReturn IpoptAlgorithm::Optimize()
     127  SolverReturn IpoptAlgorithm::Optimize()
    127128  {
    128129    DBG_START_METH("IpoptAlgorithm::Optimize", dbg_verbosity);
     
    171172        return SUCCESS;
    172173      }
     174      if (conv_status == ConvergenceCheck::CONVERGED_TO_ACCEPTABLE_POINT) {
     175        return STOP_AT_ACCEPTABLE_POINT;
     176      }
    173177      else if (conv_status == ConvergenceCheck::MAXITER_EXCEEDED) {
    174178        return MAXITER_EXCEEDED;
     
    177181    catch(TINY_STEP_DETECTED& exc) {
    178182      exc.ReportException(Jnlst());
    179 
    180183      return STOP_AT_TINY_STEP;
    181184    }
    182185    catch(ACCEPTABLE_POINT_REACHED& exc) {
    183186      exc.ReportException(Jnlst());
    184 
    185187      return STOP_AT_ACCEPTABLE_POINT;
    186188    }
    187 
    188     return FAILED;
     189    catch(LOCALLY_INFEASIBLE& exc) {
     190      exc.ReportException(Jnlst());
     191      return LOCAL_INFEASIBILITY;
     192    }
     193    catch(RESTORATION_FAILED& exc) {
     194      exc.ReportException(Jnlst());
     195      return RESTORATION_FAILURE;
     196    }
     197    catch(FEASIBILITY_PROBLEM_SOLVED& exc) {
     198      exc.ReportException(Jnlst());
     199      return SUCCESS;
     200    }
     201
     202    DBG_ASSERT(false && "Unknown return code in the algorithm");
    189203  }
    190204
  • trunk/Algorithm/IpIpoptAlg.hpp

    r352 r424  
    2222#include "IpIterationOutput.hpp"
    2323#include "IpIpoptType.hpp"
     24#include "IpAlgTypes.hpp"
    2425
    2526namespace Ipopt
     
    4344  {
    4445  public:
    45     /**@name Enumerations */
    46     //@{
    47     /** enum for the return from the optimize algorithm
    48      *  (obviously we need to add more) */
    49     enum SolverReturn {
    50       SUCCESS,
    51       MAXITER_EXCEEDED,
    52       FAILED,
    53       STOP_AT_TINY_STEP,
    54       STOP_AT_ACCEPTABLE_POINT
    55     };
    56     //@}
    5746
    5847    /**@name Constructors/Destructors */
  • trunk/Algorithm/IpIpoptCalculatedQuantities.cpp

    r361 r424  
    6565      trial_jac_d_cache_(1),
    6666      curr_jac_cT_times_vec_cache_(2),
     67      trial_jac_cT_times_vec_cache_(1),
    6768      curr_jac_dT_times_vec_cache_(2),
    68       trial_jac_cT_times_vec_cache_(1),
    6969      trial_jac_dT_times_vec_cache_(1),
    7070      curr_jac_c_times_vec_cache_(1),
    7171      curr_jac_d_times_vec_cache_(1),
    72       curr_exact_hessian_cache_(1),
    7372      curr_constraint_violation_cache_(2),
    7473      trial_constraint_violation_cache_(5),
     74      curr_nlp_constraint_violation_cache_(3),
     75      unscaled_curr_nlp_constraint_violation_cache_(3),
     76
     77      curr_exact_hessian_cache_(1),
    7578
    7679      curr_grad_lag_x_cache_(1),
     
    9699      curr_dual_infeasibility_cache_(3),
    97100      trial_dual_infeasibility_cache_(3),
     101      unscaled_curr_dual_infeasibility_cache_(3),
    98102      curr_complementarity_cache_(6),
    99103      trial_complementarity_cache_(6),
    100104      curr_centrality_measure_cache_(1),
    101105      curr_nlp_error_cache_(1),
     106      unscaled_curr_nlp_error_cache_(1),
    102107      curr_barrier_error_cache_(1),
    103108      curr_primal_dual_system_error_cache_(1),
     
    525530
    526531  Number
     532  IpoptCalculatedQuantities::unscaled_curr_f()
     533  {
     534    return ip_nlp_->NLP_scaling()->unapply_obj_scaling(curr_f());
     535  }
     536
     537  Number
    527538  IpoptCalculatedQuantities::trial_f()
    528539  {
     
    559570    DBG_PRINT((1,"result (trial_f) = %e\n", result));
    560571    return result;
     572  }
     573
     574  Number
     575  IpoptCalculatedQuantities::unscaled_trial_f()
     576  {
     577    return ip_nlp_->NLP_scaling()->unapply_obj_scaling(trial_f());
    561578  }
    562579
     
    774791    std::vector<Number> sdeps(1);
    775792    sdeps[0] = mu;
     793    DBG_PRINT((1,"curr_mu=%e\n",mu));
    776794
    777795    if (!curr_grad_barrier_obj_x_cache_.GetCachedResult(result, tdeps, sdeps)) {
     
    860878    std::vector<Number> sdeps(1);
    861879    sdeps[0] = mu;
     880    DBG_PRINT((1,"curr_mu=%e\n",mu));
    862881
    863882    if (!curr_grad_barrier_obj_s_cache_.GetCachedResult(result, tdeps, sdeps)) {
     
    10031022
    10041023  SmartPtr<const Vector>
     1024  IpoptCalculatedQuantities::unscaled_curr_c()
     1025  {
     1026    return ip_nlp_->NLP_scaling()->unapply_vector_scaling_c(curr_c());
     1027  }
     1028
     1029  SmartPtr<const Vector>
    10051030  IpoptCalculatedQuantities::trial_c()
    10061031  {
     
    10341059    }
    10351060    return result;
     1061  }
     1062
     1063  SmartPtr<const Vector>
     1064  IpoptCalculatedQuantities::unscaled_curr_d()
     1065  {
     1066    return ip_nlp_->NLP_scaling()->unapply_vector_scaling_d(curr_d());
    10361067  }
    10371068
     
    13271358                   dbg_verbosity);
    13281359    return trial_primal_infeasibility(constr_viol_normtype_);
     1360  }
     1361
     1362  Number
     1363  IpoptCalculatedQuantities::curr_nlp_constraint_violation
     1364  (ENormType NormType)
     1365  {
     1366    DBG_START_METH("IpoptCalculatedQuantities::curr_nlp_constraint_violation()",
     1367                   dbg_verbosity);
     1368    Number result;
     1369
     1370    SmartPtr<const Vector> x = ip_data_->curr()->x();
     1371
     1372    std::vector<const TaggedObject*> deps(1);
     1373    deps[0] = GetRawPtr(x);
     1374    std::vector<Number> sdeps(1);
     1375    sdeps[0] = (Number)NormType;
     1376
     1377    if (!curr_nlp_constraint_violation_cache_.GetCachedResult(result, deps, sdeps)) {
     1378      SmartPtr<const Vector> c = curr_c();
     1379      SmartPtr<const Vector> d = curr_d();
     1380
     1381      SmartPtr<Vector> d_viol_L = ip_nlp_->d_L()->MakeNewCopy();
     1382      ip_nlp_->Pd_L()->TransMultVector(-1., *d, 1., *d_viol_L);
     1383      SmartPtr<Vector> tmp = d_viol_L->MakeNew();
     1384      tmp->Set(0.);
     1385      d_viol_L->ElementWiseMax(*tmp);
     1386      DBG_PRINT_VECTOR(2, "d_viol_L", *d_viol_L);
     1387
     1388      SmartPtr<Vector> d_viol_U = ip_nlp_->d_U()->MakeNewCopy();
     1389      ip_nlp_->Pd_U()->TransMultVector(-1., *d, 1., *d_viol_U);
     1390      tmp = d_viol_U->MakeNew();
     1391      tmp->Set(0.);
     1392      d_viol_U->ElementWiseMin(*tmp);
     1393      DBG_PRINT_VECTOR(2, "d_viol_U", *d_viol_U);
     1394
     1395      std::vector<SmartPtr<const Vector> > vecs(3);
     1396      vecs[0] = c;
     1397      vecs[1] = GetRawPtr(d_viol_L);
     1398      vecs[2] = GetRawPtr(d_viol_U);
     1399      result = CalcNormOfType(NormType, vecs);
     1400      curr_nlp_constraint_violation_cache_.AddCachedResult(result, deps, sdeps);
     1401    }
     1402
     1403    return result;
     1404  }
     1405
     1406  Number
     1407  IpoptCalculatedQuantities::unscaled_curr_nlp_constraint_violation
     1408  (ENormType NormType)
     1409  {
     1410    DBG_START_METH("IpoptCalculatedQuantities::unscaled_curr_nlp_constraint_violation()",
     1411                   dbg_verbosity);
     1412    Number result;
     1413
     1414    SmartPtr<const Vector> x = ip_data_->curr()->x();
     1415
     1416    std::vector<const TaggedObject*> deps(1);
     1417    deps[0] = GetRawPtr(x);
     1418    std::vector<Number> sdeps(1);
     1419    sdeps[0] = (Number)NormType;
     1420
     1421    if (!unscaled_curr_nlp_constraint_violation_cache_.GetCachedResult(result, deps, sdeps)) {
     1422      SmartPtr<const Vector> c = unscaled_curr_c();
     1423
     1424      SmartPtr<const Vector> d = curr_d();
     1425
     1426      SmartPtr<Vector> d_viol_L = ip_nlp_->d_L()->MakeNewCopy();
     1427      ip_nlp_->Pd_L()->TransMultVector(-1., *d, 1., *d_viol_L);
     1428      SmartPtr<Vector> tmp = d_viol_L->MakeNew();
     1429      tmp->Set(0.);
     1430      d_viol_L->ElementWiseMax(*tmp);
     1431      DBG_PRINT_VECTOR(2, "d_viol_L", *d_viol_L);
     1432
     1433      SmartPtr<Vector> d_viol_U = ip_nlp_->d_U()->MakeNewCopy();
     1434      ip_nlp_->Pd_U()->TransMultVector(-1., *d, 1., *d_viol_U);
     1435      tmp = d_viol_U->MakeNew();
     1436      tmp->Set(0.);
     1437      d_viol_U->ElementWiseMin(*tmp);
     1438      DBG_PRINT_VECTOR(2, "d_viol_U", *d_viol_U);
     1439
     1440      std::vector<SmartPtr<const Vector> > vecs(3);
     1441      vecs[0] = c;
     1442      vecs[1] = GetRawPtr(d_viol_L);
     1443      vecs[2] = GetRawPtr(d_viol_U);
     1444      result = CalcNormOfType(NormType, vecs);
     1445      unscaled_curr_nlp_constraint_violation_cache_.AddCachedResult(result, deps, sdeps);
     1446    }
     1447
     1448    return result;
    13291449  }
    13301450
     
    19152035   const Vector& vec1, const Vector& vec2)
    19162036  {
    1917     std::vector<SmartPtr<const Vector> > vecs(2);
    1918     vecs[0] = &vec1;
    1919     vecs[1] = &vec2;
    1920 
    1921     return CalcNormOfType(NormType, vecs);
     2037    switch (NormType) {
     2038      case NORM_1 :
     2039      return vec1.Asum() + vec2.Asum();
     2040      case NORM_2 :
     2041      return sqrt(vec1.Nrm2()*vec2.Nrm2());
     2042      case NORM_MAX :
     2043      return Max(vec1.Amax(), vec2.Amax());
     2044      default:
     2045      DBG_ASSERT("Unknown NormType.");
     2046      return 0.0;
     2047    }
    19222048  }
    19232049
     
    21172243      }
    21182244      trial_dual_infeasibility_cache_.AddCachedResult(result, deps, sdeps);
     2245    }
     2246
     2247    return result;
     2248  }
     2249
     2250  Number
     2251  IpoptCalculatedQuantities::unscaled_curr_dual_infeasibility
     2252  (ENormType NormType)
     2253  {
     2254    DBG_START_METH("IpoptCalculatedQuantities::unscaled_curr_dual_infeasibility()",
     2255                   dbg_verbosity);
     2256    Number result;
     2257
     2258    SmartPtr<const Vector> x = ip_data_->curr()->x();
     2259    SmartPtr<const Vector> s = ip_data_->curr()->s();
     2260    SmartPtr<const Vector> y_c = ip_data_->curr()->y_c();
     2261    SmartPtr<const Vector> y_d = ip_data_->curr()->y_d();
     2262    SmartPtr<const Vector> z_L = ip_data_->curr()->z_L();
     2263    SmartPtr<const Vector> z_U = ip_data_->curr()->z_U();
     2264    SmartPtr<const Vector> v_L = ip_data_->curr()->v_L();
     2265    SmartPtr<const Vector> v_U = ip_data_->curr()->v_U();
     2266
     2267    std::vector<const TaggedObject*> deps(8);
     2268    deps[0] = GetRawPtr(x);
     2269    deps[1] = GetRawPtr(s);
     2270    deps[2] = GetRawPtr(y_c);
     2271    deps[3] = GetRawPtr(y_d);
     2272    deps[4] = GetRawPtr(z_L);
     2273    deps[5] = GetRawPtr(z_U);
     2274    deps[6] = GetRawPtr(v_L);
     2275    deps[7] = GetRawPtr(v_U);
     2276    std::vector<Number> sdeps(1);
     2277    sdeps[0] = (Number)NormType;
     2278
     2279    if (!unscaled_curr_dual_infeasibility_cache_.GetCachedResult(result, deps, sdeps)) {
     2280      SmartPtr<const Vector> grad_lag_x =
     2281        ip_nlp_->NLP_scaling()->unapply_grad_obj_scaling(curr_grad_lag_x());
     2282
     2283      Number obj_unscal = ip_nlp_->NLP_scaling()->unapply_obj_scaling(1.);
     2284      SmartPtr<const Vector> grad_lag_s;
     2285      if (obj_unscal != 1.) {
     2286        SmartPtr<Vector> tmp =
     2287          ip_nlp_->NLP_scaling()->apply_vector_scaling_d_NonConst(ConstPtr(curr_grad_lag_s()));
     2288        tmp->Scal(obj_unscal);
     2289        grad_lag_s = ConstPtr(tmp);
     2290      }
     2291      else {
     2292        grad_lag_s = ip_nlp_->NLP_scaling()->apply_vector_scaling_d(curr_grad_lag_s());
     2293      }
     2294
     2295      result = CalcNormOfType(NormType, *grad_lag_x, *grad_lag_s);
     2296      unscaled_curr_dual_infeasibility_cache_.AddCachedResult(result, deps, sdeps);
    21192297    }
    21202298
     
    22612439
    22622440  Number
     2441  IpoptCalculatedQuantities::unscaled_curr_complementarity
     2442  (Number mu, ENormType NormType)
     2443  {
     2444    return ip_nlp_->NLP_scaling()->unapply_obj_scaling(curr_complementarity(mu, NormType));
     2445  }
     2446
     2447  Number
    22632448  IpoptCalculatedQuantities::CalcCentralityMeasure(const Vector& compl_x_L,
    22642449      const Vector& compl_x_U,
     
    24132598      DBG_PRINT((1, "s_d = %lf, s_c = %lf\n", s_d, s_c));
    24142599
    2415       // Primal infeasibility
     2600      // Dual infeasibility
    24162601      DBG_PRINT((1, "curr_dual_infeasibility(NORM_MAX) = %8.2e\n",
    24172602                 curr_dual_infeasibility(NORM_MAX)));
    24182603      result = curr_dual_infeasibility(NORM_MAX)/s_d;
    2419       // Dual infeasibility
     2604      /*
     2605      // Primal infeasibility
    24202606      DBG_PRINT((1, "curr_primal_infeasibility(NORM_MAX) = %8.2e\n",
    24212607                 curr_primal_infeasibility(NORM_MAX)));
    24222608      result = Max(result, curr_primal_infeasibility(NORM_MAX));
     2609      */
     2610      result = Max(result, curr_nlp_constraint_violation(NORM_MAX));
    24232611      // Complementarity
    24242612      DBG_PRINT((1, "curr_complementarity(0., NORM_MAX) = %8.2e\n",
     
    24272615
    24282616      curr_nlp_error_cache_.AddCachedResult(result, tdeps);
     2617    }
     2618
     2619    return result;
     2620  }
     2621
     2622  Number
     2623  IpoptCalculatedQuantities::unscaled_curr_nlp_error()
     2624  {
     2625    DBG_START_METH("IpoptCalculatedQuantities::unscaled_curr_nlp_error()",
     2626                   dbg_verbosity);
     2627    DBG_ASSERT(initialize_called_);
     2628    Number result;
     2629
     2630    SmartPtr<const Vector> x = ip_data_->curr()->x();
     2631    SmartPtr<const Vector> s = ip_data_->curr()->s();
     2632    SmartPtr<const Vector> y_c = ip_data_->curr()->y_c();
     2633    SmartPtr<const Vector> y_d = ip_data_->curr()->y_d();
     2634    SmartPtr<const Vector> z_L = ip_data_->curr()->z_L();
     2635    SmartPtr<const Vector> z_U = ip_data_->curr()->z_U();
     2636    SmartPtr<const Vector> v_L = ip_data_->curr()->v_L();
     2637    SmartPtr<const Vector> v_U = ip_data_->curr()->v_U();
     2638
     2639    std::vector<const TaggedObject*> tdeps(8);
     2640    tdeps[0] = GetRawPtr(x);
     2641    tdeps[1] = GetRawPtr(s);
     2642    tdeps[2] = GetRawPtr(y_c);
     2643    tdeps[3] = GetRawPtr(y_d);
     2644    tdeps[4] = GetRawPtr(z_L);
     2645    tdeps[5] = GetRawPtr(z_U);
     2646    tdeps[6] = GetRawPtr(v_L);
     2647    tdeps[7] = GetRawPtr(v_U);
     2648
     2649    if (!unscaled_curr_nlp_error_cache_.GetCachedResult(result, tdeps)) {
     2650
     2651      // Dual infeasibility
     2652      result = unscaled_curr_dual_infeasibility(NORM_MAX);
     2653      // Constraint violation
     2654      result = Max(result, unscaled_curr_nlp_constraint_violation(NORM_MAX));
     2655      // Complementarity (ToDo use unscaled?)
     2656      DBG_PRINT((1, "curr_complementarity(0., NORM_MAX) = %8.2e\n",
     2657                 curr_complementarity(0., NORM_MAX)));
     2658      result = Max(result, unscaled_curr_complementarity(0., NORM_MAX));
     2659
     2660      unscaled_curr_nlp_error_cache_.AddCachedResult(result, tdeps);
    24292661    }
    24302662
     
    29533185  Number IpoptCalculatedQuantities::curr_gradBarrTDelta()
    29543186  {
     3187    DBG_START_METH("IpoptCalculatedQuantities::curr_gradBarrTDelta()",
     3188                   dbg_verbosity);
     3189
    29553190    Number result;
    29563191
     
    29643199    tdeps[2] = GetRawPtr(delta_x);
    29653200    tdeps[3] = GetRawPtr(delta_s);
    2966 
    2967     if (!curr_gradBarrTDelta_cache_.GetCachedResult(result, tdeps)) {
     3201    Number mu = ip_data_->curr_mu();
     3202    std::vector<Number> sdeps(1);
     3203    sdeps[0] = mu;
     3204    DBG_PRINT((1,"curr_mu=%e\n",mu));
     3205
     3206    if (!curr_gradBarrTDelta_cache_.GetCachedResult(result, tdeps, sdeps)) {
    29683207      result = curr_grad_barrier_obj_x()->Dot(*delta_x) +
    29693208               curr_grad_barrier_obj_s()->Dot(*delta_s);
    29703209
    2971       curr_gradBarrTDelta_cache_.AddCachedResult(result, tdeps);
    2972     }
    2973     return result;
     3210      curr_gradBarrTDelta_cache_.AddCachedResult(result, tdeps, sdeps);
     3211    }
     3212    return result;
     3213  }
     3214
     3215  bool IpoptCalculatedQuantities::IsSquareProblem() const
     3216  {
     3217    return (ip_data_->curr()->x()->Dim() == ip_data_->curr()->y_c()->Dim());
    29743218  }
    29753219
  • trunk/Algorithm/IpIpoptCalculatedQuantities.hpp

    r343 r424  
    7777    /** Value of objective function (at current point) */
    7878    Number curr_f();
     79    /** Unscaled value of the objective function (at the current point) */
     80    Number unscaled_curr_f();
    7981    /** Value of objective function (at trial point) */
    8082    Number trial_f();
     83    /** Unscaled value of the objective function (at the trial point) */
     84    Number unscaled_trial_f();
    8185    /** Gradient of objective function (at current point) */
    8286    SmartPtr<const Vector> curr_grad_f();
     
    115119    /** c(x) (at current point) */
    116120    SmartPtr<const Vector> curr_c();
     121    /** unscaled c(x) (at current point) */
     122    SmartPtr<const Vector> unscaled_curr_c();
    117123    /** c(x) (at trial point) */
    118124    SmartPtr<const Vector> trial_c();
    119125    /** d(x) (at current point) */
    120126    SmartPtr<const Vector> curr_d();
     127    /** unscaled d(x) (at current point) */
     128    SmartPtr<const Vector> unscaled_curr_d();
    121129    /** d(x) (at trial point) */
    122130    SmartPtr<const Vector> trial_d();
     
    171179     *  What type of norm is used depends on constr_viol_normtype */
    172180    Number trial_constraint_violation();
     181    /** Real constraint violation in a given norm (at current
     182     *  iterate).  This considers the inequality constraints without
     183     *  slacks. */
     184    Number curr_nlp_constraint_violation(ENormType NormType);
     185    /** Unscaled real constraint violation in a given norm (at current
     186     *  iterate).  This considers the inequality constraints without
     187     *  slacks. */
     188    Number unscaled_curr_nlp_constraint_violation(ENormType NormType);
    173189    //@}
    174190
     
    232248    /** Dual infeasibility in a given norm (at trial iterate) */
    233249    Number trial_dual_infeasibility(ENormType NormType);
     250    /** Unscaled dual infeasibility in a given norm (at current iterate) */
     251    Number unscaled_curr_dual_infeasibility(ENormType NormType);
    234252
    235253    /** Complementarity (for all complementarity conditions together)
     
    239257     *  in a given norm (at trial iterate) */
    240258    Number trial_complementarity(Number mu, ENormType NormType);
     259    /** Complementarity (for all complementarity conditions together)
     260     *  in a given norm (at current iterate) without NLP scaling. */
     261    Number unscaled_curr_complementarity(Number mu, ENormType NormType);
    241262
    242263    /** Centrality measure (in spirit of the -infinity-neighborhood. */
     
    248269    Number curr_centrality_measure();
    249270
    250     /** Scaled total optimality error for the original NLP at the
    251      *  current iterate. */
     271    /** Total optimality error for the original NLP at the current
     272     *  iterate, using scaling factors based on multipliers.  Note
     273     *  that here the constraint violation is measured without slacks
     274     *  (nlp_constraint_violation) */
    252275    Number curr_nlp_error();
    253 
    254     /** Scaled total optimality error for the barrier problem at the
    255      *  current iterate. */
     276    /** Total optimality error for the original NLP at the current
     277     *  iterate, but using no scaling based on multipliers, and no
     278     *  scaling for the NLP.  Note that here the constraint violation
     279     *  is measured without slacks (nlp_constraint_violation) */
     280    Number unscaled_curr_nlp_error();
     281
     282    /** Total optimality error for the barrier problem at the
     283     *  current iterate, using scaling factors based on multipliers. */
    256284    Number curr_barrier_error();
    257285
     
    336364    }
    337365
     366    /** Method returning true if this is a square problem */
     367    bool IsSquareProblem() const;
     368
    338369    /** Methods for IpoptType */
    339370    //@{
     
    439470    CachedResults<Number> curr_constraint_violation_cache_;
    440471    CachedResults<Number> trial_constraint_violation_cache_;
     472    CachedResults<Number> curr_nlp_constraint_violation_cache_;
     473    CachedResults<Number> unscaled_curr_nlp_constraint_violation_cache_;
    441474    //@}
    442475
     
    468501    CachedResults<Number> curr_dual_infeasibility_cache_;
    469502    CachedResults<Number> trial_dual_infeasibility_cache_;
     503    CachedResults<Number> unscaled_curr_dual_infeasibility_cache_;
    470504    CachedResults<Number> curr_complementarity_cache_;
    471505    CachedResults<Number> trial_complementarity_cache_;
    472506    CachedResults<Number> curr_centrality_measure_cache_;
    473507    CachedResults<Number> curr_nlp_error_cache_;
     508    CachedResults<Number> unscaled_curr_nlp_error_cache_;
    474509    CachedResults<Number> curr_barrier_error_cache_;
    475510    CachedResults<Number> curr_primal_dual_system_error_cache_;
  • trunk/Algorithm/IpIpoptData.cpp

    r343 r424  
    2626      initialize_called_(false),
    2727      have_prototypes_(false),
     28
    2829      free_mu_mode_(false),
    2930      tiny_step_flag_(false),
    3031
     32      info_regu_x_(0.),
    3133      info_alpha_primal_(0.),
    3234      info_alpha_primal_char_(' '),
    3335      info_alpha_dual_(0.),
    34       info_regu_x_(0.),
    3536      info_ls_count_(0),
    3637      info_skip_output_(false)
     
    5354      "implementation paper).  [Some other algorithmic features also use "
    5455      "this quantity.]");
    55     reg_options->AddLowerBoundedNumberOption(
    56       "dual_inf_tol",
    57       "Acceptance threshold for the unscaled dual infeasibility.",
    58       0.0, true, 1e-2,
    59       "Absolute tolerance on the dual infesaibility.  Successful termination "
    60       "requires that the (unscaled) dual infeasibility is less than this "
    61       "threshold.");
    62     reg_options->AddLowerBoundedNumberOption(
    63       "primal_inf_tol",
    64       "Acceptance threshold for the unscaled primal infeasibility.",
    65       0.0, true, 1e-2,
    66       "Absolute tolerance on the dual infesaibility.  Successful termination "
    67       "requires that the (unscaled) primal infeasibility is less than this "
    68       "threshold.");
    69     reg_options->AddLowerBoundedNumberOption(
    70       "compl_inf_tol",
    71       "Acceptance threshold for the complementarity conditions.",
    72       0.0, true, 1e-2,
    73       "Absolute tolerance on the complementarity.  Successful termination "
    74       "requires that the (unscaled) complementarity is less than this "
    75       "threshold.");
    7656  }
    7757
     
    8161  {
    8262    options.GetNumericValue("tol", tol_, prefix);
    83     options.GetNumericValue("dual_inf_tol", dual_inf_tol_, prefix);
    84     options.GetNumericValue("primal_inf_tol", primal_inf_tol_, prefix);
    85     options.GetNumericValue("compl_inf_tol", compl_inf_tol_, prefix);
    8663
    8764    iter_count_=0;
     
    10077      bool want_y_d,
    10178      bool want_z_L,
    102       bool want_z_U,
    103       bool want_v_L,
    104       bool want_v_U)
     79      bool want_z_U)
    10580  {
    10681    DBG_ASSERT(initialize_called_);
     
    126101                                  new_z_L, want_z_L,
    127102                                  new_z_U, want_z_U,
    128                                   new_v_L, want_v_L,
    129                                   new_v_U, want_v_U);
     103                                  new_v_L, new_v_U);
    130104    if (!retValue) {
    131105      return false;
     
    133107
    134108    new_s = new_y_d->MakeNew(); // same dimension as d
    135     //new_s_->Set(0.0);
    136109
    137110    iterates_space_ = new IteratesVectorSpace(*(new_x->OwnerSpace()), *(new_s->OwnerSpace()),
  • trunk/Algorithm/IpIpoptData.hpp

    r321 r424  
    5454                                  bool want_y_d,
    5555                                  bool want_z_L,
    56                                   bool want_z_U,
    57                                   bool want_v_L,
    58                                   bool want_v_U);
     56                                  bool want_z_U);
    5957
    6058    /** This method must be called to initialize the global
     
    295293      return tol_;
    296294    }
    297     Number dual_inf_tol() const
    298     {
    299       DBG_ASSERT(initialize_called_);
    300       return dual_inf_tol_;
    301     }
    302     Number primal_inf_tol() const
    303     {
    304       DBG_ASSERT(initialize_called_);
    305       return primal_inf_tol_;
    306     }
    307     Number compl_inf_tol() const
    308     {
    309       DBG_ASSERT(initialize_called_);
    310       return compl_inf_tol_;
    311     }
    312295    //@}
    313296
     
    459442    /** Overall convergence tolerance */
    460443    Number tol_;
    461     /** Tolerance on unscaled dual infeasibility */
    462     Number dual_inf_tol_;
    463     /** Tolerance on unscaled primal infeasibility */
    464     Number primal_inf_tol_;
    465     /** Tolerance on unscaled complementarity */
    466     Number compl_inf_tol_;
     444    //@}
     445
     446    /** @name Status data **/
     447    //@{
     448    /** flag indicating whether the algorithm is in the free mu mode */
     449    bool free_mu_mode_;
     450    /** flag indicating if a tiny step has been detected */
     451    bool tiny_step_flag_;
    467452    //@}
    468453
     
    486471    //@}
    487472
    488     /** @name Status data **/
    489     //@{
    490     /** flag indicating whether the algorithm is in the free mu mode */
    491     bool free_mu_mode_;
    492     /** flag indicating if a tiny step has been detected */
    493     bool tiny_step_flag_;
    494     //@}
    495 
    496473    /** VectorSpace for all the iterates */
    497474    SmartPtr<IteratesVectorSpace> iterates_space_;
  • trunk/Algorithm/IpIpoptNLP.hpp

    r353 r424  
    4646    virtual bool Initialize(const Journalist& jnlst,
    4747                            const OptionsList& options,
    48                             const std::string& prefix) = 0;
     48                            const std::string& prefix)
     49    {
     50      bool ret = true;
     51      if (IsValid(nlp_scaling_)) {
     52        ret = nlp_scaling_->Initialize(jnlst, options, prefix);
     53      }
     54      return ret;
     55    }
    4956
    5057    /**@name Possible Exceptions */
     
    6673                                      bool init_z_U,
    6774                                      SmartPtr<Vector>& v_L,
    68                                       bool init_v_L,
    69                                       SmartPtr<Vector>& v_U,
    70                                       bool init_v_U
     75                                      SmartPtr<Vector>& v_U
    7176                                     ) = 0;
    7277
     
    198203    //@}
    199204
    200   protected:
    201     /** Returns the scaling strategy object - may be NULL */
     205    /** solution routines */
     206    //@{
     207    virtual void FinalizeSolution(SolverReturn status,
     208                                  const Vector& x, const Vector& z_L, const Vector& z_U,
     209                                  const Vector& c, const Vector& d,
     210                                  const Vector& y_c, const Vector& y_d,
     211                                  Number obj_value)=0;
     212    //@}
     213
     214    /** Returns the scaling strategy object */
    202215    SmartPtr<NLPScalingObject> NLP_scaling() const
    203216    {
  • trunk/Algorithm/IpIteratesVector.cpp

    • Property svn:keywords set to Author Date Id Revision
    r321 r424  
    33// This code is published under the Common Public License.
    44//
    5 // $Id: IpCompoundVector.cpp,v 1.7 2005/03/21 16:08:47 andreasw Exp $
     5// $Id$
    66//
    77// Authors:  Carl Laird, Andreas Waechter     IBM    2004-08-13
  • trunk/Algorithm/IpIteratesVector.hpp

    • Property svn:keywords set to Author Date Id Revision
    r321 r424  
    33// This code is published under the Common Public License.
    44//
    5 // $Id: IpCompoundVector.hpp,v 1.7 2005/03/21 16:08:47 andreasw Exp $
     5// $Id$
    66//
    77// Authors:  Carl Laird, Andreas Waechter     IBM    2004-08-13
     
    1717  class IteratesVectorSpace;
    1818
    19   /** Specialized CompoundVector class specifically for the algorithm iterates.
    20    *  This class inherits from CompoundVector and is a specialized class for
    21    *  handling the iterates of the Ipopt Algorithm, that is,
    22    *  x, s, y_c, y_d, z_L, z_U, v_L, and v_U. It inherits from CompoundVector
    23    *  so it can behave like a CV in most calculations, but it has fixed
    24    *  dimensions and cannot be customized
     19  /** Specialized CompoundVector class specifically for the algorithm
     20   *  iterates.  This class inherits from CompoundVector and is a
     21   *  specialized class for handling the iterates of the Ipopt
     22   *  Algorithm, that is, x, s, y_c, y_d, z_L, z_U, v_L, and v_U. It
     23   *  inherits from CompoundVector so it can behave like a CV in most
     24   *  calculations, but it has fixed dimensions and cannot be
     25   *  customized
    2526   */
    2627  class IteratesVector : public CompoundVector
     
    3637    /** Make New methods */
    3738    //@{
    38     /** Use this method to create a new iterates vector. The MakeNew method
    39      *  on the Vector class also works, but it does not give the create_new
    40      *  option.
     39    /** Use this method to create a new iterates vector. The MakeNew
     40     *  method on the Vector class also works, but it does not give
     41     *  the create_new option.
    4142     */
    4243    SmartPtr<IteratesVector> MakeNewIteratesVector(bool create_new = true) const;
    4344
    44     /** Use this method to create a new iterates vector with a copy of all
    45      *  the data.
     45    /** Use this method to create a new iterates vector with a copy of
     46     *  all the data.
    4647     */
    4748    SmartPtr<IteratesVector> MakeNewIteratesVectorCopy() const
     
    5253    }
    5354
    54     /** Use this method to create a new iterates vector container. This creates a
    55      *  new NonConst container, but the elements inside the iterates vector may be
    56      *  const. Therefore, the container can be modified to point to new entries, but
    57      *  the existing entries may or may not be modifiable.
     55    /** Use this method to create a new iterates vector
     56     *  container. This creates a new NonConst container, but the
     57     *  elements inside the iterates vector may be const. Therefore,
     58     *  the container can be modified to point to new entries, but the
     59     *  existing entries may or may not be modifiable.
    5860     */
    5961    SmartPtr<IteratesVector> MakeNewContainer() const;
     
    6971
    7072    /** Get the x iterate (non-const) - this can only be called if the
    71      *  vector was created intenally, or the Set_x_NonConst method was used. */
     73     *  vector was created intenally, or the Set_x_NonConst method was
     74     *  used. */
    7275    SmartPtr<Vector> x_NonConst()
    7376    {
     
    7881    SmartPtr<Vector> create_new_x();
    7982
    80     /** Create a new vector in the x entry and copy the current values into it. */
     83    /** Create a new vector in the x entry and copy the current values
     84     *  into it. */
    8185    SmartPtr<Vector> create_new_x_copy()
    8286    {
     
    8791    }
    8892
    89     /** Set the x iterate (const). Sets the pointer, does NOT copy data. */
     93    /** Set the x iterate (const). Sets the pointer, does NOT copy
     94     *  data. */
    9095    void Set_x(const Vector& vec)
    9196    {
     
    9398    }
    9499
    95     /** Set the x iterate (non-const). Sets the pointer, does NOT copy data. */
     100    /** Set the x iterate (non-const). Sets the pointer, does NOT copy
     101     *  data. */
    96102    void Set_x_NonConst(Vector& vec)
    97103    {
     
    106112
    107113    /** Get the s iterate (non-const) - this can only be called if the
    108      *  vector was created intenally, or the Set_s_NonConst method was used. */
     114     *  vector was created intenally, or the Set_s_NonConst method was
     115     *  used. */
    109116    SmartPtr<Vector> s_NonConst()
    110117    {
     
    115122    SmartPtr<Vector> create_new_s();
    116123
    117     /** Create a new vector in the s entry and copy the current values into it. */
     124    /** Create a new vector in the s entry and copy the current values
     125     *  into it. */
    118126    SmartPtr<Vector> create_new_s_copy()
    119127    {
     
    124132    }
    125133
    126     /** Set the s iterate (const). Sets the pointer, does NOT copy data. */
     134    /** Set the s iterate (const). Sets the pointer, does NOT copy
     135     *  data. */
    127136    void Set_s(const Vector& vec)
    128137    {
     
    130139    }
    131140
    132     /** Set the s iterate (non-const). Sets the pointer, does NOT copy data. */
     141    /** Set the s iterate (non-const). Sets the pointer, does NOT copy
     142     *  data. */
    133143    void Set_s_NonConst(Vector& vec)
    134144    {
     
    142152    }
    143153
    144     /** Get the y_c iterate (non-const) - this can only be called if the
    145      *  vector was created intenally, or the Set_y_c_NonConst method was used. */
     154    /** Get the y_c iterate (non-const) - this can only be called if
     155     *  the vector was created intenally, or the Set_y_c_NonConst
     156     *  method was used. */
    146157    SmartPtr<Vector> y_c_NonConst()
    147158    {
     
    152163    SmartPtr<Vector> create_new_y_c();
    153164
    154     /** Create a new vector in the y_c entry and copy the current values into it. */
     165    /** Create a new vector in the y_c entry and copy the current
     166     *  values into it. */
    155167    SmartPtr<Vector> create_new_y_c_copy()
    156168    {
     
    161173    }
    162174
    163     /** Set the y_c iterate (const). Sets the pointer, does NOT copy data. */
     175    /** Set the y_c iterate (const). Sets the pointer, does NOT copy
     176     *  data. */
    164177    void Set_y_c(const Vector& vec)
    165178    {
     
    167180    }
    168181
    169     /** Set the y_c iterate (non-const). Sets the pointer, does NOT copy data. */
     182    /** Set the y_c iterate (non-const). Sets the pointer, does NOT
     183     *  copy data. */
    170184    void Set_y_c_NonConst(Vector& vec)
    171185    {
     
    179193    }
    180194
    181     /** Get the y_d iterate (non-const) - this can only be called if the
    182      *  vector was created intenally, or the Set_y_d_NonConst method was used. */
     195    /** Get the y_d iterate (non-const) - this can only be called if
     196     *  the vector was created intenally, or the Set_y_d_NonConst
     197     *  method was used. */
    183198    SmartPtr<Vector> y_d_NonConst()
    184199    {
     
    189204    SmartPtr<Vector> create_new_y_d();
    190205
    191     /** Create a new vector in the y_d entry and copy the current values into it. */
     206    /** Create a new vector in the y_d entry and copy the current
     207     *  values into it. */
    192208    SmartPtr<Vector> create_new_y_d_copy()
    193209    {
     
    198214    }
    199215
    200     /** Set the y_d iterate (const). Sets the pointer, does NOT copy data. */
     216    /** Set the y_d iterate (const). Sets the pointer, does NOT copy
     217     *  data. */
    201218    void Set_y_d(const Vector& vec)
    202219    {
     
    204221    }
    205222
    206     /** Set the y_d iterate (non-const). Sets the pointer, does NOT copy data. */
     223    /** Set the y_d iterate (non-const). Sets the pointer, does NOT
     224     *  copy data. */
    207225    void Set_y_d_NonConst(Vector& vec)
    208226    {
     
    216234    }
    217235
    218     /** Get the z_L iterate (non-const) - this can only be called if the
    219      *  vector was created intenally, or the Set_z_L_NonConst method was used. */
     236    /** Get the z_L iterate (non-const) - this can only be called if
     237     *  the vector was created intenally, or the Set_z_L_NonConst
     238     *  method was used. */
    220239    SmartPtr<Vector> z_L_NonConst()
    221240    {
     
    226245    SmartPtr<Vector> create_new_z_L() ;
    227246
    228     /** Create a new vector in the z_L entry and copy the current values into it. */
     247    /** Create a new vector in the z_L entry and copy the current
     248     *  values into it. */
    229249    SmartPtr<Vector> create_new_z_L_copy()
    230250    {
     
    235255    }
    236256
    237     /** Set the z_L iterate (const). Sets the pointer, does NOT copy data. */
     257    /** Set the z_L iterate (const). Sets the pointer, does NOT copy
     258     *  data. */
    238259    void Set_z_L(const Vector& vec)
    239260    {
     
    241262    }
    242263
    243     /** Set the z_L iterate (non-const). Sets the pointer, does NOT copy data. */
     264    /** Set the z_L iterate (non-const). Sets the pointer, does NOT
     265     *  copy data. */
    244266    void Set_z_L_NonConst(Vector& vec)
    245267    {
     
    253275    }
    254276
    255     /** Get the z_U iterate (non-const) - this can only be called if the
    256      *  vector was created intenally, or the Set_z_U_NonConst method was used. */
     277    /** Get the z_U iterate (non-const) - this can only be called if
     278     *  the vector was created intenally, or the Set_z_U_NonConst
     279     *  method was used. */
    257280    SmartPtr<Vector> z_U_NonConst()
    258281    {
     
    263286    SmartPtr<Vector> create_new_z_U();
    264287
    265     /** Create a new vector in the z_U entry and copy the current values into it. */
     288    /** Create a new vector in the z_U entry and copy the current
     289     *  values into it. */
    266290    SmartPtr<Vector> create_new_z_U_copy()
    267291    {
     
    272296    }
    273297
    274     /** Set the z_U iterate (const). Sets the pointer, does NOT copy data. */
     298    /** Set the z_U iterate (const). Sets the pointer, does NOT copy
     299     *  data. */
    275300    void Set_z_U(const Vector& vec)
    276301    {
     
    278303    }
    279304
    280     /** Set the z_U iterate (non-const). Sets the pointer, does NOT copy data. */
     305    /** Set the z_U iterate (non-const). Sets the pointer, does NOT
     306     *  copy data. */
    281307    void Set_z_U_NonConst(Vector& vec)
    282308    {
     
    290316    }
    291317
    292     /** Get the v_L iterate (non-const) - this can only be called if the
    293      *  vector was created intenally, or the Set_v_L_NonConst method was used. */
     318    /** Get the v_L iterate (non-const) - this can only be called if
     319     *  the vector was created intenally, or the Set_v_L_NonConst
     320     *  method was used. */
    294321    SmartPtr<Vector> v_L_NonConst()
    295322    {
     
    300327    SmartPtr<Vector> create_new_v_L();
    301328
    302     /** Create a new vector in the v_L entry and copy the current values into it. */
     329    /** Create a new vector in the v_L entry and copy the current
     330     *  values into it. */
    303331    SmartPtr<Vector> create_new_v_L_copy()
    304332    {
     
    309337    }
    310338
    311     /** Set the v_L iterate (const). Sets the pointer, does NOT copy data. */
     339    /** Set the v_L iterate (const). Sets the pointer, does NOT copy
     340     *  data. */
    312341    void Set_v_L(const Vector& vec)
    313342    {
     
    315344    }
    316345
    317     /** Set the v_L iterate (non-const). Sets the pointer, does NOT copy data. */
     346    /** Set the v_L iterate (non-const). Sets the pointer, does NOT
     347     *  copy data. */
    318348    void Set_v_L_NonConst(Vector& vec)
    319349    {
     
    327357    }
    328358
    329     /** Get the v_U iterate (non-const) - this can only be called if the
    330      *  vector was created intenally, or the Set_v_U_NonConst method was used. */
     359    /** Get the v_U iterate (non-const) - this can only be called if
     360     *  the vector was created intenally, or the Set_v_U_NonConst
     361     *  method was used. */
    331362    SmartPtr<Vector> v_U_NonConst()
    332363    {
     
    337368    SmartPtr<Vector> create_new_v_U();
    338369
    339     /** Create a new vector in the v_U entry and copy the current values into it. */
     370    /** Create a new vector in the v_U entry and copy the current
     371     *  values into it. */
    340372    SmartPtr<Vector> create_new_v_U_copy()
    341373    {
     
    346378    }
    347379
    348     /** Set the v_U iterate (const). Sets the pointer, does NOT copy data. */
     380    /** Set the v_U iterate (const). Sets the pointer, does NOT copy
     381     *  data. */
    349382    void Set_v_U(const Vector& vec)
    350383    {
     
    352385    }
    353386
    354     /** Set the v_U iterate (non-const). Sets the pointer, does NOT copy data. */
     387    /** Set the v_U iterate (non-const). Sets the pointer, does NOT
     388     *  copy data. */
    355389    void Set_v_U_NonConst(Vector& vec)
    356390    {
     
    358392    }
    359393
    360     /** Set the primal variables all in one shot. Sets the pointers, does NOT copy data */
     394    /** Set the primal variables all in one shot. Sets the pointers,
     395     *  does NOT copy data */
    361396    void Set_primal(const Vector& x, const Vector& s)
    362397    {
     
    370405    }
    371406
    372     /** Set the eq multipliers all in one shot. Sets the pointers, does not copy data. */
     407    /** Set the eq multipliers all in one shot. Sets the pointers,
     408     *  does not copy data. */
    373409    void Set_eq_mult(const Vector& y_c, const Vector& y_d)
    374410    {
     
    382418    }
    383419
    384     /** Set the bound multipliers all in one shot. Sets the pointers, does not copy data. */
     420    /** Set the bound multipliers all in one shot. Sets the pointers,
     421     *  does not copy data. */
    385422    void Set_bound_mult(const Vector& z_L, const Vector& z_U, const Vector& v_L, const Vector& v_U)
    386423    {
     
    398435    }
    399436
    400     /** Get a sum of the tags of the contained items. There is no guarantee that
    401      *  this is unique, but there is a high chance it is unique and it can
    402      *  be used for debug checks relatively reliably.
     437    /** Get a sum of the tags of the contained items. There is no
     438     *  guarantee that this is unique, but there is a high chance it
     439     *  is unique and it can be used for debug checks relatively
     440     *  reliably.
    403441     */
    404442    TaggedObject::Tag GetTagSum() const
     
    436474
    437475  private:
    438     /**@name Default Compiler Generated Methods
    439      * (Hidden to avoid implicit creation/calling).
    440      * These methods are not implemented and
    441      * we do not want the compiler to implement
    442      * them for us, so we declare them private
    443      * and do not define them. This ensures that
     476    /**@name Default Compiler Generated Methods (Hidden to avoid
     477     * implicit creation/calling).  These methods are not implemented
     478     * and we do not want the compiler to implement them for us, so we
     479     * declare them private and do not define them. This ensures that
    444480     * they will not be implicitly created/called.
    445481     */
     
    457493    const IteratesVectorSpace* owner_space_;
    458494
    459     /** private method to return the const element from the compound vector.
    460      *  This method will return NULL if none is currently set.
     495    /** private method to return the const element from the compound
     496     *  vector.  This method will return NULL if none is currently
     497     *  set.
    461498     */
    462499    SmartPtr<const Vector> GetIterateFromComp(Index i) const
     
    468505    }
    469506
    470     /** private method to return the non-const element from the compound vector.
    471      *  This method will return NULL if none is currently set.
     507    /** private method to return the non-const element from the
     508     *  compound vector.  This method will return NULL if none is
     509     *  currently set.
    472510     */
    473511    SmartPtr<Vector> GetNonConstIterateFromComp(Index i)
     
    481519  };
    482520
    483   /** Vector Space for the IteratesVector class.
    484    *  This is a specialized vector space for the IteratesVector class.
     521  /** Vector Space for the IteratesVector class.  This is a
     522   *  specialized vector space for the IteratesVector class.
    485523   */
    486524  class IteratesVectorSpace : public CompoundVectorSpace
     
    503541    /** Method for creating vectors . */
    504542    //@{
    505     /** Use this to create a new IteratesVector. You can pass-in create_new = false
    506      *  if you only want a container and do not want vectors allocated.
     543    /** Use this to create a new IteratesVector. You can pass-in
     544     *  create_new = false if you only want a container and do not
     545     *  want vectors allocated.
    507546     */
    508547    virtual IteratesVector* MakeNewIteratesVector(bool create_new = true) const
     
    532571
    533572
    534     /** This method overloads ComooundVectorSpace::MakeNewCompoundVector to make
    535      *  sure that we get a vector of the correct type
     573    /** This method overloads
     574     *  ComooundVectorSpace::MakeNewCompoundVector to make sure that
     575     *  we get a vector of the correct type
    536576     */
    537577    virtual CompoundVector* MakeNewCompoundVector(bool create_new = true) const
     
    540580    }
    541581
    542     /** This method creates a new vector (and allocates space in all the
    543      *  contained vectors. This is really only used for code that does not
    544      *  know what type of vector it is dealing with - for example, this method
    545      *  is called from Vector::MakeNew()
     582    /** This method creates a new vector (and allocates space in all
     583     *  the contained vectors. This is really only used for code that
     584     *  does not know what type of vector it is dealing with - for
     585     *  example, this method is called from Vector::MakeNew()
    546586     */
    547587    virtual Vector* MakeNew() const
     
    551591    //@}
    552592
    553     /** This method hides the CompoundVectorSpace::SetCompSpace method since
    554      *  the components of the Iterates are fixed at construction.
     593    /** This method hides the CompoundVectorSpace::SetCompSpace method
     594     *  since the components of the Iterates are fixed at
     595     *  construction.
    555596     */
    556597    virtual void SetCompSpace(Index icomp, const VectorSpace& vec_space)
     
    560601
    561602  private:
    562     /**@name Default Compiler Generated Methods
    563     * (Hidden to avoid implicit creation/calling).
    564     * These methods are not implemented and
    565     * we do not want the compiler to implement
    566     * them for us, so we declare them private
    567     * and do not define them. This ensures that
     603    /**@name Default Compiler Generated Methods (Hidden to avoid
     604    * implicit creation/calling).  These methods are not implemented
     605    * and we do not want the compiler to implement them for us, so we
     606    * declare them private and do not define them. This ensures that
    568607    * they will not be implicitly created/called. */
    569608    //@{
  • trunk/Algorithm/IpMonotoneMuUpdate.cpp

    r321 r424  
    4242  void MonotoneMuUpdate::RegisterOptions(SmartPtr<RegisteredOptions> roptions)
    4343  {
    44     roptions->AddLowerBoundedNumberOption("mu0", "initial value for the barrier parameter, mu",
    45                                           0.0, true, 0.1);
    46     //     roptions->AddLowerBoundedNumberOption("kappa_epsilon", "???",
    47     //                                    0.0, true, 10.0);
    48     //     roptions->AddBoundedNumberOption("kappa_mu", "???",
    49     //                               0.0, true, 1.0, true, 0.2);
    50     //     roptions->AddBoundedNumberOption("theta_mu", "???",
    51     //                               1.0, true, 2.0, true, 1.5);
    52     //     roptions->AddBoundedNumberOption("tau_min", "???",
    53     //                                       0.0, true, 1.0, true, 0.99);
     44    roptions->AddLowerBoundedNumberOption(
     45      "mu_init", "Initial value for the barrier parameter.",
     46      0.0, true,
     47      0.1,
     48      "This option determines the initial value for the barrier parameter "
     49      "(mu).  It is only relevant in the monotone, Fiacco-McCormick "
     50      "version of the algorithm., i.e., if \"mu_strategy\" is chosen "
     51      "as \"monotone\"");
     52    roptions->AddLowerBoundedNumberOption(
     53      "barrier_tol_factor",
     54      "Factor for mu in barrier stop test.",
     55      0.0, true,
     56      10.0,
     57      "The convergence tolerance for each barrier problem in the montone mode "
     58      "is the value of the barrier parameter times this value. This option is "
     59      "also used in class \"AdaptiveMuUpdate\" during the monotone mode. "
     60      "(This is kappa_epsilon in implementation paper).");
     61    roptions->AddBoundedNumberOption(
     62      "mu_linear_decrease_factor",
     63      "Determines linear decrease rate of barrier parameter.",
     64      0.0, true, 1.0, true,
     65      0.2,
     66      "For the Fiacco-McCormick update procedure the new barrier parameter mu "
     67      "is obtained by taking the minimum of mu*\"mu_linear_decrease_factor\" "
     68      "and mu^\"superlinear_decrease_power\".  (This is kappa_mu in "
     69      "implementation paper.) [This option is also used in "
     70      "AdaptiveMuUpdate]");
     71    roptions->AddBoundedNumberOption(
     72      "mu_superlinear_decrease_power",
     73      "Determines superlinear decrease rate of barrier parameter.",
     74      1.0, true, 2.0, true,
     75      1.5,
     76      "For the Fiacco-McCormick update procedure the new barrier parameter mu "
     77      "is obtained by taking the minimum of mu*\"mu_linear_decrease_factor\" "
     78      "and mu^\"superlinear_decrease_power\".  (This is theta_mu in "
     79      "implementation paper.) [This option is also used in "
     80      "AdaptiveMuUpdate]");
     81    roptions->AddBoundedNumberOption(
     82      "tau_min",
     83      "Lower bound on fraction-to-the-boundary parameter tau.",
     84      0.0, true, 1.0, true,
     85      0.99,
     86      "(This is tau_min in implementation paper.) [This option is also "
     87      "used in AdaptiveMuUpdate]");
    5488  }
    5589
     
    5791                                        const std::string& prefix)
    5892  {
    59     options.GetNumericValue("mu0", mu0_, prefix);
    60     options.GetNumericValue("kappa_epsilon", kappa_epsilon_, prefix);
    61     options.GetNumericValue("kappa_mu", kappa_mu_, prefix);
    62     options.GetNumericValue("theta_mu", theta_mu_, prefix);
     93    options.GetNumericValue("mu_init", mu_init_, prefix);
     94    options.GetNumericValue("barrier_tol_factor", barrier_tol_factor_, prefix);
     95    options.GetNumericValue("mu_linear_decrease_factor", mu_linear_decrease_factor_, prefix);
     96    options.GetNumericValue("mu_superlinear_decrease_power", mu_superlinear_decrease_power_, prefix);
    6397    options.GetNumericValue("tau_min", tau_min_, prefix);
     98    options.GetNumericValue("compl_inf_tol", compl_inf_tol_, prefix);
    6499
    65     IpData().Set_mu(mu0_);
    66     Number tau = Max(tau_min_, 1.0 - mu0_);
     100    IpData().Set_mu(mu_init_);
     101    Number tau = Max(tau_min_, 1.0 - mu_init_);
    67102    IpData().Set_tau(tau);
    68103
     
    90125                   "Optimality Error for Barrier Sub-problem = %e\n",
    91126                   sub_problem_error);
    92     Number kappa_eps_mu = kappa_epsilon_ * mu;
     127    Number kappa_eps_mu = barrier_tol_factor_ * mu;
    93128
    94129    bool done = false;
     
    126161      else {
    127162        sub_problem_error = IpCq().curr_barrier_error();
    128         kappa_eps_mu = kappa_epsilon_ * mu;
     163        kappa_eps_mu = barrier_tol_factor_ * mu;
    129164        done = (sub_problem_error > kappa_eps_mu);
    130165      }
     
    148183    Number mu = IpData().curr_mu();
    149184    Number tol = IpData().tol();
    150     Number compl_inf_tol = IpData().compl_inf_tol();
    151185
    152     new_mu = Min( kappa_mu_*mu, pow(mu, theta_mu_) );
    153     new_mu = Max(new_mu, Min(tol, compl_inf_tol)/10.);
     186    // Here we need the complementarity tolerance that is posed to the
     187    // scaled problem
     188    Number compl_inf_tol =
     189      IpNLP().NLP_scaling()->apply_obj_scaling(compl_inf_tol_);
     190
     191    new_mu = Min( mu_linear_decrease_factor_*mu,
     192                  pow(mu, mu_superlinear_decrease_power_) );
     193    new_mu = Max(new_mu, Min(tol, compl_inf_tol)/(barrier_tol_factor_+1.));
    154194
    155195    // update the fraction to the boundary parameter
  • trunk/Algorithm/IpMonotoneMuUpdate.hpp

    r321 r424  
    7474    /** @name Algorithmic parameters */
    7575    //@{
    76     /** Initial value of the barrier parameter (TODO: NEED GOOD NAMES HERE)*/
    77     Number mu0_;
    78     Number kappa_epsilon_;
    79     Number kappa_mu_;
    80     Number theta_mu_;
     76    /** Initial value of the barrier parameter */
     77    Number mu_init_;
     78    Number barrier_tol_factor_;
     79    Number mu_linear_decrease_factor_;
     80    Number mu_superlinear_decrease_power_;
    8181    /** Tau_min for fraction to boundary rule */
    8282    Number tau_min_;
     83    Number compl_inf_tol_;
    8384    //@}
    8485
  • trunk/Algorithm/IpMuUpdate.hpp

    r310 r424  
    1212#include "IpUtils.hpp"
    1313#include "IpAlgStrategy.hpp"
    14 #include "IpIpoptType.hpp"
    1514
    1615namespace Ipopt
    1716{
    18 
    19   DeclareIpoptType(MuUpdate);
    20 
    2117  /** Abstract Base Class for classes that implement methods for computing
    2218   *  the barrier and fraction-to-the-boundary rule parameter for the
     
    2925    //@{
    3026    /** Default Constructor */
    31     MuUpdate();
     27    MuUpdate()
     28    {}
    3229
    3330    /** Default destructor */
     
    4643     *  for setting the fraction-to-the-boundary parameter tau */
    4744    virtual void UpdateBarrierParameter() = 0;
    48 
    49     /** Methods for IpoptType */
    50     //@{
    51     static void RegisterOptions(SmartPtr<RegisteredOptions> roptions);
    52     //@}
    5345
    5446  private:
  • trunk/Algorithm/IpNLPScaling.cpp

    • Property svn:keywords set to Author Date Id Revision
    r340 r424  
    33// This code is published under the Common Public License.
    44//
    5 // $Id: IpOrigIpoptNLP.cpp 321 2005-06-20 21:53:55Z andreasw $
     5// $Id$
    66//
    77// Authors:  Carl Laird, Andreas Waechter     IBM    2004-08-13
     
    1212{
    1313
    14   SmartPtr<Vector> NLPScalingObject::apply_vector_scaling_x_L_NonConst(SmartPtr<Matrix> Px_L, const SmartPtr<const Vector>& l, const SmartPtr<const VectorSpace> x_space)
    15   {
    16     SmartPtr<Vector> tmp_x = x_space->MakeNew();
    17 
    18     // move to full x space
    19     Px_L->MultVector(1.0, *l, 0.0, *tmp_x);
    20 
    21     // scale in full x space
    22     tmp_x = apply_vector_scaling_x_NonConst(ConstPtr(tmp_x));
    23 
    24     // move back to x_L space
    25     SmartPtr<Vector> scaled_x_L = l->MakeNew();
    26     Px_L->TransMultVector(1.0, *tmp_x, 0.0, *scaled_x_L);
    27 
    28     return scaled_x_L;
    29   }
    30 
    31   SmartPtr<Vector> NLPScalingObject::apply_vector_scaling_x_U_NonConst(SmartPtr<Matrix> Px_U, const SmartPtr<const Vector>& u, const SmartPtr<const VectorSpace> x_space)
    32   {
    33     SmartPtr<Vector> tmp_x = x_space->MakeNew();
    34 
    35     // move to full x space
    36     Px_U->MultVector(1.0, *u, 0.0, *tmp_x);
    37 
    38     // scale in full x space
    39     tmp_x = apply_vector_scaling_x_NonConst(ConstPtr(tmp_x));
    40 
    41     // move back to x_L space
    42     SmartPtr<Vector> scaled_x_U = u->MakeNew();
    43     Px_U->TransMultVector(1.0, *tmp_x, 0.0, *scaled_x_U);
    44 
    45     return scaled_x_U;
    46   }
    47 
    48   SmartPtr<Vector> NLPScalingObject::apply_vector_scaling_d_L_NonConst(SmartPtr<Matrix> Pd_L, const SmartPtr<const Vector>& l, const SmartPtr<const VectorSpace> d_space)
    49   {
    50     SmartPtr<Vector> tmp_d = d_space->MakeNew();
    51 
    52     // move to full d space
    53     Pd_L->MultVector(1.0, *l, 0.0, *tmp_d);
    54 
    55     // scale in full d space
    56     tmp_d = apply_vector_scaling_d_NonConst(ConstPtr(tmp_d));
    57 
    58     // move back to d_L space
    59     SmartPtr<Vector> scaled_d_L = l->MakeNew();
    60     Pd_L->TransMultVector(1.0, *tmp_d, 0.0, *scaled_d_L);
    61 
    62     return scaled_d_L;
    63   }
    64 
    65   SmartPtr<Vector> NLPScalingObject::apply_vector_scaling_d_U_NonConst(SmartPtr<Matrix> Pd_U, const SmartPtr<const Vector>& u, const SmartPtr<const VectorSpace> d_space)
    66   {
    67     SmartPtr<Vector> tmp_d = d_space->MakeNew();
    68 
    69     // move to full d space
    70     Pd_U->MultVector(1.0, *u, 0.0, *tmp_d);
    71 
    72     // scale in full d space
    73     tmp_d = apply_vector_scaling_d_NonConst(ConstPtr(tmp_d));
    74 
    75     // move back to d_L space
    76     SmartPtr<Vector> scaled_d_U = u->MakeNew();
    77     Pd_U->TransMultVector(1.0, *tmp_d, 0.0, *scaled_d_U);
    78 
    79     return scaled_d_U;
    80   }
    81 
    82   SmartPtr<Vector> NLPScalingObject::apply_grad_obj_scaling_NonConst(const SmartPtr<const Vector>& v)
     14  DBG_SET_VERBOSITY(0);
     15
     16  SmartPtr<Vector> NLPScalingObject::apply_vector_scaling_x_LU_NonConst(
     17    const Matrix& Px_LU,
     18    const SmartPtr<const Vector>& lu,
     19    const VectorSpace& x_space)
     20  {
     21    SmartPtr<Vector> scaled_x_LU = lu->MakeNew();
     22    if (have_x_scaling()) {
     23      SmartPtr<Vector> tmp_x = x_space.MakeNew();
     24
     25      // move to full x space
     26      Px_LU.MultVector(1.0, *lu, 0.0, *tmp_x);
     27
     28      // scale in full x space
     29      tmp_x = apply_vector_scaling_x_NonConst(ConstPtr(tmp_x));
     30
     31      // move back to x_L space
     32      Px_LU.TransMultVector(1.0, *tmp_x, 0.0, *scaled_x_LU);
     33    }
     34    else {
     35      scaled_x_LU->Copy(*lu);
     36    }
     37
     38    return scaled_x_LU;
     39  }
     40
     41  SmartPtr<const Vector> NLPScalingObject::apply_vector_scaling_x_LU(
     42    const Matrix& Px_LU,
     43    const SmartPtr<const Vector>& lu,
     44    const VectorSpace& x_space)
     45  {
     46    if (have_x_scaling()) {
     47      return ConstPtr(apply_vector_scaling_x_LU_NonConst(Px_LU, lu, x_space));
     48    }
     49    else {
     50      return lu;
     51    }
     52  }
     53
     54  SmartPtr<Vector> NLPScalingObject::apply_vector_scaling_d_LU_NonConst(
     55    const Matrix& Pd_LU,
     56    const SmartPtr<const Vector>& lu,
     57    const VectorSpace& d_space)
     58  {
     59    SmartPtr<Vector> scaled_d_LU = lu->MakeNew();
     60    if (have_d_scaling()) {
     61      SmartPtr<Vector> tmp_d = d_space.MakeNew();
     62
     63      // move to full d space
     64      Pd_LU.MultVector(1.0, *lu, 0.0, *tmp_d);
     65
     66      // scale in full x space
     67      tmp_d = apply_vector_scaling_d_NonConst(ConstPtr(tmp_d));
     68
     69      // move back to x_L space
     70      Pd_LU.TransMultVector(1.0, *tmp_d, 0.0, *scaled_d_LU);
     71    }
     72    else {
     73      scaled_d_LU->Copy(*lu);
     74    }
     75
     76    return scaled_d_LU;
     77  }
     78
     79  SmartPtr<const Vector> NLPScalingObject::apply_vector_scaling_d_LU(
     80    const Matrix& Pd_LU,
     81    const SmartPtr<const Vector>& lu,
     82    const VectorSpace& d_space)
     83  {
     84    if (have_d_scaling()) {
     85      return ConstPtr(apply_vector_scaling_d_LU_NonConst(Pd_LU, lu, d_space));
     86    }
     87    else {
     88      return lu;
     89    }
     90  }
     91
     92  SmartPtr<Vector> NLPScalingObject::unapply_vector_scaling_d_LU_NonConst(
     93    const Matrix& Pd_LU,
     94    const SmartPtr<const Vector>& lu,
     95    const VectorSpace& d_space)
     96  {
     97    SmartPtr<Vector> unscaled_d_LU = lu->MakeNew();
     98    if (have_d_scaling()) {
     99      SmartPtr<Vector> tmp_d = d_space.MakeNew();
     100
     101      // move to full d space
     102      Pd_LU.MultVector(1.0, *lu, 0.0, *tmp_d);
     103
     104      // scale in full x space
     105      tmp_d = unapply_vector_scaling_d_NonConst(ConstPtr(tmp_d));
     106
     107      // move back to x_L space
     108      Pd_LU.TransMultVector(1.0, *tmp_d, 0.0, *unscaled_d_LU);
     109    }
     110    else {
     111      unscaled_d_LU->Copy(*lu);
     112    }
     113
     114    return unscaled_d_LU;
     115  }
     116
     117  SmartPtr<const Vector> NLPScalingObject::unapply_vector_scaling_d_LU(
     118    const Matrix& Pd_LU,
     119    const SmartPtr<const Vector>& lu,
     120    const VectorSpace& d_space)
     121  {
     122    if (have_d_scaling()) {
     123      return ConstPtr(unapply_vector_scaling_d_LU_NonConst(Pd_LU, lu, d_space));
     124    }
     125    else {
     126      return lu;
     127    }
     128  }
     129
     130  SmartPtr<Vector> NLPScalingObject::apply_grad_obj_scaling_NonConst(
     131    const SmartPtr<const Vector>& v)
    83132  {
    84133    SmartPtr<Vector> scaled_v = unapply_vector_scaling_x_NonConst(v);
    85134    Number df = apply_obj_scaling(1.0);
    86     scaled_v->Scal(df);
     135    if (df != 1.) {
     136      scaled_v->Scal(df);
     137    }
    87138    return scaled_v;
    88139  }
    89140
    90   SmartPtr<const Vector> NLPScalingObject::apply_grad_obj_scaling(const SmartPtr<const Vector>& v)
    91   {
    92     SmartPtr<Vector> scaled_v = apply_grad_obj_scaling_NonConst(v);
    93     return ConstPtr(scaled_v);
    94   }
    95 
    96   SmartPtr<Vector> NLPScalingObject::unapply_grad_obj_scaling_NonConst(const SmartPtr<const Vector>& v)
     141  SmartPtr<const Vector> NLPScalingObject::apply_grad_obj_scaling(
     142    const SmartPtr<const Vector>& v)
     143  {
     144    Number df = apply_obj_scaling(1.);
     145    if (df != 1.) {
     146      SmartPtr<Vector> scaled_v = apply_grad_obj_scaling_NonConst(v);
     147      return ConstPtr(scaled_v);
     148    }
     149    else {
     150      SmartPtr<const Vector> scaled_v = unapply_vector_scaling_x(v);
     151      return scaled_v;
     152    }
     153  }
     154
     155  SmartPtr<Vector> NLPScalingObject::unapply_grad_obj_scaling_NonConst(
     156    const SmartPtr<const Vector>& v)
    97157  {
    98158    SmartPtr<Vector> unscaled_v = apply_vector_scaling_x_NonConst(v);
    99     Number df = unapply_obj_scaling(1.0);
    100     unscaled_v->Scal(df);
     159    Number df = unapply_obj_scaling(1.);
     160    if (df != 1.) {
     161      unscaled_v->Scal(df);
     162    }
    101163    return unscaled_v;
    102164  }
    103165
    104   SmartPtr<const Vector> NLPScalingObject::unapply_grad_obj_scaling(const SmartPtr<const Vector>& v)
    105   {
    106     SmartPtr<Vector> unscaled_v = unapply_grad_obj_scaling_NonConst(v);
    107     return ConstPtr(unscaled_v);
    108   }
    109 
    110   void NoNLPScalingObject::DetermineScaling(const SmartPtr<const VectorSpace> x_space,
    111       const SmartPtr<const VectorSpace> c_space,
    112       const SmartPtr<const VectorSpace> d_space,
    113       const SmartPtr<const MatrixSpace> jac_c_space,
    114       const SmartPtr<const MatrixSpace> jac_d_space,
    115       const SmartPtr<const SymMatrixSpace> h_space,
    116       SmartPtr<const MatrixSpace>& new_jac_c_space,
    117       SmartPtr<const MatrixSpace>& new_jac_d_space,
    118       SmartPtr<const SymMatrixSpace>& new_h_space)
    119   {
    120     new_jac_c_space = jac_c_space;
    121     new_jac_d_space = jac_d_space;
    122     new_h_space = h_space;
    123   }
     166  SmartPtr<const Vector> NLPScalingObject::unapply_grad_obj_scaling(
     167    const SmartPtr<const Vector>& v)
     168  {
     169    Number df = unapply_obj_scaling(1.);
     170    if (df != 1.) {
     171      SmartPtr<Vector> unscaled_v = unapply_grad_obj_scaling_NonConst(v);
     172      return ConstPtr(unscaled_v);
     173    }
     174    else {
     175      SmartPtr<const Vector> scaled_v = apply_vector_scaling_x(v);
     176      return scaled_v;
     177    }
     178  }
     179
     180  DefineIpoptType(StandardScalingBase);
     181
     182  void
     183  StandardScalingBase::RegisterOptions(SmartPtr<RegisteredOptions> roptions)
     184  {
     185    roptions->AddNumberOption(
     186      "obj_scaling_factor",
     187      "Scaling factor for the objective function.",
     188      1.,
     189      "This option allows to set a scaling factor for the objective function "
     190      "that is used to scale the problem that is seen internally by "
     191      "Ipopt. If additional scaling parameters are computed "
     192      "(e.g. user-scaling or gradient-based), this factor is multiplied "
     193      "in addition. If this value is chosen to be negative, Ipopt will "
     194      "maximize the objective function.");
     195  }
     196
     197  bool StandardScalingBase::InitializeImpl(const OptionsList& options,
     198      const std::string& prefix)
     199  {
     200    options.GetNumericValue("obj_scaling_factor", obj_scaling_factor_, prefix);
     201    return true;
     202  }
     203
     204  void StandardScalingBase::DetermineScaling(
     205    const SmartPtr<const VectorSpace> x_space,
     206    const SmartPtr<const VectorSpace> c_space,
     207    const SmartPtr<const VectorSpace> d_space,
     208    const SmartPtr<const MatrixSpace> jac_c_space,
     209    const SmartPtr<const MatrixSpace> jac_d_space,
     210    const SmartPtr<const SymMatrixSpace> h_space,
     211    SmartPtr<const MatrixSpace>& new_jac_c_space,
     212    SmartPtr<const MatrixSpace>& new_jac_d_space,
     213    SmartPtr<const SymMatrixSpace>& new_h_space)
     214  {
     215    SmartPtr<Vector> dc;
     216    SmartPtr<Vector> dd;
     217    DetermineScalingParametersImpl(x_space, c_space, d_space,
     218                                   jac_c_space, jac_d_space,
     219                                   h_space, df_, dx_, dc, dd);
     220
     221    df_ *= obj_scaling_factor_;
     222
     223    if (Jnlst().ProduceOutput(J_VECTOR, J_MAIN)) {
     224      Jnlst().Printf(J_VECTOR, J_MAIN, "objective scaling factor = %g\n", df_);
     225      if (IsValid(dx_)) {
     226        Jnlst().PrintVector(J_VECTOR, J_MAIN, "x scaling vector", *dx_);
     227      }
     228      else {
     229        Jnlst().Printf(J_VECTOR, J_MAIN, "No x scaling provided\n");
     230      }
     231      if (IsValid(dc)) {
     232        Jnlst().PrintVector(J_VECTOR, J_MAIN, "c scaling vector", *dc);
     233      }
     234      else {
     235        Jnlst().Printf(J_VECTOR, J_MAIN, "No c scaling provided\n");
     236      }
     237      if (IsValid(dd)) {
     238        Jnlst().PrintVector(J_VECTOR, J_MAIN, "d scaling vector", *dd);
     239      }
     240      else {
     241        Jnlst().Printf(J_VECTOR, J_MAIN, "No d scaling provided\n");
     242      }
     243    }
     244
     245    // create the scaling matrix spaces
     246    if (IsValid(dx_) || IsValid(dc)) {
     247      scaled_jac_c_space_ =
     248        new ScaledMatrixSpace(ConstPtr(dc), false, jac_c_space,
     249                              ConstPtr(dx_), true);
     250      new_jac_c_space = GetRawPtr(scaled_jac_c_space_);
     251    }
     252    else {
     253      scaled_jac_c_space_ = NULL;
     254      new_jac_c_space = jac_c_space;
     255    }
     256
     257    if (IsValid(dx_) || IsValid(dc)) {
     258      scaled_jac_d_space_ =
     259        new ScaledMatrixSpace(ConstPtr(dd), false, jac_d_space,
     260                              ConstPtr(dx_), true);
     261      new_jac_d_space = GetRawPtr(scaled_jac_d_space_);
     262    }
     263    else {
     264      scaled_jac_d_space_ = NULL;
     265      new_jac_d_space =jac_d_space ;
     266    }
     267
     268    if (IsValid(dx_)) {
     269      scaled_h_space_ = new SymScaledMatrixSpace(ConstPtr(dx_), true, h_space);
     270      new_h_space = GetRawPtr(scaled_h_space_);
     271    }
     272    else {
     273      scaled_h_space_ = NULL;
     274      new_h_space = h_space;
     275    }
     276  }
     277
     278  Number StandardScalingBase::apply_obj_scaling(const Number& f)
     279  {
     280    return df_*f;
     281  }
     282
     283  Number StandardScalingBase::unapply_obj_scaling(const Number& f)
     284  {
     285    return f/df_;
     286  }
     287
     288  SmartPtr<Vector> StandardScalingBase::apply_vector_scaling_x_NonConst(
     289    const SmartPtr<const Vector>& v)
     290  {
     291    DBG_START_METH("StandardScalingBase::apply_vector_scaling_x_NonConst",
     292                   dbg_verbosity);
     293    SmartPtr<Vector> scaled_x = v->MakeNewCopy();
     294    if (IsValid(dx_)) {
     295      scaled_x->ElementWiseMultiply(*dx_);
     296    }
     297    else {
     298      DBG_PRINT((1, "Creating copy in apply_vector_scaling_x_NonConst!"));
     299    }
     300    return scaled_x;
     301  };
     302
     303  SmartPtr<const Vector> StandardScalingBase::apply_vector_scaling_x(
     304    const SmartPtr<const Vector>& v)
     305  {
     306    if (IsValid(dx_)) {
     307      return ConstPtr(apply_vector_scaling_x_NonConst(v));
     308    }
     309    else {
     310      return v;
     311    }
     312  }
     313
     314  SmartPtr<Vector> StandardScalingBase::unapply_vector_scaling_x_NonConst(
     315    const SmartPtr<const Vector>& v)
     316  {
     317    DBG_START_METH("StandardScalingBase::unapply_vector_scaling_x_NonConst",
     318                   dbg_verbosity);
     319    SmartPtr<Vector> unscaled_x = v->MakeNewCopy();
     320    if (IsValid(dx_)) {
     321      unscaled_x->ElementWiseDivide(*dx_);
     322    }
     323    else {
     324      DBG_PRINT((1, "Creating copy in unapply_vector_scaling_x_NonConst!"));
     325    }
     326    return unscaled_x;
     327  }
     328
     329  SmartPtr<const Vector> StandardScalingBase::unapply_vector_scaling_x(
     330    const SmartPtr<const Vector>& v)
     331  {
     332    if (IsValid(dx_)) {
     333      return ConstPtr(unapply_vector_scaling_x_NonConst(v));
     334    }
     335    else {
     336      return v;
     337    }
     338  }
     339
     340  SmartPtr<Vector> StandardScalingBase::apply_vector_scaling_c_NonConst(
     341    const SmartPtr<const Vector>& v)
     342  {
     343    DBG_START_METH("StandardScalingBase::apply_vector_scaling_c_NonConst",
     344                   dbg_verbosity);
     345    SmartPtr<Vector> scaled_c = v->MakeNewCopy();
     346    if (IsValid(scaled_jac_c_space_) &&
     347        IsValid(scaled_jac_c_space_->RowScaling())) {
     348      scaled_c->ElementWiseMultiply(*scaled_jac_c_space_->RowScaling());
     349    }
     350    else {
     351      DBG_PRINT((1,"Creating copy in apply_vector_scaling_c_NonConst!"));
     352    }
     353    return scaled_c;
     354  }
     355
     356  SmartPtr<const Vector> StandardScalingBase::apply_vector_scaling_c(
     357    const SmartPtr<const Vector>& v)
     358  {
     359    if (IsValid(scaled_jac_c_space_) &&
     360        IsValid(scaled_jac_c_space_->RowScaling())) {
     361      return ConstPtr(apply_vector_scaling_c_NonConst(v));
     362    }
     363    else {
     364      return v;
     365    }
     366  }
     367
     368  SmartPtr<Vector> StandardScalingBase::unapply_vector_scaling_c_NonConst(
     369    const SmartPtr<const Vector>& v)
     370  {
     371    DBG_START_METH("StandardScalingBase::unapply_vector_scaling_c_NonConst",
     372                   dbg_verbosity);
     373    SmartPtr<Vector> scaled_c = v->MakeNewCopy();
     374    if (IsValid(scaled_jac_c_space_) &&
     375        IsValid(scaled_jac_c_space_->RowScaling())) {
     376      scaled_c->ElementWiseDivide(*scaled_jac_c_space_->RowScaling());
     377    }
     378    else {
     379      DBG_PRINT((1,"Creating copy in unapply_vector_scaling_c_NonConst!"));
     380    }
     381    return scaled_c;
     382  }
     383
     384  SmartPtr<const Vector> StandardScalingBase::unapply_vector_scaling_c(
     385    const SmartPtr<const Vector>& v)
     386  {
     387    if (IsValid(scaled_jac_c_space_) &&
     388        IsValid(scaled_jac_c_space_->RowScaling())) {
     389      return ConstPtr(unapply_vector_scaling_c_NonConst(v));
     390    }
     391    else {
     392      return v;
     393    }
     394  }
     395
     396  SmartPtr<Vector> StandardScalingBase::apply_vector_scaling_d_NonConst(
     397    const SmartPtr<const Vector>& v)
     398  {
     399    DBG_START_METH("StandardScalingBase::apply_vector_scaling_d_NonConst",
     400                   dbg_verbosity);
     401    SmartPtr<Vector> scaled_d = v->MakeNewCopy();
     402    if (IsValid(scaled_jac_d_space_) &&
     403        IsValid(scaled_jac_d_space_->RowScaling())) {
     404      scaled_d->ElementWiseMultiply(*scaled_jac_d_space_->RowScaling());
     405    }
     406    else {
     407      DBG_PRINT((1,"Creating copy in apply_vector_scaling_d_NonConst!"));
     408    }
     409    return scaled_d;
     410  }
     411
     412  SmartPtr<const Vector> StandardScalingBase::apply_vector_scaling_d(
     413    const SmartPtr<const Vector>& v)
     414  {
     415    if (IsValid(scaled_jac_d_space_) &&
     416        IsValid(scaled_jac_d_space_->RowScaling())) {
     417      return ConstPtr(apply_vector_scaling_d_NonConst(v));
     418    }
     419    else {
     420      return v;
     421    }
     422  }
     423
     424  SmartPtr<Vector> StandardScalingBase::unapply_vector_scaling_d_NonConst(
     425    const SmartPtr<const Vector>& v)
     426  {
     427    DBG_START_METH("StandardScalingBase::unapply_vector_scaling_d_NonConst",
     428                   dbg_verbosity);
     429    SmartPtr<Vector> scaled_d = v->MakeNewCopy();
     430    if (IsValid(scaled_jac_d_space_) &&
     431        IsValid(scaled_jac_d_space_->RowScaling())) {
     432      scaled_d->ElementWiseDivide(*scaled_jac_d_space_->RowScaling());
     433    }
     434    else {
     435      DBG_PRINT((1,"Creating copy in unapply_vector_scaling_d_NonConst!"));
     436    }
     437    return scaled_d;
     438  }
     439
     440  SmartPtr<const Vector> StandardScalingBase::unapply_vector_scaling_d(
     441    const SmartPtr<const Vector>& v)
     442  {
     443    if (IsValid(scaled_jac_d_space_) &&
     444        IsValid(scaled_jac_d_space_->RowScaling())) {
     445      return ConstPtr(unapply_vector_scaling_d_NonConst(v));
     446    }
     447    else {
     448      return v;
     449    }
     450  }
     451
     452  // ToDo: matrix not passed by reference, so setting to NULL doesn't make difference
     453  SmartPtr<const Matrix> StandardScalingBase::apply_jac_c_scaling(
     454    SmartPtr<const Matrix> matrix)
     455  {
     456    if (IsValid(scaled_jac_c_space_)) {
     457      SmartPtr<ScaledMatrix> ret = scaled_jac_c_space_->MakeNewScaledMatrix(false);
     458      ret->SetUnscaledMatrix(matrix);
     459      return GetRawPtr(ret);
     460    }
     461    else {
     462      SmartPtr<const Matrix> ret = matrix;
     463      matrix = NULL;
     464      return ret;
     465    }
     466  }
     467
     468  SmartPtr<const Matrix> StandardScalingBase::apply_jac_d_scaling(
     469    SmartPtr<const Matrix> matrix)
     470  {
     471    if (IsValid(scaled_jac_d_space_)) {
     472      SmartPtr<ScaledMatrix> ret = scaled_jac_d_space_->MakeNewScaledMatrix(false);
     473      ret->SetUnscaledMatrix(matrix);
     474      return GetRawPtr(ret);
     475    }
     476    else {
     477      SmartPtr<const Matrix> ret = matrix;
     478      matrix = NULL;
     479      return ret;
     480    }
     481  }
     482
     483  SmartPtr<const SymMatrix> StandardScalingBase::apply_hessian_scaling(
     484    SmartPtr<const SymMatrix> matrix)
     485  {
     486    if (IsValid(scaled_h_space_)) {
     487      SmartPtr<SymScaledMatrix> ret = scaled_h_space_->MakeNewSymScaledMatrix(false);
     488      ret->SetUnscaledMatrix(matrix);
     489      return GetRawPtr(ret);
     490    }
     491    else {
     492      SmartPtr<const SymMatrix> ret = matrix;
     493      matrix = NULL;
     494      return ret;
     495    }
     496  }
     497
     498  bool StandardScalingBase::have_x_scaling()
     499  {
     500    return IsValid(dx_);
     501  }
     502
     503  bool StandardScalingBase::have_c_scaling()
     504  {
     505    return (IsValid(scaled_jac_c_space_) &&
     506            IsValid(scaled_jac_c_space_->RowScaling()));
     507  }
     508
     509  bool StandardScalingBase::have_d_scaling()
     510  {
     511    return (IsValid(scaled_jac_d_space_) &&
     512            IsValid(scaled_jac_d_space_->RowScaling()));
     513  }
     514
     515  void NoNLPScalingObject::DetermineScalingParametersImpl(
     516    const SmartPtr<const VectorSpace> x_space,
     517    const SmartPtr<const VectorSpace> c_space,
     518    const SmartPtr<const VectorSpace> d_space,
     519    const SmartPtr<const MatrixSpace> jac_c_space,
     520    const SmartPtr<const MatrixSpace> jac_d_space,
     521    const SmartPtr<const SymMatrixSpace> h_space,
     522    Number& df,
     523    SmartPtr<Vector>& dx,
     524    SmartPtr<Vector>& dc,
     525    SmartPtr<Vector>& dd)
     526  {
     527    df = 1.;
     528    dx = NULL;
     529    dc = NULL;
     530    dd = NULL;
     531  }
     532
    124533} // namespace Ipopt
  • trunk/Algorithm/IpNLPScaling.hpp

    • Property svn:keywords set to Author Date Id Revision
    r346 r424  
    33// This code is published under the Common Public License.
    44//
    5 // $Id: IpNLPScaling.hpp,v 1.3 2004/11/19 00:55:18 andreasw Exp $
     5// $Id$
    66//
    77// Authors:  Carl Laird, Andreas Waechter     IBM    2004-08-13
     
    1414#include "IpMatrix.hpp"
    1515#include "IpSymMatrix.hpp"
     16#include "IpScaledMatrix.hpp"
     17#include "IpSymScaledMatrix.hpp"
     18#include "IpOptionsList.hpp"
     19#include "IpIpoptType.hpp"
    1620
    1721namespace Ipopt
     
    3539    //@}
    3640
     41    /** Method to initialize the options */
     42    bool Initialize(const Journalist& jnlst,
     43                    const OptionsList& options,
     44                    const std::string& prefix)
     45    {
     46      jnlst_ = &jnlst;
     47      return InitializeImpl(options, prefix);
     48    }
     49
    3750    /** Methods to map scaled and unscaled matrices */
    3851    //@{
     
    4255    virtual Number unapply_obj_scaling(const Number& f)=0;
    4356    /** Returns an x-scaled version of the given vector */
    44     virtual SmartPtr<Vector> apply_vector_scaling_x_NonConst(const SmartPtr<const Vector>& v)=0;
     57    virtual SmartPtr<Vector>
     58    apply_vector_scaling_x_NonConst(const SmartPtr<const Vector>& v)=0;
    4559    /** Returns an x-scaled version of the given vector */
    46     virtual SmartPtr<const Vector> apply_vector_scaling_x(const SmartPtr<const Vector>& v)=0;
     60    virtual SmartPtr<const Vector>
     61    apply_vector_scaling_x(const SmartPtr<const Vector>& v)=0;
    4762    /** Returns an x-unscaled version of the given vector */
    48     virtual SmartPtr<Vector> unapply_vector_scaling_x_NonConst(const SmartPtr<const Vector>& v)=0;
     63    virtual SmartPtr<Vector>
     64    unapply_vector_scaling_x_NonConst(const SmartPtr<const Vector>& v)=0;
    4965    /** Returns an x-unscaled version of the given vector */
    50     virtual SmartPtr<const Vector> unapply_vector_scaling_x(const SmartPtr<const Vector>& v)=0;
     66    virtual SmartPtr<const Vector>
     67    unapply_vector_scaling_x(const SmartPtr<const Vector>& v)=0;
    5168    /** Returns an c-scaled version of the given vector */
    52     virtual SmartPtr<const Vector> apply_vector_scaling_c(const SmartPtr<const Vector>& v)=0;
     69    virtual SmartPtr<const Vector>
     70    apply_vector_scaling_c(const SmartPtr<const Vector>& v)=0;
    5371    /** Returns an c-unscaled version of the given vector */
    54     virtual SmartPtr<const Vector> unapply_vector_scaling_c(const SmartPtr<const Vector>& v)=0;
     72    virtual SmartPtr<const Vector>
     73    unapply_vector_scaling_c(const SmartPtr<const Vector>& v)=0;
    5574    /** Returns an c-scaled version of the given vector */
    56     virtual SmartPtr<Vector> apply_vector_scaling_c_NonConst(const SmartPtr<const Vector>& v)=0;
     75    virtual SmartPtr<Vector>
     76    apply_vector_scaling_c_NonConst(const SmartPtr<const Vector>& v)=0;
    5777    /** Returns an c-unscaled version of the given vector */
    58     virtual SmartPtr<Vector> unapply_vector_scaling_c_NonConst(const SmartPtr<const Vector>& v)=0;
     78    virtual SmartPtr<Vector>
     79    unapply_vector_scaling_c_NonConst(const SmartPtr<const Vector>& v)=0;
    5980    /** Returns an d-scaled version of the given vector */
    60     virtual SmartPtr<const Vector> apply_vector_scaling_d(const SmartPtr<const Vector>& v)=0;
     81    virtual SmartPtr<const Vector>
     82    apply_vector_scaling_d(const SmartPtr<const Vector>& v)=0;
    6183    /** Returns an d-unscaled version of the given vector */
    62     virtual SmartPtr<const Vector> unapply_vector_scaling_d(const SmartPtr<const Vector>& v)=0;
     84    virtual SmartPtr<const Vector>
     85    unapply_vector_scaling_d(const SmartPtr<const Vector>& v)=0;
    6386    /** Returns an d-scaled version of the given vector */
    64     virtual SmartPtr<Vector> apply_vector_scaling_d_NonConst(const SmartPtr<const Vector>& v)=0;
     87    virtual SmartPtr<Vector>
     88    apply_vector_scaling_d_NonConst(const SmartPtr<const Vector>& v)=0;
    6589    /** Returns an d-unscaled version of the given vector */
    66     virtual SmartPtr<Vector> unapply_vector_scaling_d_NonConst(const SmartPtr<const Vector>& v)=0;
    67     /** Returns a scaled version of the jacobian for c.
    68      *  If the overloaded method does not make a new matrix, make sure to set the matrix
    69      *  ptr passed in to NULL.
    70      */
    71     virtual SmartPtr<const Matrix> apply_jac_c_scaling(SmartPtr<const Matrix> matrix)=0;
    72     /** Returns a scaled version of the jacobian for d
    73      *  If the overloaded method does not create a new matrix, make sure to set the matrix
    74      *  ptr passed in to NULL.
    75      */
    76     virtual SmartPtr<const Matrix> apply_jac_d_scaling(SmartPtr<const Matrix> matrix)=0;
    77     /** Returns a scaled version of the hessian of the lagrangian
    78      *  If the overloaded method does not create a new matrix, make sure to set the matrix
    79      *  ptr passed in to NULL.
    80      */
    81     virtual SmartPtr<const SymMatrix> apply_hessian_scaling(SmartPtr<const SymMatrix> matrix)=0;
     90    virtual SmartPtr<Vector>
     91    unapply_vector_scaling_d_NonConst(const SmartPtr<const Vector>& v)=0;
     92    /** Returns a scaled version of the jacobian for c.  If the
     93     *  overloaded method does not make a new matrix, make sure to set
     94     *  the matrix ptr passed in to NULL.
     95     */
     96    virtual SmartPtr<const Matrix>
     97    apply_jac_c_scaling(SmartPtr<const Matrix> matrix)=0;
     98    /** Returns a scaled version of the jacobian for d If the
     99     *  overloaded method does not create a new matrix, make sure to
     100     *  set the matrix ptr passed in to NULL.
     101     */
     102    virtual SmartPtr<const Matrix>
     103    apply_jac_d_scaling(SmartPtr<const Matrix> matrix)=0;
     104    /** Returns a scaled version of the hessian of the lagrangian If
     105     *  the overloaded method does not create a new matrix, make sure
     106     *  to set the matrix ptr passed in to NULL.
     107     */
     108    virtual SmartPtr<const SymMatrix>
     109    apply_hessian_scaling(SmartPtr<const SymMatrix> matrix)=0;
    82110    //@}
    83111
    84112    /** Methods for scaling bounds - these wrap those above */
    85113    //@{
    86     /** Returns an x-scaled vector in the x_L space */
    87     SmartPtr<Vector> apply_vector_scaling_x_L_NonConst(SmartPtr<Matrix> Px_L, const SmartPtr<const Vector>& l, const SmartPtr<const VectorSpace> x_space);
    88     /** Returns an x-scaled vector in the x_U space */
    89     SmartPtr<Vector> apply_vector_scaling_x_U_NonConst(SmartPtr<Matrix> Px_U, const SmartPtr<const Vector>& u, const SmartPtr<const VectorSpace> x_space);
    90     /** Returns an d-scaled vector in the d_L space */
    91     SmartPtr<Vector> apply_vector_scaling_d_L_NonConst(SmartPtr<Matrix> Pd_L, const SmartPtr<const Vector>& l, const SmartPtr<const VectorSpace> d_space);
    92     /** Returns an d-scaled vector in the d_U space */
    93     SmartPtr<Vector> apply_vector_scaling_d_U_NonConst(SmartPtr<Matrix> Pd_U, const SmartPtr<const Vector>& u, const SmartPtr<const VectorSpace> d_space);
     114    /** Returns an x-scaled vector in the x_L or x_U space */
     115    SmartPtr<Vector> apply_vector_scaling_x_LU_NonConst(
     116      const Matrix& Px_LU,
     117      const SmartPtr<const Vector>& lu,
     118      const VectorSpace& x_space);
     119    /** Returns an x-scaled vector in the x_L or x_U space */
     120    SmartPtr<const Vector> apply_vector_scaling_x_LU(
     121      const Matrix& Px_LU,
     122      const SmartPtr<const Vector>& lu,
     123      const VectorSpace& x_space);
     124    /** Returns an d-scaled vector in the d_L or d_U space */
     125    SmartPtr<Vector> apply_vector_scaling_d_LU_NonConst(
     126      const Matrix& Pd_LU,
     127      const SmartPtr<const Vector>& lu,
     128      const VectorSpace& d_space);
     129    /** Returns an d-scaled vector in the d_L or d_U space */
     130    SmartPtr<const Vector> apply_vector_scaling_d_LU(
     131      const Matrix& Pd_LU,
     132      const SmartPtr<const Vector>& lu,
     133      const VectorSpace& d_space);
     134    /** Returns an d-unscaled vector in the d_L or d_U space */
     135    SmartPtr<Vector> unapply_vector_scaling_d_LU_NonConst(
     136      const Matrix& Pd_LU,
     137      const SmartPtr<const Vector>& lu,
     138      const VectorSpace& d_space);
     139    /** Returns an d-unscaled vector in the d_L or d_U space */
     140    SmartPtr<const Vector> unapply_vector_scaling_d_LU(
     141      const Matrix& Pd_LU,
     142      const SmartPtr<const Vector>& lu,
     143      const VectorSpace& d_space);
    94144    //@}
    95145
     
    99149    //@{
    100150    /** Returns a grad_f scaled version (d_f * D_x^{-1}) of the given vector */
    101     virtual SmartPtr<Vector> apply_grad_obj_scaling_NonConst(const SmartPtr<const Vector>& v);
     151    virtual SmartPtr<Vector>
     152    apply_grad_obj_scaling_NonConst(const SmartPtr<const Vector>& v);
    102153    /** Returns a grad_f scaled version (d_f * D_x^{-1}) of the given vector */
    103     virtual SmartPtr<const Vector> apply_grad_obj_scaling(const SmartPtr<const Vector>& v);
     154    virtual SmartPtr<const Vector>
     155    apply_grad_obj_scaling(const SmartPtr<const Vector>& v);
    104156    /** Returns a grad_f unscaled version (d_f * D_x^{-1}) of the
    105157     *  given vector */
    106     virtual SmartPtr<Vector> unapply_grad_obj_scaling_NonConst(const SmartPtr<const Vector>& v);
    107     /** Returns a grad_f unscaled version (d_f * D_x^{-1}) of the given vector */
    108     virtual SmartPtr<const Vector> unapply_grad_obj_scaling(const SmartPtr<const Vector>& v);
     158    virtual SmartPtr<Vector>
     159    unapply_grad_obj_scaling_NonConst(const SmartPtr<const Vector>& v);
     160    /** Returns a grad_f unscaled version (d_f * D_x^{-1}) of the
     161     *  given vector */
     162    virtual SmartPtr<const Vector>
     163    unapply_grad_obj_scaling(const SmartPtr<const Vector>& v);
     164    //@}
     165
     166    /** @name Methods for determining whether scaling for entities is
     167     *  done */
     168    //@{
     169    /** Returns true if the primal x variables are scaled. */
     170    virtual bool have_x_scaling()=0;
     171    /** Returns true if the equality constraints are scaled. */
     172    virtual bool have_c_scaling()=0;
     173    /** Returns true if the inequality constraints are scaled. */
     174    virtual bool have_d_scaling()=0;
    109175    //@}
    110176
     
    121187                                  SmartPtr<const MatrixSpace>& new_jac_d_space,
    122188                                  SmartPtr<const SymMatrixSpace>& new_h_space)=0;
     189  protected:
     190    /** Implementation of the initialization method that has to be
     191     *  overloaded by for each derived class. */
     192    virtual bool InitializeImpl(const OptionsList& options,
     193                                const std::string& prefix)=0;
     194
     195    /** Accessor method for the journalist */
     196    const Journalist& Jnlst() const
     197    {
     198      return *jnlst_;
     199    }
    123200  private:
    124201
     
    138215    void operator=(const NLPScalingObject&);
    139216    //@}
     217
     218    SmartPtr<const Journalist> jnlst_;
    140219  };
    141220
    142 
    143   class NoNLPScalingObject : public NLPScalingObject
     221  DeclareIpoptType(StandardScalingBase);
     222
     223  /** This is a base class for many standard scaling
     224   *  techniques. The overloaded classes only need to
     225   *  provide the scaling parameters
     226   */
     227  class StandardScalingBase : public NLPScalingObject
    144228  {
    145229  public:
    146230    /**@name Constructors/Destructors */
    147231    //@{
    148     NoNLPScalingObject()
     232    StandardScalingBase()
    149233    {}
    150234
    151235    /** Default destructor */
    152     virtual ~NoNLPScalingObject()
    153     {}
    154     //@}
    155 
    156     /** Methods to map scaled and unscaled matrices - overloaded from NLPScalingObject */
     236    virtual ~StandardScalingBase()
     237    {}
     238    //@}
     239
     240    /** Methods to map scaled and unscaled matrices */
    157241    //@{
    158242    /** Returns an obj-scaled version of the given scalar */
    159     virtual Number apply_obj_scaling(const Number& f)
    160     {
    161       return f;
    162     }
     243    virtual Number apply_obj_scaling(const Number& f);
    163244    /** Returns an obj-unscaled version of the given scalar */
    164     virtual Number unapply_obj_scaling(const Number& f)
    165     {
    166       return f;
    167     }
     245    virtual Number unapply_obj_scaling(const Number& f);
    168246    /** Returns an x-scaled version of the given vector */
    169     virtual SmartPtr<Vector> apply_vector_scaling_x_NonConst(const SmartPtr<const Vector>& v)
    170     {
    171       DBG_ASSERT(true);
    172       return v->MakeNewCopy();
    173     }
     247    virtual SmartPtr<Vector>
     248    apply_vector_scaling_x_NonConst(const SmartPtr<const Vector>& v);
    174249    /** Returns an x-scaled version of the given vector */
    175     virtual SmartPtr<const Vector> apply_vector_scaling_x(const SmartPtr<const Vector>& v)
    176     {
    177       return v;
    178     }
     250    virtual SmartPtr<const Vector>
     251    apply_vector_scaling_x(const SmartPtr<const Vector>& v);
    179252    /** Returns an x-unscaled version of the given vector */
    180     virtual SmartPtr<Vector> unapply_vector_scaling_x_NonConst(const SmartPtr<const Vector>& v)
    181     {
    182       DBG_ASSERT(true);
    183       return v->MakeNewCopy();
    184     }
     253    virtual SmartPtr<Vector>
     254    unapply_vector_scaling_x_NonConst(const SmartPtr<const Vector>& v);
    185255    /** Returns an x-unscaled version of the given vector */
    186     virtual SmartPtr<const Vector> unapply_vector_scaling_x(const SmartPtr<const Vector>& v)
    187     {
    188       return v;
    189     }
     256    virtual SmartPtr<const Vector>
     257    unapply_vector_scaling_x(const SmartPtr<const Vector>& v);
    190258    /** Returns an c-scaled version of the given vector */
    191     virtual SmartPtr<const Vector> apply_vector_scaling_c(const SmartPtr<const Vector>& v)
    192     {
    193       return v;
    194     }
     259    virtual SmartPtr<const Vector>
     260    apply_vector_scaling_c(const SmartPtr<const Vector>& v);
     261    /** Returns an c-unscaled version of the given vector */
     262    virtual SmartPtr<const Vector>
     263    unapply_vector_scaling_c(const SmartPtr<const Vector>& v);
    195264    /** Returns an c-scaled version of the given vector */
    196     virtual SmartPtr<Vector> apply_vector_scaling_c_NonConst(const SmartPtr<const Vector>& v)
    197     {
    198       DBG_ASSERT(true);
    199       return v->MakeNewCopy();
    200     }
     265    virtual SmartPtr<Vector>
     266    apply_vector_scaling_c_NonConst(const SmartPtr<const Vector>& v);
    201267    /** Returns an c-unscaled version of the given vector */
    202     virtual SmartPtr<const Vector> unapply_vector_scaling_c(const SmartPtr<const Vector>& v)
    203     {
    204       return v;
    205     }
    206     /** Returns an c-unscaled version of the given vector */
    207     virtual SmartPtr<Vector> unapply_vector_scaling_c_NonConst(const SmartPtr<const Vector>& v)
    208     {
    209       DBG_ASSERT(true);
    210       return v->MakeNewCopy();
    211     }
     268    virtual SmartPtr<Vector>
     269    unapply_vector_scaling_c_NonConst(const SmartPtr<const Vector>& v);
    212270    /** Returns an d-scaled version of the given vector */
    213     virtual SmartPtr<const Vector> apply_vector_scaling_d(const SmartPtr<const Vector>& v)
    214     {
    215       return v;
    216     }
     271    virtual SmartPtr<const Vector>
     272    apply_vector_scaling_d(const SmartPtr<const Vector>& v);
     273    /** Returns an d-unscaled version of the given vector */
     274    virtual SmartPtr<const Vector>
     275    unapply_vector_scaling_d(const SmartPtr<const Vector>& v);
    217276    /** Returns an d-scaled version of the given vector */
    218     virtual SmartPtr<Vector> apply_vector_scaling_d_NonConst(const SmartPtr<const Vector>& v)
    219     {
    220       DBG_ASSERT(true);
    221       return v->MakeNewCopy();
    222     }
     277    virtual SmartPtr<Vector>
     278    apply_vector_scaling_d_NonConst(const SmartPtr<const Vector>& v);
    223279    /** Returns an d-unscaled version of the given vector */
    224     virtual SmartPtr<const Vector> unapply_vector_scaling_d(const SmartPtr<const Vector>& v)
    225     {
    226       return v;
    227     }
    228     /** Returns an d-unscaled version of the given vector */
    229     virtual SmartPtr<Vector> unapply_vector_scaling_d_NonConst(const SmartPtr<const Vector>& v)
    230     {
    231       DBG_ASSERT(true);
    232       return v->MakeNewCopy();
    233     }
    234     /** Returns a scaled version of the jacobian for c.
    235      *  If the overloaded method does not create a new matrix, make sure to set the matrix
    236      *  ptr passed in to NULL.
    237      */
    238     virtual SmartPtr<const Matrix> apply_jac_c_scaling(SmartPtr<const Matrix> matrix)
    239     {
    240       SmartPtr<const Matrix> ret = matrix;
    241       matrix = NULL;
    242       return ret;
    243     }
    244     /** Returns a scaled version of the jacobian for d
    245      *  If the overloaded method does not create a new matrix, make sure to set the matrix
    246      *  ptr passed in to NULL.
    247      */
    248     virtual SmartPtr<const Matrix> apply_jac_d_scaling(SmartPtr<const Matrix> matrix)
    249     {
    250       SmartPtr<const Matrix> ret = matrix;
    251       matrix = NULL;
    252       return ret;
    253     }
    254     /** Returns a scaled version of the hessian of the lagrangian
    255      *  If the overloaded method does not create a new matrix, make sure to set the matrix
    256      *  ptr passed in to NULL.
    257      */
    258     virtual SmartPtr<const SymMatrix> apply_hessian_scaling(SmartPtr<const SymMatrix> matrix)
    259     {
    260       SmartPtr<const SymMatrix> ret = matrix;
    261       matrix = NULL;
    262       return ret;
    263     }
    264     //@}
    265 
    266     /** Method to determine the scaling - overloaded from NLPScaling.
    267      *  Since no scaling is done, we do not need to determine anything
     280    virtual SmartPtr<Vector>
     281    unapply_vector_scaling_d_NonConst(const SmartPtr<const Vector>& v);
     282    /** Returns a scaled version of the jacobian for c.  If the
     283     *  overloaded method does not make a new matrix, make sure to set
     284     *  the matrix ptr passed in to NULL.
     285     */
     286    virtual SmartPtr<const Matrix>
     287    apply_jac_c_scaling(SmartPtr<const Matrix> matrix);
     288    /** Returns a scaled version of the jacobian for d If the
     289     *  overloaded method does not create a new matrix, make sure to
     290     *  set the matrix ptr passed in to NULL.
     291     */
     292    virtual SmartPtr<const Matrix>
     293    apply_jac_d_scaling(SmartPtr<const Matrix> matrix);
     294    /** Returns a scaled version of the hessian of the lagrangian If
     295     *  the overloaded method does not create a new matrix, make sure
     296     *  to set the matrix ptr passed in to NULL.
     297     */
     298    virtual SmartPtr<const SymMatrix>
     299    apply_hessian_scaling(SmartPtr<const SymMatrix> matrix);
     300    //@}
     301
     302    /** @name Methods for determining whether scaling for entities is
     303     *  done */
     304    //@{
     305    virtual bool have_x_scaling();
     306    virtual bool have_c_scaling();
     307    virtual bool have_d_scaling();
     308    //@}
     309
     310    /** This method is called by the IpoptNLP's at a convenient time to
     311     *  compute and/or read scaling factors
    268312     */
    269313    virtual void DetermineScaling(const SmartPtr<const VectorSpace> x_space,
     
    277321                                  SmartPtr<const SymMatrixSpace>& new_h_space);
    278322
    279     /** Methods for scaling the gradient of the objective - made efficient for
    280      *  the NoScaling implementation.
    281      */
    282     //@{
    283     /** Returns a grad_f scaled version (d_f * D_x^{-1}) of the given vector */
    284     virtual SmartPtr<Vector> apply_grad_obj_scaling_NonConst(const SmartPtr<const Vector>& v)
    285     {
    286       DBG_ASSERT(true);
    287       return v->MakeNewCopy();
    288     }
    289     /** Returns a grad_f scaled version (d_f * D_x^{-1}) of the given vector */
    290     virtual SmartPtr<const Vector> apply_grad_obj_scaling(const SmartPtr<const Vector>& v)
    291     {
    292       return v;
    293     }
    294     /** Returns a grad_f unscaled version (d_f * D_x^{-1}) of the given vector */
    295     virtual SmartPtr<Vector> unapply_grad_obj_scaling_NonConst(const SmartPtr<const Vector>& v)
    296     {
    297       DBG_ASSERT(true);
    298       return v->MakeNewCopy();
    299     }
    300     /** Returns a grad_f unscaled version (d_f * D_x^{-1}) of the given vector */
    301     virtual SmartPtr<const Vector> unapply_grad_obj_scaling(const SmartPtr<const Vector>& v)
    302     {
    303       return v;
    304     }
    305     //@}
     323    /** Methods for IpoptType */
     324    //@{
     325    static void RegisterOptions(SmartPtr<RegisteredOptions> roptions);
     326    //@}
     327
     328  protected:
     329    /** Overloaded initialization method */
     330    virtual bool InitializeImpl(const OptionsList& options,
     331                                const std::string& prefix);
     332
     333    /** This is the method that has to be overloaded by a particular
     334     *  scaling method that somehow computes the scaling vectors dx,
     335     *  dc, and dd.  The pointers to those vectors can be NULL, in
     336     *  which case no scaling for that item will be done later. */
     337    virtual void DetermineScalingParametersImpl(
     338      const SmartPtr<const VectorSpace> x_space,
     339      const SmartPtr<const VectorSpace> c_space,
     340      const SmartPtr<const VectorSpace> d_space,
     341      const SmartPtr<const MatrixSpace> jac_c_space,
     342      const SmartPtr<const MatrixSpace> jac_d_space,
     343      const SmartPtr<const SymMatrixSpace> h_space,
     344      Number& df,
     345      SmartPtr<Vector>& dx,
     346      SmartPtr<Vector>& dc,
     347      SmartPtr<Vector>& dd)=0;
    306348
    307349  private:
     
    317359
    318360    /** Copy Constructor */
     361    StandardScalingBase(const StandardScalingBase&);
     362
     363    /** Overloaded Equals Operator */
     364    void operator=(const StandardScalingBase&);
     365    //@}
     366
     367    /** Scaling parameters - we only need to keep copies of
     368     *  the objective scaling and the x scaling - the others we can
     369     *  get from the scaled matrix spaces.
     370     */
     371    //@{
     372    /** objective scaling parameter */
     373    Number df_;
     374    /** x scaling */
     375    SmartPtr<Vector> dx_;
     376    //@}
     377
     378    /** Scaled Matrix Spaces */
     379    //@{
     380    /** Scaled jacobian of c space */
     381    SmartPtr<ScaledMatrixSpace> scaled_jac_c_space_;
     382    /** Scaled jacobian of d space */
     383    SmartPtr<ScaledMatrixSpace> scaled_jac_d_space_;
     384    /** Scaled hessian of lagrangian spacea */
     385    SmartPtr<SymScaledMatrixSpace> scaled_h_space_;
     386    //@}
     387
     388    /** @name Algorithmic parameters */
     389    //@{
     390    /** Additional scaling value for the objective function */
     391    Number obj_scaling_factor_;
     392    //@}
     393  };
     394
     395  /** Class implementing the scaling object that doesn't to any scaling */
     396  class NoNLPScalingObject : public StandardScalingBase
     397  {
     398  public:
     399    /**@name Constructors/Destructors */
     400    //@{
     401    NoNLPScalingObject()
     402    {}
     403
     404    /** Default destructor */
     405    virtual ~NoNLPScalingObject()
     406    {}
     407    //@}
     408
     409
     410  protected:
     411    /** Overloaded from StandardScalingBase */
     412    virtual void DetermineScalingParametersImpl(
     413      const SmartPtr<const VectorSpace> x_space,
     414      const SmartPtr<const VectorSpace> c_space,
     415      const SmartPtr<const VectorSpace> d_space,
     416      const SmartPtr<const MatrixSpace> jac_c_space,
     417      const SmartPtr<const MatrixSpace> jac_d_space,
     418      const SmartPtr<const SymMatrixSpace> h_space,
     419      Number& df,
     420      SmartPtr<Vector>& dx,
     421      SmartPtr<Vector>& dc,
     422      SmartPtr<Vector>& dd);
     423
     424  private:
     425
     426    /**@name Default Compiler Generated Methods
     427     * (Hidden to avoid implicit creation/calling).
     428     * These methods are not implemented and
     429     * we do not want the compiler to implement
     430     * them for us, so we declare them private
     431     * and do not define them. This ensures that
     432     * they will not be implicitly created/called. */
     433    //@{
     434
     435    /** Copy Constructor */
    319436    NoNLPScalingObject(const NoNLPScalingObject&);
    320437
  • trunk/Algorithm/IpOptErrorConvCheck.cpp

    r321 r424  
    2323  void OptimalityErrorConvergenceCheck::RegisterOptions(SmartPtr<RegisteredOptions> roptions)
    2424  {
    25     roptions->AddLowerBoundedIntegerOption("maxiter", "Maximum number of iterations to allow",
    26                                            0, 1000);
     25    roptions->AddLowerBoundedIntegerOption(
     26      "max_iter",
     27      "Maximum number of iterations.",
     28      0, 3000,
     29      "The algorithm terminates with an error message if the number of "
     30      "iterations exceeded this number. [Also used in RestoFilterConvCheck]");
     31    roptions->AddLowerBoundedNumberOption(
     32      "dual_inf_tol",
     33      "Acceptance threshold for the dual infeasibility.",
     34      0.0, true, 1e-4,
     35      "Absolute tolerance on the dual infesaibility. Successful termination "
     36      "requires that the (unscaled) dual infeasibility is less than this "
     37      "threshold.");
     38    roptions->AddLowerBoundedNumberOption(
     39      "constr_viol_tol",
     40      "Acceptance threshold for the constraint violation.",
     41      0.0, true, 1e-4,
     42      "Absolute tolerance on the constraint violation. Successful termination "
     43      "requires that the (unscaled) constraint violation is less than this "
     44      "threshold.");
     45    roptions->AddLowerBoundedNumberOption(
     46      "compl_inf_tol",
     47      "Acceptance threshold for the complementarity conditions.",
     48      0.0, true, 1e-4,
     49      "Absolute tolerance on the complementarity. Successful termination "
     50      "requires that the (unscaled) complementarity is less than this "
     51      "threshold.");
     52    roptions->AddLowerBoundedIntegerOption(
     53      "acceptable_iter",
     54      "Number of acceptable iterates to trigger termination.",
     55      0, 15,
     56      "If the algorithm encounters so many successive acceptable iterates "
     57      "(see \"acceptable_tol\"), it terminates, assuming that the problem "
     58      "has been solved to best possible accuracy given round-off.  If it is "
     59      "set to zero, this heuristic is disabled.");
     60    roptions->AddLowerBoundedNumberOption(
     61      "acceptable_tol",
     62      "Acceptable convergence tolerance (relative).",
     63      0.0, true,  1e-6,
     64      "Determines which (scaled) overall optimality error is considered to be "
     65      "acceptable. If the algorithm encounters \"acceptable_iter\" iterations "
     66      "in a row that are considered \"acceptable\", it will terminate before "
     67      "the desired convergence tolerance (\"tol\", \"dual_inf_tol\", etc) is "
     68      "met.");
     69    roptions->AddLowerBoundedNumberOption(
     70      "acceptable_dual_inf_tol",
     71      "Acceptance threshold for the dual infeasibility.",
     72      0.0, true, 1e-2,
     73      "Absolute tolerance on the dual infesaibility. Acceptable termination "
     74      "requires that the (unscaled) dual infeasibility is less than this "
     75      "threshold.");
     76    roptions->AddLowerBoundedNumberOption(
     77      "acceptable_constr_viol_tol",
     78      "Acceptance threshold for the constraint violation.",
     79      0.0, true, 1e-2,
     80      "Absolute tolerance on the constraint violation. Acceptable termination "
     81      "requires that the (unscaled) constraint violation is less than this "
     82      "threshold.");
     83    roptions->AddLowerBoundedNumberOption(
     84      "acceptable_compl_inf_tol",
     85      "Acceptance threshold for the complementarity conditions.",
     86      0.0, true, 1e-2,
     87      "Absolute tolerance on the complementarity. Acceptable termination "
     88      "requires that the (unscaled) complementarity is less than this "
     89      "threshold.");
    2790  }
    2891
     
    3194      const std::string& prefix)
    3295  {
    33     options.GetIntegerValue("maxiter", max_iterations_, prefix);
     96    options.GetIntegerValue("max_iter", max_iterations_, prefix);
     97    options.GetNumericValue("dual_inf_tol", dual_inf_tol_, prefix);
     98    options.GetNumericValue("constr_viol_tol", constr_viol_tol_, prefix);
     99    options.GetNumericValue("compl_inf_tol", compl_inf_tol_, prefix);
     100    options.GetIntegerValue("acceptable_iter", acceptable_iter_, prefix);
     101    options.GetNumericValue("acceptable_tol", acceptable_tol_, prefix);
     102    options.GetNumericValue("acceptable_dual_inf_tol", acceptable_dual_inf_tol_, prefix);
     103    options.GetNumericValue("acceptable_constr_viol_tol", acceptable_constr_viol_tol_, prefix);
     104    options.GetNumericValue("acceptable_compl_inf_tol", acceptable_compl_inf_tol_, prefix);
     105    acceptable_counter_ = 0;
    34106
    35107    return true;
    36108  }
    37109
    38   ConvergenceCheck::ConvergenceStatus OptimalityErrorConvergenceCheck::CheckConvergence()
     110  ConvergenceCheck::ConvergenceStatus
     111  OptimalityErrorConvergenceCheck::CheckConvergence()
    39112  {
    40113    DBG_START_METH("OptimalityErrorConvergenceCheck::CheckConvergence", dbg_verbosity);
    41     // maybe we should throw exceptions here instead?
    42114
    43115    if (IpData().iter_count() >= max_iterations_) {
     
    46118
    47119    Number overall_error = IpCq().curr_nlp_error();
    48     Number dual_inf = IpCq().curr_dual_infeasibility(NORM_MAX);
    49     Number primal_inf = IpCq().curr_primal_infeasibility(NORM_MAX);
    50     Number compl_inf = IpCq().curr_complementarity(0., NORM_MAX);
    51     DBG_PRINT((1,"overall_error = %8.2e\n",overall_error));
     120    Number dual_inf = IpCq().unscaled_curr_dual_infeasibility(NORM_MAX);
     121    Number constr_viol = IpCq().unscaled_curr_nlp_constraint_violation(NORM_MAX);
     122    Number compl_inf = IpCq().unscaled_curr_complementarity(0., NORM_MAX);
     123
    52124    if (overall_error <= IpData().tol() &&
    53         dual_inf <= IpData().dual_inf_tol() &&
    54         primal_inf <= IpData().primal_inf_tol() &&
    55         compl_inf <= IpData().compl_inf_tol()) {
     125        dual_inf <= dual_inf_tol_ &&
     126        constr_viol <= constr_viol_tol_ &&
     127        compl_inf <= compl_inf_tol_) {
    56128      return ConvergenceCheck::CONVERGED;
     129    }
     130
     131    if (acceptable_iter_>0 && CurrentIsAcceptable()) {
     132      IpData().Append_info_string("A");
     133      acceptable_counter_++;
     134      if (acceptable_counter_ >= acceptable_iter_) {
     135        return ConvergenceCheck::CONVERGED_TO_ACCEPTABLE_POINT;
     136      }
     137    }
     138    else {
     139      acceptable_counter_ = 0;
    57140    }
    58141
    59142    return ConvergenceCheck::CONTINUE;
    60143  }
     144
     145  bool OptimalityErrorConvergenceCheck::CurrentIsAcceptable()
     146  {
     147    DBG_START_METH("OptimalityErrorConvergenceCheck::CurrentIsAcceptable",
     148                   dbg_verbosity);
     149
     150    Number overall_error = IpCq().curr_nlp_error();
     151    Number dual_inf = IpCq().unscaled_curr_dual_infeasibility(NORM_MAX);
     152    Number constr_viol = IpCq().unscaled_curr_nlp_constraint_violation(NORM_MAX);
     153    Number compl_inf = IpCq().unscaled_curr_complementarity(0., NORM_MAX);
     154
     155    DBG_PRINT((1, "overall_error = %e\n", overall_error));
     156    DBG_PRINT((1, "dual_inf = %e\n", dual_inf));
     157    DBG_PRINT((1, "constr_viol = %e\n", constr_viol));
     158    DBG_PRINT((1, "compl_inf = %e\n", compl_inf));
     159
     160    DBG_PRINT((1, "acceptable_tol_ = %e\n", acceptable_tol_));
     161    DBG_PRINT((1, "acceptable_dual_inf_tol_ = %e\n", acceptable_dual_inf_tol_));
     162    DBG_PRINT((1, "acceptable_constr_viol_tol_ = %e\n", acceptable_constr_viol_tol_));
     163    DBG_PRINT((1, "acceptable_compl_inf_tol_ = %e\n", acceptable_compl_inf_tol_));
     164
     165    return (overall_error <= acceptable_tol_ &&
     166            dual_inf <= acceptable_dual_inf_tol_ &&
     167            constr_viol <= acceptable_constr_viol_tol_ &&
     168            compl_inf <= acceptable_compl_inf_tol_);
     169  }
     170
     171
    61172} // namespace Ipopt
  • trunk/Algorithm/IpOptErrorConvCheck.hpp

    r310 r424  
    4040    virtual ConvergenceStatus CheckConvergence();
    4141
     42    /** Auxilliary function for testing whether current iterate
     43     *  satisfies the acceptable level of optimality */
     44    virtual bool CurrentIsAcceptable();
     45
    4246    /** Methods for IpoptType */
    4347    //@{
    4448    static void RegisterOptions(SmartPtr<RegisteredOptions> roptions);
     49    //@}
     50
     51  protected:
     52    /** @name Algorithmic parameters */
     53    //@{
     54    /** Maximal number of iterations */
     55    Index max_iterations_;
     56    /** Tolerance on unscaled dual infeasibility */
     57    Number dual_inf_tol_;
     58    /** Tolerance on unscaled constraint violation */
     59    Number constr_viol_tol_;
     60    /** Tolerance on unscaled complementarity */
     61    Number compl_inf_tol_;
     62    /** Number of iterations with acceptable level of accuracy, after
     63     *  which the algorithm terminates.  If 0, this heuristic is
     64     *  disabled. */
     65    Index acceptable_iter_;
     66    /** Acceptable tolerance for the problem to terminate earlier if
     67     *  algorithm seems stuck or cycling */
     68    Number acceptable_tol_;
     69    /** Acceptable tolerance on unscaled dual infeasibility */
     70    Number acceptable_dual_inf_tol_;
     71    /** Acceptable tolerance on unscaled constraint violation */
     72    Number acceptable_constr_viol_tol_;
     73    /** Acceptable tolerance on unscaled complementarity */
     74    Number acceptable_compl_inf_tol_;
    4575    //@}
    4676
     
    5989    //@}
    6090
    61     /** @name Algorithmic parameters */
    62     //@{
    63     /** Maximal number of iterations */
    64     Index max_iterations_;
    65     //@}
     91    /** Counter for successive iterations in which acceptability
     92     *  criteria are met. */
     93    Index acceptable_counter_;
    6694  };
    6795
  • trunk/Algorithm/IpOrigIpoptNLP.cpp

    r340 r424  
    4949  void OrigIpoptNLP::RegisterOptions(SmartPtr<RegisteredOptions> roptions)
    5050  {
    51     roptions->AddLowerBoundedNumberOption("bound_relax_factor","factor for initial relaxation of the bounds",
    52                                           0, false, 1e-8);
     51    roptions->AddLowerBoundedNumberOption(
     52      "bound_relax_factor",
     53      "Factor for initial relaxation of the bounds.",
     54      0, false,
     55      1e-8,
     56      "Before start of the optimization, the bounds given by the user are "
     57      "relaxed.  This option determines the factor by how much.  If it "
     58      "is set to zero, then this option is disabled.  (See Eqn.(35) in "
     59      "implmentation paper.)");
    5360  }
    5461
     
    6471
    6572    initialized_ = true;
    66     return true;
     73    return IpoptNLP::Initialize(jnlst, options, prefix);
    6774  }
    6875
     
    7885                                          bool init_z_U,
    7986                                          SmartPtr<Vector>& v_L,
    80                                           bool init_v_L,
    81                                           SmartPtr<Vector>& v_U,
    82                                           bool init_v_U
     87                                          SmartPtr<Vector>& v_U
    8388                                         )
    8489  {
     
    123128
    124129    // Create the bounds structures
    125     x_L_ = x_l_space_->MakeNew();
    126     Px_L_ = px_l_space_->MakeNew();
    127     x_U_ = x_u_space_->MakeNew();
    128     Px_U_ = px_u_space_->MakeNew();
    129     d_L_ = d_l_space_->MakeNew();
    130     Pd_L_ = pd_l_space_->MakeNew();
    131     d_U_ = d_u_space_->MakeNew();
    132     Pd_U_ = pd_u_space_->MakeNew();
    133 
    134     retValue = nlp_->GetBoundsInformation(*Px_L_, *x_L_, *Px_U_, *x_U_,
    135                                           *Pd_L_, *d_L_, *Pd_U_, *d_U_);
     130    SmartPtr<Vector> x_L = x_l_space_->MakeNew();
     131    SmartPtr<Matrix> Px_L = px_l_space_->MakeNew();
     132    SmartPtr<Vector> x_U = x_u_space_->MakeNew();
     133    SmartPtr<Matrix> Px_U = px_u_space_->MakeNew();
     134    SmartPtr<Vector> d_L = d_l_space_->MakeNew();
     135    SmartPtr<Matrix> Pd_L = pd_l_space_->MakeNew();
     136    SmartPtr<Vector> d_U = d_u_space_->MakeNew();
     137    SmartPtr<Matrix> Pd_U = pd_u_space_->MakeNew();
     138
     139    retValue = nlp_->GetBoundsInformation(*Px_L, *x_L, *Px_U, *x_U,
     140                                          *Pd_L, *d_L, *Pd_U, *d_U);
    136141
    137142    if (!retValue) {
     
    139144    }
    140145
    141     relax_bounds(-bound_relax_factor_, *x_L_);
    142     relax_bounds( bound_relax_factor_, *x_U_);
    143     relax_bounds(-bound_relax_factor_, *d_L_);
    144     relax_bounds( bound_relax_factor_, *d_U_);
     146    relax_bounds(-bound_relax_factor_, *x_L);
     147    relax_bounds( bound_relax_factor_, *x_U);
     148    relax_bounds(-bound_relax_factor_, *d_L);
     149    relax_bounds( bound_relax_factor_, *d_U);
     150
     151    x_L_ = ConstPtr(x_L);
     152    Px_L_ = ConstPtr(Px_L);
     153    x_U_ = ConstPtr(x_U);
     154    Px_U_ = ConstPtr(Px_U);
     155    d_L_ = ConstPtr(d_L);
     156    Pd_L_ = ConstPtr(Pd_L);
     157    d_U_ = ConstPtr(d_U);
     158    Pd_U_ = ConstPtr(Pd_U);
    145159
    146160    // now create and store the scaled bounds
    147     x_L_ = NLP_scaling()->apply_vector_scaling_x_L_NonConst(Px_L_, ConstPtr(x_L_), x_space_);
    148     x_U_ = NLP_scaling()->apply_vector_scaling_x_U_NonConst(Px_U_, ConstPtr(x_U_), x_space_);
    149     d_L_ = NLP_scaling()->apply_vector_scaling_d_L_NonConst(Pd_L_, ConstPtr(d_L_), d_space_);
    150     d_U_ = NLP_scaling()->apply_vector_scaling_d_U_NonConst(Pd_U_, ConstPtr(d_U_), d_space_);
     161    x_L_ = NLP_scaling()->apply_vector_scaling_x_LU(*Px_L_, x_L_, *x_space_);
     162    x_U_ = NLP_scaling()->apply_vector_scaling_x_LU(*Px_U_, x_U_, *x_space_);
     163    d_L_ = NLP_scaling()->apply_vector_scaling_d_LU(*Pd_L_, d_L_, *d_space_);
     164    d_U_ = NLP_scaling()->apply_vector_scaling_d_LU(*Pd_U_, d_U_, *d_space_);
    151165
    152166    // Create the iterates structures
     
    159173    v_U = d_u_space_->MakeNew();
    160174
    161     retValue = nlp_->GetStartingPoint(*x, init_x,
    162                                       *y_c, init_y_c,
    163                                       *y_d, init_y_d,
    164                                       *z_L, init_z_L,
    165                                       *z_U, init_z_U,
    166                                       *v_L, init_v_L,
    167                                       *v_U, init_v_U);
    168 
    169     if (init_x) {
    170       x = NLP_scaling()->apply_vector_scaling_x_NonConst(ConstPtr(x));
    171     }
    172     if (init_y_c) {
    173       y_c = NLP_scaling()->apply_vector_scaling_c_NonConst(ConstPtr(y_c));
    174     }
    175     if (init_y_d) {
    176       y_d = NLP_scaling()->apply_vector_scaling_d_NonConst(ConstPtr(y_d));
    177     }
    178     if (init_z_L) {
    179       z_L = NLP_scaling()->apply_vector_scaling_x_L_NonConst(Px_L_, ConstPtr(z_L), x_space_);
    180     }
    181     if (init_z_U) {
    182       z_U = NLP_scaling()->apply_vector_scaling_x_U_NonConst(Px_U_, ConstPtr(z_U), x_space_);
    183     }
    184     if (init_v_L) {
    185       v_L = NLP_scaling()->apply_vector_scaling_d_L_NonConst(Pd_L_, ConstPtr(v_L), d_space_);
    186     }
    187     if (init_v_U) {
    188       v_U = NLP_scaling()->apply_vector_scaling_d_U_NonConst(Pd_U_, ConstPtr(v_U), d_space_);
    189     }
     175    retValue = nlp_->GetStartingPoint(GetRawPtr(x), init_x,
     176                                      GetRawPtr(y_c), init_y_c,
     177                                      GetRawPtr(y_d), init_y_d,
     178                                      GetRawPtr(z_L), init_z_L,
     179                                      GetRawPtr(z_U), init_z_U);
    190180
    191181    if (!retValue) {
    192182      return false;
     183    }
     184
     185    Number obj_scal = NLP_scaling()->apply_obj_scaling(1.);
     186    if (init_x) {
     187      if (NLP_scaling()->have_x_scaling()) {
     188        x = NLP_scaling()->apply_vector_scaling_x_NonConst(ConstPtr(x));
     189      }
     190    }
     191    if (init_y_c) {
     192      if (NLP_scaling()->have_c_scaling()) {
     193        y_c = NLP_scaling()->unapply_vector_scaling_c_NonConst(ConstPtr(y_c));
     194      }
     195      if (obj_scal!=1.) {
     196        y_c->Scal(obj_scal);
     197      }
     198    }
     199    if (init_y_d) {
     200      if (NLP_scaling()->have_d_scaling()) {
     201        y_d = NLP_scaling()->unapply_vector_scaling_d_NonConst(ConstPtr(y_d));
     202      }
     203      if (obj_scal!=1.) {
     204        y_d->Scal(obj_scal);
     205      }
     206    }
     207    if (init_z_L) {
     208      if (NLP_scaling()->have_x_scaling()) {
     209        z_L = NLP_scaling()->apply_vector_scaling_x_LU_NonConst(*Px_L_, ConstPtr(z_L), *x_space_);
     210      }
     211      if (obj_scal!=1.) {
     212        z_L->Scal(obj_scal);
     213      }
     214    }
     215    if (init_z_U) {
     216      if (NLP_scaling()->have_x_scaling()) {
     217        z_U = NLP_scaling()->apply_vector_scaling_x_LU_NonConst(*Px_U_, ConstPtr(z_U), *x_space_);
     218      }
     219      if (obj_scal!=1.) {
     220        z_U->Scal(obj_scal);
     221      }
    193222    }
    194223
     
    219248    DBG_START_METH("OrigIpoptNLP::f", dbg_verbosity);
    220249    Number ret = 0.0;
    221     DBG_PRINT((1, "x.Tag = %d\n", x.GetTag()));
     250    DBG_PRINT((2, "x.Tag = %d\n", x.GetTag()));
    222251    if (!f_cache_.GetCachedResult1Dep(ret, &x)) {
    223252      f_evals_++;
    224253      SmartPtr<const Vector> unscaled_x = NLP_scaling()->unapply_vector_scaling_x(&x);
    225254      bool success = nlp_->Eval_f(*unscaled_x, ret);
     255      DBG_PRINT((1, "success = %d ret = %e\n", success, ret));
    226256      ASSERT_EXCEPTION(success && FiniteNumber(ret), Eval_Error,
    227257                       "Error evaluating the objective function");
     
    252282  }
    253283
    254 
    255284  /** Equality constraint residual */
    256285  SmartPtr<const Vector> OrigIpoptNLP::c(const Vector& x)
    257286  {
    258     SmartPtr<Vector> unscaled_c;
    259287    SmartPtr<const Vector> retValue;
    260     if (!c_cache_.GetCachedResult1Dep(retValue, &x)) {
     288    SmartPtr<const Vector> dep;
     289    if (c_space_->Dim()==0) {
     290      // We do this caching of an empty vector so that the returned
     291      // Vector has always the same tag (this might make a difference
     292      // in cases where only the constraints are supposed to change...
     293      dep = NULL;
     294      if (!c_cache_.GetCachedResult1Dep(retValue, GetRawPtr(dep))) {
     295        retValue = c_space_->MakeNew();
     296        c_cache_.AddCachedResult1Dep(retValue, GetRawPtr(dep));
     297      }
     298    }
     299    else {
     300      dep = &x;
     301    }
     302    if (!c_cache_.GetCachedResult1Dep(retValue, GetRawPtr(dep))) {
     303      SmartPtr<Vector> unscaled_c = c_space_->MakeNew();
    261304      c_evals_++;
    262       unscaled_c = c_space_->MakeNew();
    263305      SmartPtr<const Vector> unscaled_x = NLP_scaling()->unapply_vector_scaling_x(&x);
    264306      bool success = nlp_->Eval_c(*unscaled_x, *unscaled_c);
     
    266308                       Eval_Error, "Error evaluating the equality constraints");
    267309      retValue = NLP_scaling()->apply_vector_scaling_c(ConstPtr(unscaled_c));
    268       c_cache_.AddCachedResult1Dep(retValue, &x);
     310      c_cache_.AddCachedResult1Dep(retValue, GetRawPtr(dep));
    269311    }
    270312
     
    272314  }
    273315
    274 
    275316  SmartPtr<const Vector> OrigIpoptNLP::d(const Vector& x)
    276317  {
    277318    DBG_START_METH("OrigIpoptNLP::d", dbg_verbosity);
    278     SmartPtr<Vector> unscaled_d;
    279319    SmartPtr<const Vector> retValue;
    280     if (!d_cache_.GetCachedResult1Dep(retValue, &x)) {
     320    SmartPtr<const Vector> dep;
     321    if (d_space_->Dim()==0) {
     322      // We do this caching of an empty vector so that the returned
     323      // Vector has always the same tag (this might make a difference
     324      // in cases where only the constraints are supposed to change...
     325      dep = NULL;
     326      if (!d_cache_.GetCachedResult1Dep(retValue, GetRawPtr(dep))) {
     327        retValue = d_space_->MakeNew();
     328        d_cache_.AddCachedResult1Dep(retValue, GetRawPtr(dep));
     329      }
     330    }
     331    else {
     332      dep = &x;
     333    }
     334    if (!d_cache_.GetCachedResult1Dep(retValue, GetRawPtr(dep))) {
    281335      d_evals_++;
    282       unscaled_d = d_space_->MakeNew();
     336      SmartPtr<Vector> unscaled_d = d_space_->MakeNew();
    283337
    284338      DBG_PRINT_VECTOR(2, "scaled_x", x);
     
    289343                       Eval_Error, "Error evaluating the inequality constraints");
    290344      retValue = NLP_scaling()->apply_vector_scaling_d(ConstPtr(unscaled_d));
    291       d_cache_.AddCachedResult1Dep(retValue, &x);
     345      d_cache_.AddCachedResult1Dep(retValue, GetRawPtr(dep));
    292346    }
    293347
     
    297351  SmartPtr<const Matrix> OrigIpoptNLP::jac_c(const Vector& x)
    298352  {
    299     SmartPtr<Matrix> unscaled_jac_c;
    300353    SmartPtr<const Matrix> retValue;
    301     if (!jac_c_cache_.GetCachedResult1Dep(retValue, &x)) {
     354    SmartPtr<const Vector> dep;
     355    if (c_space_->Dim()==0) {
     356      // We do this caching of an empty vector so that the returned
     357      // Matrix has always the same tag
     358      dep = NULL;
     359      if (!jac_c_cache_.GetCachedResult1Dep(retValue, GetRawPtr(dep))) {
     360        retValue = jac_c_space_->MakeNew();
     361        jac_c_cache_.AddCachedResult1Dep(retValue, GetRawPtr(dep));
     362      }
     363    }
     364    else {
     365      dep = &x;
     366    }
     367    if (!jac_c_cache_.GetCachedResult1Dep(retValue, GetRawPtr(dep))) {
    302368      jac_c_evals_++;
    303       unscaled_jac_c = jac_c_space_->MakeNew();
     369      SmartPtr<Matrix> unscaled_jac_c = jac_c_space_->MakeNew();
    304370
    305371      SmartPtr<const Vector> unscaled_x = NLP_scaling()->unapply_vector_scaling_x(&x);
     
    307373      ASSERT_EXCEPTION(success, Eval_Error, "Error evaluating the jacobian of the equality constraints");
    308374      retValue = NLP_scaling()->apply_jac_c_scaling(ConstPtr(unscaled_jac_c));
    309       jac_c_cache_.AddCachedResult1Dep(retValue, &x);
     375      jac_c_cache_.AddCachedResult1Dep(retValue, GetRawPtr(dep));
    310376    }
    311377
     
    315381  SmartPtr<const Matrix> OrigIpoptNLP::jac_d(const Vector& x)
    316382  {
    317     SmartPtr<Matrix> unscaled_jac_d;
    318383    SmartPtr<const Matrix> retValue;
    319     if (!jac_d_cache_.GetCachedResult1Dep(retValue, &x)) {
     384    SmartPtr<const Vector> dep;
     385    if (d_space_->Dim()==0) {
     386      // We do this caching of an empty vector so that the returned
     387      // Matrix has always the same tag
     388      dep = NULL;
     389      if (!jac_d_cache_.GetCachedResult1Dep(retValue, GetRawPtr(dep))) {
     390        retValue = jac_d_space_->MakeNew();
     391        jac_d_cache_.AddCachedResult1Dep(retValue, GetRawPtr(dep));
     392      }
     393    }
     394    else {
     395      dep = &x;
     396    }
     397    if (!jac_d_cache_.GetCachedResult1Dep(retValue, GetRawPtr(dep))) {
    320398      jac_d_evals_++;
    321       unscaled_jac_d = jac_d_space_->MakeNew();
     399      SmartPtr<Matrix> unscaled_jac_d = jac_d_space_->MakeNew();
    322400
    323401      SmartPtr<const Vector> unscaled_x = NLP_scaling()->unapply_vector_scaling_x(&x);
     
    325403      ASSERT_EXCEPTION(success, Eval_Error, "Error evaluating the jacobian of the inequality constraints");
    326404      retValue = NLP_scaling()->apply_jac_d_scaling(ConstPtr(unscaled_jac_d));
    327       jac_d_cache_.AddCachedResult1Dep(retValue, &x);
     405      jac_d_cache_.AddCachedResult1Dep(retValue, GetRawPtr(dep));
    328406    }
    329407
     
    411489  }
    412490
     491  void OrigIpoptNLP::FinalizeSolution(SolverReturn status,
     492                                      const Vector& x, const Vector& z_L, const Vector& z_U,
     493                                      const Vector& c, const Vector& d,
     494                                      const Vector& y_c, const Vector& y_d,
     495                                      Number obj_value)
     496  {
     497    // need to submit the unscaled solution back to the nlp
     498    SmartPtr<const Vector> unscaled_x =
     499      NLP_scaling()->unapply_vector_scaling_x(&x);
     500    SmartPtr<const Vector> unscaled_c =
     501      NLP_scaling()->unapply_vector_scaling_c(&c);
     502    SmartPtr<const Vector> unscaled_d =
     503      NLP_scaling()->unapply_vector_scaling_d(&d);
     504    const Number unscaled_obj = NLP_scaling()->unapply_obj_scaling(obj_value);
     505
     506    SmartPtr<const Vector> unscaled_z_L;
     507    SmartPtr<const Vector> unscaled_z_U;
     508    SmartPtr<const Vector> unscaled_y_c;
     509    SmartPtr<const Vector> unscaled_y_d;
     510
     511    // The objective function scaling factor also appears in the constraints
     512    Number obj_unscale_factor = NLP_scaling()->unapply_obj_scaling(1.);
     513    if (obj_unscale_factor!=1.) {
     514      SmartPtr<Vector> tmp = NLP_scaling()->apply_vector_scaling_x_LU_NonConst(*Px_L_, &z_L, *x_space_);
     515      tmp->Scal(obj_unscale_factor);
     516      unscaled_z_L = ConstPtr(tmp);
     517
     518      tmp = NLP_scaling()->apply_vector_scaling_x_LU_NonConst(*Px_U_, &z_U, *x_space_);
     519      tmp->Scal(obj_unscale_factor);
     520      unscaled_z_U = ConstPtr(tmp);
     521
     522      tmp = NLP_scaling()->apply_vector_scaling_c_NonConst(&y_c);
     523      tmp->Scal(obj_unscale_factor);
     524      unscaled_y_c = ConstPtr(tmp);
     525
     526      tmp = NLP_scaling()->apply_vector_scaling_d_NonConst(&y_d);
     527      tmp->Scal(obj_unscale_factor);
     528      unscaled_y_d = ConstPtr(tmp);
     529    }
     530    else {
     531      unscaled_z_L = NLP_scaling()->apply_vector_scaling_x_LU(*Px_L_, &z_L, *x_space_);
     532      unscaled_z_U = NLP_scaling()->apply_vector_scaling_x_LU(*Px_U_, &z_U, *x_space_);
     533      unscaled_y_c = NLP_scaling()->apply_vector_scaling_c(&y_c);
     534      unscaled_y_d = NLP_scaling()->apply_vector_scaling_d(&y_d);
     535    }
     536
     537    nlp_->FinalizeSolution(status, *unscaled_x,
     538                           *unscaled_z_L, *unscaled_z_U,
     539                           *unscaled_c, *unscaled_d,
     540                           *unscaled_y_c, *unscaled_y_d,
     541                           unscaled_obj);
     542  }
     543
    413544  void OrigIpoptNLP::AdjustVariableBounds(const Vector& new_x_L, const Vector& new_x_U,
    414545                                          const Vector& new_d_L, const Vector& new_d_U)
    415546  {
    416     x_L_->Copy(new_x_L);
    417     x_U_->Copy(new_x_U);
    418     d_L_->Copy(new_d_L);
    419     d_U_->Copy(new_d_U);
     547    x_L_ = new_x_L.MakeNewCopy();
     548    x_U_ = new_x_U.MakeNewCopy();
     549    d_L_ = new_d_L.MakeNewCopy();
     550    d_U_ = new_d_U.MakeNewCopy();
    420551  }
    421552
  • trunk/Algorithm/IpOrigIpoptNLP.hpp

    r343 r424  
    1212#include "IpIpoptNLP.hpp"
    1313#include "IpException.hpp"
    14 #include"IpIpoptType.hpp"
     14#include "IpIpoptType.hpp"
    1515
    1616namespace Ipopt
     
    5656                                      bool init_z_U,
    5757                                      SmartPtr<Vector>& v_L,
    58                                       bool init_v_L,
    59                                       SmartPtr<Vector>& v_U,
    60                                       bool init_v_U
     58                                      SmartPtr<Vector>& v_U
    6159                                     );
    6260
     
    9290    virtual SmartPtr<const Vector> x_L()
    9391    {
    94       return ConstPtr(x_L_);
     92      return x_L_;
    9593    }
    9694
     
    9896    virtual SmartPtr<const Matrix> Px_L()
    9997    {
    100       return ConstPtr(Px_L_);
     98      return Px_L_;
    10199    }
    102100
     
    104102    virtual SmartPtr<const Vector> x_U()
    105103    {
    106       return ConstPtr(x_U_);
     104      return x_U_;
    107105    }
    108106
     
    110108    virtual SmartPtr<const Matrix> Px_U()
    111109    {
    112       return ConstPtr(Px_U_);
     110      return Px_U_;
    113111    }
    114112
     
    116114    virtual SmartPtr<const Vector> d_L()
    117115    {
    118       return ConstPtr(d_L_);
     116      return d_L_;
    119117    }
    120118
     
    122120    virtual SmartPtr<const Matrix> Pd_L()
    123121    {
    124       return ConstPtr(Pd_L_);
     122      return Pd_L_;
    125123    }
    126124
     
    128126    virtual SmartPtr<const Vector> d_U()
    129127    {
    130       return ConstPtr(d_U_);
     128      return d_U_;
    131129    }
    132130
     
    134132    virtual SmartPtr<const Matrix> Pd_U()
    135133    {
    136       return ConstPtr(Pd_U_);
     134      return Pd_U_;
    137135    }
    138136    //@}
     
    193191    //@}
    194192
     193    /** Solution Routines - overloaded from IpoptNLP*/
     194    //@{
     195    void FinalizeSolution(SolverReturn status,
     196                          const Vector& x, const Vector& z_L, const Vector& z_U,
     197                          const Vector& c, const Vector& d,
     198                          const Vector& y_c, const Vector& y_d,
     199                          Number obj_value);
     200    //@}
     201
    195202    /** Methods for IpoptType */
    196203    //@{
     
    200207
    201208  private:
     209    /** journalist */
     210    SmartPtr<const Journalist> jnlst_;
     211
    202212    /** Pointer to the NLP */
    203213    SmartPtr<NLP> nlp_;
    204 
    205     /** journalist */
    206     SmartPtr<const Journalist> jnlst_;
    207214
    208215    /** Necessary Vector/Matrix spaces */
     
    255262
    256263    /** Lower bounds on x */
    257     SmartPtr<Vector> x_L_;
     264    SmartPtr<const Vector> x_L_;
    258265
    259266    /** Permutation matrix (x_L_ -> x) */
    260     SmartPtr<Matrix> Px_L_;
     267    SmartPtr<const Matrix> Px_L_;
    261268
    262269    /** Upper bounds on x */
    263     SmartPtr<Vector> x_U_;
     270    SmartPtr<const Vector> x_U_;
    264271
    265272    /** Permutation matrix (x_U_ -> x */
    266     SmartPtr<Matrix> Px_U_;
     273    SmartPtr<const Matrix> Px_U_;
    267274
    268275    /** Lower bounds on d */
    269     SmartPtr<Vector> d_L_;
     276    SmartPtr<const Vector> d_L_;
    270277
    271278    /** Permutation matrix (d_L_ -> d) */
    272     SmartPtr<Matrix> Pd_L_;
     279    SmartPtr<const Matrix> Pd_L_;
    273280
    274281    /** Upper bounds on d */
    275     SmartPtr<Vector> d_U_;
     282    SmartPtr<const Vector> d_U_;
    276283
    277284    /** Permutation matrix (d_U_ -> d */
    278     SmartPtr<Matrix> Pd_U_;
     285    SmartPtr<const Matrix> Pd_U_;
    279286    //@}
    280287
     
    312319    //@}
    313320
    314     /** Flag indicating if initialization method has been called */
    315     bool initialized_;
    316 
    317321    /** @name Counters for the function evaluations */
    318322    //@{
     
    324328    Index jac_d_evals_;
    325329    Index h_evals_;
     330
     331    /** Flag indicating if initialization method has been called */
     332    bool initialized_;
    326333    //@}
    327334  };
  • trunk/Algorithm/IpOrigIterationOutput.cpp

    r311 r424  
    6464      dnrm = 0.;
    6565    }
    66     Number f = IpCq().curr_f();
     66    Number unscaled_f = IpCq().unscaled_curr_f();
    6767
    6868    // Retrieve some information set in the different parts of the algorithm
     
    8888      Jnlst().Printf(J_SUMMARY, J_MAIN,
    8989                     "%5d%c %14.7e %7.2e %7.2e %5.1f %7.2e %5s %7.2e %7.2e%c%3d %s\n",
    90                      iter, info_iter, f, inf_pr, inf_du, log10(mu), dnrm, regu_x_ptr,
     90                     iter, info_iter, unscaled_f, inf_pr, inf_du, log10(mu), dnrm, regu_x_ptr,
    9191                     alpha_dual, alpha_primal, alpha_primal_char,
    9292                     ls_count, info_string.c_str());
     
    9898    //////////////////////////////////////////////////////////////////////
    9999
    100     Jnlst().Printf(J_DETAILED, J_MAIN,
    101                    "\n**************************************************\n");
    102     Jnlst().Printf(J_DETAILED, J_MAIN,
    103                    "*** Beginning Iteration %d from the following point:",
    104                    IpData().iter_count());
    105     Jnlst().Printf(J_DETAILED, J_MAIN,
    106                    "\n**************************************************\n\n");
    107 
    108     Jnlst().Printf(J_DETAILED, J_MAIN,
    109                    "Current barrier parameter mu = %21.16e\n", IpData().curr_mu());
    110     Jnlst().Printf(J_DETAILED, J_MAIN,
    111                    "Current fraction-to-the-boundary parameter tau = %21.16e\n\n",
    112                    IpData().curr_tau());
    113 
    114100    if (Jnlst().ProduceOutput(J_DETAILED, J_MAIN)) {
     101      Jnlst().Printf(J_DETAILED, J_MAIN,
     102                     "\n**************************************************\n");
     103      Jnlst().Printf(J_DETAILED, J_MAIN,
     104                     "*** Beginning Iteration %d from the following point:",
     105                     IpData().iter_count());
     106      Jnlst().Printf(J_DETAILED, J_MAIN,
     107                     "\n**************************************************\n\n");
     108
     109      Jnlst().Printf(J_DETAILED, J_MAIN,
     110                     "Current barrier parameter mu = %21.16e\n", IpData().curr_mu());
     111      Jnlst().Printf(J_DETAILED, J_MAIN,
     112                     "Current fraction-to-the-boundary parameter tau = %21.16e\n\n",
     113                     IpData().curr_tau());
    115114      Jnlst().Printf(J_DETAILED, J_MAIN,
    116115                     "||curr_x||_inf   = %.16e\n", IpData().curr()->x()->Amax());
     
    196195    }
    197196
    198     Jnlst().Printf(J_DETAILED, J_MAIN,
    199                    "\n\n***Current NLP Values for Iteration %d:\n",
    200                    IpData().iter_count());
    201     Jnlst().Printf(J_DETAILED, J_MAIN, "Objective = %.16e\n", IpCq().curr_f());
    202     Jnlst().PrintVector(J_MOREVECTOR, J_MAIN, "grad_f", *IpCq().curr_grad_f());
    203     Jnlst().PrintVector(J_VECTOR, J_MAIN, "curr_c", *IpCq().curr_c());
    204     Jnlst().PrintVector(J_VECTOR, J_MAIN, "curr_d", *IpCq().curr_d());
    205     Jnlst().PrintVector(J_VECTOR, J_MAIN,
    206                         "curr_d - curr_s", *IpCq().curr_d_minus_s());
    207 
    208     Jnlst().PrintMatrix(J_MATRIX, J_MAIN, "jac_c", *IpCq().curr_jac_c());
    209     Jnlst().PrintMatrix(J_MATRIX, J_MAIN, "jac_d", *IpCq().curr_jac_d());
    210     Jnlst().PrintMatrix(J_MATRIX, J_MAIN, "h", *IpCq().curr_exact_hessian());
     197    if (Jnlst().ProduceOutput(J_DETAILED, J_MAIN)) {
     198      Jnlst().Printf(J_DETAILED, J_MAIN,
     199                     "\n\n***Current NLP Values for Iteration %d:\n",
     200                     IpData().iter_count());
     201      Jnlst().Printf(J_DETAILED, J_MAIN, "\n                                   (scaled)                 (unscaled)\n");
     202      Jnlst().Printf(J_DETAILED, J_MAIN, "Objective...............: %24.16e  %24.16e\n", IpCq().curr_f(), IpCq().unscaled_curr_f());
     203      Jnlst().Printf(J_DETAILED, J_MAIN, "Dual infeasibility......: %24.16e  %24.16e\n", IpCq().curr_dual_infeasibility(NORM_MAX), IpCq().unscaled_curr_dual_infeasibility(NORM_MAX));
     204      Jnlst().Printf(J_DETAILED, J_MAIN, "Constraint violation....: %24.16e  %24.16e\n", IpCq().curr_nlp_constraint_violation(NORM_MAX), IpCq().unscaled_curr_nlp_constraint_violation(NORM_MAX));
     205      Jnlst().Printf(J_DETAILED, J_MAIN, "Complementarity.........: %24.16e  %24.16e\n", IpCq().curr_complementarity(0., NORM_MAX), IpCq().unscaled_curr_complementarity(0., NORM_MAX));
     206      Jnlst().Printf(J_DETAILED, J_MAIN, "Overall NLP error.......: %24.16e  %24.16e\n\n", IpCq().curr_nlp_error(), IpCq().unscaled_curr_nlp_error());
     207    }
     208    if (Jnlst().ProduceOutput(J_VECTOR, J_MAIN)) {
     209      Jnlst().PrintVector(J_VECTOR, J_MAIN, "grad_f", *IpCq().curr_grad_f());
     210      Jnlst().PrintVector(J_VECTOR, J_MAIN, "curr_c", *IpCq().curr_c());
     211      Jnlst().PrintVector(J_VECTOR, J_MAIN, "curr_d", *IpCq().curr_d());
     212      Jnlst().PrintVector(J_VECTOR, J_MAIN,
     213                          "curr_d - curr_s", *IpCq().curr_d_minus_s());
     214    }
     215
     216    if (Jnlst().ProduceOutput(J_MATRIX, J_MAIN)) {
     217      Jnlst().PrintMatrix(J_MATRIX, J_MAIN, "jac_c", *IpCq().curr_jac_c());
     218      Jnlst().PrintMatrix(J_MATRIX, J_MAIN, "jac_d", *IpCq().curr_jac_d());
     219      Jnlst().PrintMatrix(J_MATRIX, J_MAIN, "h", *IpCq().curr_exact_hessian());
     220    }
     221
    211222    Jnlst().Printf(J_DETAILED, J_MAIN, "\n\n");
    212223
  • trunk/Algorithm/IpQualityFunctionMuOracle.cpp

    • Property svn:keywords set to Author Date Id Revision
    r329 r424  
    33// This code is published under the Common Public License.
    44//
    5 // $Id: IpOptProbingMuOracle.cpp 321 2005-06-20 21:53:55Z andreasw $
     5// $Id$
    66//
    77// Authors:  Andreas Waechter            IBM    2004-11-12
     
    4343      tmp_z_U_(NULL),
    4444      tmp_v_L_(NULL),
    45       tmp_v_U_(NULL),
    46 
    47       tmp_dual_inf_x_(NULL),
    48       tmp_dual_inf_s_(NULL)
     45      tmp_v_U_(NULL)
    4946  {
    5047    DBG_ASSERT(IsValid(pd_solver_));
     
    5653  void QualityFunctionMuOracle::RegisterOptions(SmartPtr<RegisteredOptions> roptions)
    5754  {
    58     roptions->AddLowerBoundedNumberOption("sigma_max", "???",
    59                                           0.0, true, 1e2);
    60     roptions->AddStringOption4("quality_function_norm_type", "norm to be used for the quality function", "2-norm-squared",
    61                                "1-norm", "use the 1-norm (abs sum)",
    62                                "2-norm-squared", "use the 2-norm squared (sum of squares)",
    63                                "max-norm", "use the infinity norm (max)",
    64                                "2-norm", "use 2-norm");
    65     roptions->AddStringOption2("quality_function_normalized", "set whether or not the quality function should be normalized", "no",
    66                                "no", "do not normalize the quality function",
    67                                "yes", "normalize the quality function");
    68     roptions->AddStringOption4("quality_function_centrality", "???", "none",
    69                                "none", "???",
    70                                "log", "compute the centrality as the complementarity * the log of the centrality measure",
    71                                "reciprocal", "compute the centrality as the complementarity * the reciprocal of the centrality measure",
    72                                "cubed-reciprocal", "compute the centrality as the complementarity * the reciprocal of the centrality measure cubed");
    73     roptions->AddStringOption2("quality_function_dual_inf", "???", "type-1",
    74                                "type-1", "???",
    75                                "type-2", "???");
    76     roptions->AddStringOption2("quality_function_balancing_term", "???", "none",
    77                                "none", "no balancing term",
    78                                "standard", "standard cubic balancing term");
    79     roptions->AddLowerBoundedIntegerOption("max_bisection_steps", "No Range given before...", 0, 4);
    80 
    81 
    82     roptions->AddBoundedNumberOption("bisection_tol", "tolerance for the bisection step length search",
    83                                      0.0, true, 1.0, true, 1e-3);
    84     roptions->AddStringOption2("dual_alpha_for_y", "???", "no",
    85                                "no", "???",
    86                                "yes", "???");
    87 
     55    roptions->AddLowerBoundedNumberOption(
     56      "sigma_max",
     57      "Maximal value of centering parameter.",
     58      0.0, true, 1e2);
     59    roptions->AddStringOption4(
     60      "quality_function_norm_type",
     61      "Norm used for components of the quality function.",
     62      "2-norm-squared",
     63      "1-norm", "use the 1-norm (abs sum)",
     64      "2-norm-squared", "use the 2-norm squared (sum of squares)",
     65      "max-norm", "use the infinity norm (max)",
     66      "2-norm", "use 2-norm");
     67    roptions->AddStringOption4(
     68      "quality_function_centrality",
     69      "Determines whether a penalty term for centrality is included quality function.",
     70      "none",
     71      "none", "no penalty term is added",
     72      "log", "complementarity * the log of the centrality measure",
     73      "reciprocal", "complementarity * the reciprocal of the centrality measure",
     74      "cubed-reciprocal", "complementarity * the reciprocal of the centrality measure cubed",
     75      "This determines whether a term penalizing deviation from centrality "
     76      "with respect to complementarity is added the quality function.  The "
     77      "complementarity measure here is the xi in the Loqo update rule.");
     78    roptions->AddStringOption2(
     79      "quality_function_balancing_term",
     80      "Determines whether a balancing term for centrality is included in quality function.",
     81      "none",
     82      "none", "no balancing term is added",
     83      "cubic", "Max(0,Max(dual_ing,primal_inf)-compl)^3",
     84      "This determines whether a term penalizing stuations there the "
     85      "complementality is much smaller than dual and primal "
     86      "infeasibilities is added to the quality function.");
     87    roptions->AddLowerBoundedIntegerOption(
     88      "max_bisection_steps",
     89      "Maximal number of search steps during direct search procedure determining optimal centering parameter.",
     90      0, 4);
     91    roptions->AddBoundedNumberOption(
     92      "bisection_tol",
     93      "Tolerance for the bisection search procedure determining optimal centering parameter.",
     94      0.0, true, 1.0, true,
     95      1e-3);
    8896  }
    8997
     
    97105
    98106    options.GetEnumValue("quality_function_norm_type", enum_int, prefix);
    99     quality_function_norm_ = NonmonotoneMuUpdate::NormEnum(enum_int);
    100     options.GetBoolValue("quality_function_normalized", quality_function_normalized_, prefix);
     107    quality_function_norm_ = NormEnum(enum_int);
    101108    options.GetEnumValue("quality_function_centrality", enum_int, prefix);
    102     quality_function_centrality_ = NonmonotoneMuUpdate::CentralityEnum(enum_int);
    103     options.GetEnumValue("quality_function_dual_inf", enum_int, prefix);
    104     quality_function_dual_inf_ = QualityFunctionDualInfeasibilityTypeEnum(enum_int);
     109    quality_function_centrality_ = CentralityEnum(enum_int);
    105110    options.GetEnumValue("quality_function_balancing_term", enum_int, prefix);
    106     quality_function_balancing_term_ = NonmonotoneMuUpdate::BalancingTermEnum(enum_int);
     111    quality_function_balancing_term_ = BalancingTermEnum(enum_int);
    107112    options.GetIntegerValue("max_bisection_steps", max_bisection_steps_, prefix);
    108113    options.GetNumericValue("bisection_tol", bisection_tol_, prefix);
    109 
    110     if (quality_function_normalized_) {
    111       // Set the scaling values to a negative value to indicate that
    112       // they still have to be determined
    113       dual_inf_scal_ = -1.;
    114       primal_inf_scal_ = -1.;
    115       compl_inf_scal_ = -1.;
    116     }
    117     else {
    118       dual_inf_scal_ = 1.;
    119       primal_inf_scal_ = 1.;
    120       compl_inf_scal_ = 1.;
    121     }
    122 
    123     options.GetBoolValue("dual_alpha_for_y", dual_alpha_for_y_, prefix);
    124114
    125115    return true;
     
    152142    tmp_v_L_ = IpNLP().d_L()->MakeNew();
    153143    tmp_v_U_ = IpNLP().d_U()->MakeNew();
    154 
    155     if (quality_function_dual_inf_==TYPE2) {
    156       tmp_dual_inf_x_ = IpCq().curr_grad_lag_x()->MakeNew();
    157       tmp_dual_inf_s_ = IpCq().curr_grad_lag_s()->MakeNew();
    158     }
    159144
    160145    /////////////////////////////////////
     
    200185    // First get the right hand side
    201186    SmartPtr<IteratesVector> rhs_cen = IpData().curr()->MakeNewIteratesVector(true);
    202     rhs_cen->x_NonConst()->Copy(*IpCq().grad_kappa_times_damping_x());
    203     rhs_cen->x_NonConst()->Scal(-avrg_compl);
    204     rhs_cen->s_NonConst()->Copy(*IpCq().grad_kappa_times_damping_s());
    205     rhs_cen->s_NonConst()->Scal(-avrg_compl);
     187    rhs_cen->x_NonConst()->AddOneVector(-avrg_compl,
     188                                        *IpCq().grad_kappa_times_damping_x(),
     189                                        0.);
     190    rhs_cen->s_NonConst()->AddOneVector(-avrg_compl,
     191                                        *IpCq().grad_kappa_times_damping_s(),
     192                                        0.);
    206193
    207194    rhs_cen->y_c_NonConst()->Set(0.);
     
    212199    rhs_cen->v_U_NonConst()->Set(avrg_compl);
    213200
    214     // Get space for the affine scaling step
     201    // Get space for the centering step
    215202    SmartPtr<IteratesVector> step_cen = IpData().curr()->MakeNewIteratesVector(true);
    216203
     
    247234    IpNLP().Pd_U()->TransMultVector(-1., *step_cen->s(), 0., *step_cen_s_U);
    248235
    249     // If necessary, compute some products with the constraint Jacobian
    250     SmartPtr<const Vector> jac_cT_times_step_aff_y_c;
    251     SmartPtr<const Vector> jac_dT_times_step_aff_y_d;
    252     SmartPtr<const Vector> jac_cT_times_step_cen_y_c;
    253     SmartPtr<const Vector> jac_dT_times_step_cen_y_d;
    254     if (quality_function_dual_inf_==TYPE2 && dual_alpha_for_y_) {
    255       jac_cT_times_step_aff_y_c = IpCq().curr_jac_cT_times_vec(*step_aff->y_c());
    256       jac_dT_times_step_aff_y_d = IpCq().curr_jac_dT_times_vec(*step_aff->y_d());
    257       jac_cT_times_step_cen_y_c = IpCq().curr_jac_cT_times_vec(*step_cen->y_c());
    258       jac_dT_times_step_cen_y_d = IpCq().curr_jac_dT_times_vec(*step_cen->y_d());
    259     }
    260 
    261236    Number sigma;
    262237    if (max_bisection_steps_>0) {
     
    268243      Number sigma_lo = 1e-9/avrg_compl;
    269244      sigma = PerformGoldenBisection(sigma_up, sigma_lo, bisection_tol_,
     245                                     //sigma = PerformGoldenBisectionLog(sigma_up, sigma_lo, bisection_tol_,
    270246                                     *step_aff_x_L,
    271247                                     *step_aff_x_U,
     
    287263                                     *step_cen->z_U(),
    288264                                     *step_cen->v_L(),
    289                                      *step_cen->v_U(),
    290                                      jac_cT_times_step_aff_y_c,
    291                                      jac_dT_times_step_aff_y_d,
    292                                      jac_cT_times_step_cen_y_c,
    293                                      jac_dT_times_step_cen_y_d);
     265                                     *step_cen->v_U());
    294266
    295267      if (sigma_max_ > 1. && sigma >= 1.-2*bisection_tol_) {
     
    298270        sigma_lo = sigma;
    299271        sigma = PerformGoldenBisection(sigma_up, sigma_lo, bisection_tol_,
     272                                       //sigma = PerformGoldenBisectionLog(sigma_up, sigma_lo, bisection_tol_,
    300273                                       *step_aff_x_L,
    301274                                       *step_aff_x_U,
     
    317290                                       *step_cen->z_U(),
    318291                                       *step_cen->v_L(),
    319                                        *step_cen->v_U(),
    320                                        jac_cT_times_step_aff_y_c,
    321                                        jac_dT_times_step_aff_y_d,
    322                                        jac_cT_times_step_cen_y_c,
    323                                        jac_dT_times_step_cen_y_d);
     292                                       *step_cen->v_U());
    324293      }
    325294
     
    352321                                 *step_cen->z_U(),
    353322                                 *step_cen->v_L(),
    354                                  *step_cen->v_U(),
    355                                  jac_cT_times_step_aff_y_c,
    356                                  jac_dT_times_step_aff_y_d,
    357                                  jac_cT_times_step_cen_y_c,
    358                                  jac_dT_times_step_cen_y_d);
     323                                 *step_cen->v_U());
    359324      }
    360325#endif
     
    390355                                        *step_cen->z_U(),
    391356                                        *step_cen->v_L(),
    392                                         *step_cen->v_U(),
    393                                         jac_cT_times_step_aff_y_c,
    394                                         jac_dT_times_step_aff_y_d,
    395                                         jac_cT_times_step_cen_y_c,
    396                                         jac_dT_times_step_cen_y_d);
     357                                        *step_cen->v_U());
    397358      Index l_min = (Index)(-(log(avrg_compl)-log(1e-9))/log(base))-1;
    398359      for (; l>=l_min; l--) {
     
    418379                                            *step_cen->z_U(),
    419380                                            *step_cen->v_L(),
    420                                             *step_cen->v_U(),
    421                                             jac_cT_times_step_aff_y_c,
    422                                             jac_dT_times_step_aff_y_d,
    423                                             jac_cT_times_step_cen_y_c,
    424                                             jac_dT_times_step_cen_y_d);
     381                                            *step_cen->v_U());
    425382        if (q<=q_best) {
    426383          q_best = q;
     
    469426    tmp_v_L_ = NULL;
    470427    tmp_v_U_ = NULL;
    471 
    472     if (quality_function_dual_inf_==TYPE2) {
    473       tmp_dual_inf_x_ = NULL;
    474       tmp_dual_inf_s_ = NULL;
    475     }
    476428
    477429    // DELETEME
     
    509461   const Vector& step_cen_z_U,
    510462   const Vector& step_cen_v_L,
    511    const Vector& step_cen_v_U,
    512    SmartPtr<const Vector> jac_cT_times_step_aff_y_c,
    513    SmartPtr<const Vector> jac_dT_times_step_aff_y_d,
    514    SmartPtr<const Vector> jac_cT_times_step_cen_y_c,
    515    SmartPtr<const Vector> jac_dT_times_step_cen_y_d
     463   const Vector& step_cen_v_U
    516464  )
    517465  {
     
    523471    Index n_comp = IpData().curr()->z_L()->Dim() + IpData().curr()->z_U()->Dim() +
    524472                   IpData().curr()->v_L()->Dim() + IpData().curr()->v_U()->Dim();
    525 
    526     // The scaling values have not yet been determined, compute them now
    527     if (dual_inf_scal_ < 0.) {
    528       // DELETEME
    529       assert(false && "Scaling in quality function not supported.");
    530       DBG_ASSERT(primal_inf_scal_ < 0.);
    531       DBG_ASSERT(compl_inf_scal_ < 0.);
    532 
    533       switch (quality_function_norm_) {
    534         case NonmonotoneMuUpdate::NM_NORM_1:
    535         dual_inf_scal_ = Max(1., IpCq().curr_grad_lag_x()->Asum() +
    536                              IpCq().curr_grad_lag_s()->Asum());
    537 
    538         primal_inf_scal_ = Max(1., IpCq().curr_c()->Asum() +
    539                                IpCq().curr_d_minus_s()->Asum());
    540 
    541         compl_inf_scal_ = Max(1., IpCq().curr_compl_x_L()->Asum() +
    542                               IpCq().curr_compl_x_U()->Asum() +
    543                               IpCq().curr_compl_s_L()->Asum() +
    544                               IpCq().curr_compl_s_U()->Asum());
    545 
    546         dual_inf_scal_ /= n_dual;
    547         if (n_pri>0) {
    548           primal_inf_scal_ /= n_pri;
    549         }
    550         DBG_ASSERT(n_comp>0);
    551         compl_inf_scal_ /= n_comp;
    552         break;
    553         case NonmonotoneMuUpdate::NM_NORM_2_SQUARED:
    554         dual_inf_scal_ = Max(1., pow(IpCq().curr_grad_lag_x()->Nrm2(), 2) +
    555                              pow(IpCq().curr_grad_lag_s()->Nrm2(), 2));
    556 
    557         primal_inf_scal_ = Max(1., pow(IpCq().curr_c()->Nrm2(), 2) +
    558                                pow(IpCq().curr_d_minus_s()->Nrm2(), 2));
    559 
    560         compl_inf_scal_ = Max(1., pow(IpCq().curr_compl_x_L()->Nrm2(), 2) +
    561                               pow(IpCq().curr_compl_x_U()->Nrm2(), 2) +
    562                               pow(IpCq().curr_compl_s_L()->Nrm2(), 2) +
    563                               pow(IpCq().curr_compl_s_U()->Nrm2(), 2));
    564 
    565         dual_inf_scal_ /= n_dual;
    566         if (n_pri>0) {
    567           primal_inf_scal_ /= n_pri;
    568         }
    569         DBG_ASSERT(n_comp>0);
    570         compl_inf_scal_ /= n_comp;
    571         break;
    572         case NonmonotoneMuUpdate::NM_NORM_MAX:
    573         dual_inf_scal_ = Max(1., IpCq().curr_grad_lag_x()->Amax(),
    574                              IpCq().curr_grad_lag_s()->Amax());
    575 
    576         primal_inf_scal_ = Max(1., IpCq().curr_c()->Amax(),
    577                                IpCq().curr_d_minus_s()->Amax());
    578 
    579         compl_inf_scal_ = Max(1., Max(IpCq().curr_compl_x_L()->Amax(),
    580                                       IpCq().curr_compl_x_U()->Amax(),
    581                                       IpCq().curr_compl_s_L()->Amax(),
    582                                       IpCq().curr_compl_s_U()->Amax()));
    583       }
    584     }
    585473
    586474    tmp_step_x_L_->AddTwoVectors(1., step_aff_x_L, sigma, step_cen_x_L, 0.);
     
    594482
    595483    // Compute the fraction-to-the-boundary step sizes
    596     // ToDo make sure we use the correct tau
    597     //Number tau = 0.99;
    598484    Number tau = IpData().curr_tau();
    599485    Number alpha_primal = IpCq().slack_frac_to_the_bound(tau,
     
    608494                        *tmp_step_v_L_,
    609495                        *tmp_step_v_U_);
    610 
    611     if (false) {
    612       if (alpha_dual < alpha_primal) {
    613         alpha_primal = alpha_dual;
    614       }
    615       else {
    616         alpha_dual = alpha_primal;
    617       }
    618     }
    619496
    620497    Number xi; // centrality measure
     
    655532    Number compl_inf;
    656533
    657     if (quality_function_dual_inf_==TYPE2) {
    658       tmp_dual_inf_x_->AddOneVector(1.-alpha_primal,
    659                                     *IpCq().curr_grad_lag_x(), 0.);
    660       tmp_dual_inf_s_->AddOneVector(1.-alpha_primal,
    661                                     *IpCq().curr_grad_lag_s(), 0.);
    662       IpNLP().Px_L()->MultVector(alpha_primal-alpha_dual,
    663                                  *tmp_step_z_L_, 1., *tmp_dual_inf_x_);
    664       IpNLP().Px_U()->MultVector(-alpha_primal+alpha_dual,
    665                                  *tmp_step_z_U_, 1., *tmp_dual_inf_x_);
    666       IpNLP().Pd_L()->MultVector(alpha_primal-alpha_dual,
    667                                  *tmp_step_v_L_, 1., *tmp_dual_inf_s_);
    668       IpNLP().Pd_U()->MultVector(-alpha_primal+alpha_dual,
    669                                  *tmp_step_v_U_, 1., *tmp_dual_inf_s_);
    670       if (dual_alpha_for_y_) {
    671         tmp_dual_inf_x_->AddTwoVectors(alpha_dual-alpha_primal,
    672                                        *jac_cT_times_step_aff_y_c,
    673                                        sigma*(alpha_dual-alpha_primal),
    674                                        *jac_cT_times_step_cen_y_c, 1.);
    675         tmp_dual_inf_x_->AddTwoVectors(alpha_dual-alpha_primal,
    676                                        *jac_dT_times_step_aff_y_d,
    677                                        sigma*(alpha_dual-alpha_primal),
    678                                        *jac_dT_times_step_cen_y_d, 1.);
    679         tmp_dual_inf_s_->AddTwoVectors(alpha_primal-alpha_dual,
    680                                        step_aff_y_d,
    681                                        sigma*(alpha_primal-alpha_dual),
    682                                        step_cen_y_d, 1.);
    683       }
    684     }
    685 
    686534    switch (quality_function_norm_) {
    687       case NonmonotoneMuUpdate::NM_NORM_1:
    688       if (quality_function_dual_inf_==TYPE2) {
    689         dual_inf = tmp_dual_inf_x_->Asum() + tmp_dual_inf_s_->Asum();
    690       }
    691       else {
    692         dual_inf = (1.-alpha_dual)*(IpCq().curr_grad_lag_x()->Asum() +
    693                                     IpCq().curr_grad_lag_s()->Asum());
    694       }
     535      case NM_NORM_1:
     536      dual_inf = (1.-alpha_dual)*(IpCq().curr_grad_lag_x()->Asum() +
     537                                  IpCq().curr_grad_lag_s()->Asum());
    695538
    696539      primal_inf = (1.-alpha_primal)*(IpCq().curr_c()->Asum() +
     
    707550      compl_inf /= n_comp;
    708551      break;
    709       case NonmonotoneMuUpdate::NM_NORM_2_SQUARED:
    710       if (quality_function_dual_inf_==TYPE2) {
    711         dual_inf = pow(tmp_dual_inf_x_->Nrm2(), 2) + pow(tmp_dual_inf_s_->Nrm2(), 2);
    712       }
    713       else {
    714         dual_inf =
    715           pow(1.-alpha_dual, 2)*(pow(IpCq().curr_grad_lag_x()->Nrm2(), 2) +
    716                                  pow(IpCq().curr_grad_lag_s()->Nrm2(), 2));
    717       }
     552      case NM_NORM_2_SQUARED:
     553      dual_inf =
     554        pow(1.-alpha_dual, 2)*(pow(IpCq().curr_grad_lag_x()->Nrm2(), 2) +
     555                               pow(IpCq().curr_grad_lag_s()->Nrm2(), 2));
    718556      primal_inf =
    719557        pow(1.-alpha_primal, 2)*(pow(IpCq().curr_c()->Nrm2(), 2) +
     
    730568      compl_inf /= n_comp;
    731569      break;
    732       case NonmonotoneMuUpdate::NM_NORM_MAX:
    733       if (quality_function_dual_inf_==TYPE2) {
    734         dual_inf = Max(tmp_dual_inf_x_->Amax(), tmp_dual_inf_s_->Amax());
    735       }
    736       else {
    737         dual_inf =
    738           (1.-alpha_dual)*Max(IpCq().curr_grad_lag_x()->Amax(),
    739                               IpCq().curr_grad_lag_s()->Amax());
    740       }
     570      case NM_NORM_MAX:
     571      dual_inf =
     572        (1.-alpha_dual)*Max(IpCq().curr_grad_lag_x()->Amax(),
     573                            IpCq().curr_grad_lag_s()->Amax());
    741574      primal_inf =
    742         (1.-alpha_primal, 2)*Max(IpCq().curr_c()->Amax(),
    743                                  IpCq().curr_d_minus_s()->Amax());
     575        (1.-alpha_primal)*Max(IpCq().curr_c()->Amax(),
     576                              IpCq().curr_d_minus_s()->Amax());
    744577      compl_inf =
    745578        Max(tmp_slack_x_L_->Amax(), tmp_slack_x_U_->Amax(),
    746579            tmp_slack_s_L_->Amax(), tmp_slack_s_U_->Amax());
    747580      break;
    748       case NonmonotoneMuUpdate::NM_NORM_2:
    749       if (quality_function_dual_inf_==TYPE2) {
    750         dual_inf = sqrt(pow(tmp_dual_inf_x_->Nrm2(), 2) + pow(tmp_dual_inf_s_->Nrm2(), 2));
    751       }
    752       else {
    753         dual_inf =
    754           (1.-alpha_dual)*sqrt(pow(IpCq().curr_grad_lag_x()->Nrm2(), 2) +
    755                                pow(IpCq().curr_grad_lag_s()->Nrm2(), 2));
    756       }
     581      case NM_NORM_2:
     582      dual_inf =
     583        (1.-alpha_dual)*sqrt(pow(IpCq().curr_grad_lag_x()->Nrm2(), 2) +
     584                             pow(IpCq().curr_grad_lag_s()->Nrm2(), 2));
    757585      primal_inf =
    758         (1.-alpha_primal, 2)*sqrt(pow(IpCq().curr_c()->Nrm2(), 2) +
    759                                   pow(IpCq().curr_d_minus_s()->Nrm2(), 2));
     586        (1.-alpha_primal)*sqrt(pow(IpCq().curr_c()->Nrm2(), 2) +
     587                               pow(IpCq().curr_d_minus_s()->Nrm2(), 2));
    760588      compl_inf =
    761589        sqrt(pow(tmp_slack_x_L_->Nrm2(), 2) + pow(tmp_slack_x_U_->Nrm2(), 2) +
     
    773601    }
    774602
    775     // Scale the quantities
    776     dual_inf /= dual_inf_scal_;
    777     primal_inf /= primal_inf_scal_;
    778     compl_inf /= compl_inf_scal_;
    779 
    780603    Number quality_function = dual_inf + primal_inf + compl_inf;
    781604
    782605    switch (quality_function_centrality_) {
    783       case NonmonotoneMuUpdate::CEN_NONE:
     606      case CEN_NONE:
    784607      //Nothing
    785608      break;
    786       case NonmonotoneMuUpdate::CEN_LOG:
     609      case CEN_LOG:
    787610      quality_function -= compl_inf*log(xi);
    788611      break;
    789       case NonmonotoneMuUpdate::CEN_RECIPROCAL:
     612      case CEN_RECIPROCAL:
    790613      quality_function += compl_inf/xi;
    791       case NonmonotoneMuUpdate::CEN_CUBED_RECIPROCAL:
     614      case CEN_CUBED_RECIPROCAL:
    792615      quality_function += compl_inf/pow(xi,3);
    793616      break;
     
    797620
    798621    switch (quality_function_balancing_term_) {
    799       case NonmonotoneMuUpdate::BT_NONE:
     622      case BT_NONE:
    800623      //Nothing
    801624      break;
    802       case NonmonotoneMuUpdate::BT_CUBIC:
     625      case BT_CUBIC:
    803626      quality_function += pow(Max(0., Max(dual_inf,primal_inf)-compl_inf),3);
    804627      break;
     
    839662   const Vector& step_cen_z_U,
    840663   const Vector& step_cen_v_L,
    841    const Vector& step_cen_v_U,
    842    SmartPtr<const Vector> jac_cT_times_step_aff_y_c,
    843    SmartPtr<const Vector> jac_dT_times_step_aff_y_d,
    844    SmartPtr<const Vector> jac_cT_times_step_cen_y_c,
    845    SmartPtr<const Vector> jac_dT_times_step_cen_y_d
     664   const Vector& step_cen_v_U
    846665  )
    847666  {
     
    873692                                            step_cen_z_U,
    874693                                            step_cen_v_L,
    875                                             step_cen_v_U,
    876                                             jac_cT_times_step_aff_y_c,
    877                                             jac_dT_times_step_aff_y_d,
    878                                             jac_cT_times_step_cen_y_c,
    879                                             jac_dT_times_step_cen_y_d);
     694                                            step_cen_v_U);
    880695    Number qmid2 = CalculateQualityFunction(sigma_mid2,
    881696                                            step_aff_x_L,
     
    898713                                            step_cen_z_U,
    899714                                            step_cen_v_L,
    900                                             step_cen_v_U,
    901                                             jac_cT_times_step_aff_y_c,
    902                                             jac_dT_times_step_aff_y_d,
    903                                             jac_cT_times_step_cen_y_c,
    904                                             jac_dT_times_step_cen_y_d);
     715                                            step_cen_v_U);
    905716
    906717    Index nbisections = 0;
    907718    while ((sigma_up-sigma_lo)>=tol*sigma_up && nbisections<max_bisection_steps_) {
    908       //      printf("lo = %e mid1 = %e mid2 = %e up = %e\n",sigma_lo,sigma_mid1,sigma_mid2,sigma_up);
    909719      nbisections++;
    910720      if (qmid1 > qmid2) {
     
    933743                                         step_cen_z_U,
    934744                                         step_cen_v_L,
    935                                          step_cen_v_U,
    936                                          jac_cT_times_step_aff_y_c,
    937                                          jac_dT_times_step_aff_y_d,
    938                                          jac_cT_times_step_cen_y_c,
    939                                          jac_dT_times_step_cen_y_d);
     745                                         step_cen_v_U);
    940746      }
    941747      else {
     
    964770                                         step_cen_z_U,
    965771                                         step_cen_v_L,
    966                                          step_cen_v_U,
    967                                          jac_cT_times_step_aff_y_c,
    968                                          jac_dT_times_step_aff_y_d,
    969                                          jac_cT_times_step_cen_y_c,
    970                                          jac_dT_times_step_cen_y_d);
     772                                         step_cen_v_U);
    971773      }
    972774    }
     
    1002804                                             step_cen_z_U,
    1003805                                             step_cen_v_L,
    1004                                              step_cen_v_U,
    1005                                              jac_cT_times_step_aff_y_c,
    1006                                              jac_dT_times_step_aff_y_d,
    1007                                              jac_cT_times_step_cen_y_c,
    1008                                              jac_dT_times_step_cen_y_d);
     806                                             step_cen_v_U);
    1009807      if (qtmp < q) {
    1010808        sigma = sigma_up;
     
    1033831                                             step_cen_z_U,
    1034832                                             step_cen_v_L,
    1035                                              step_cen_v_U,
    1036                                              jac_cT_times_step_aff_y_c,
    1037                                              jac_dT_times_step_aff_y_d,
    1038                                              jac_cT_times_step_cen_y_c,
    1039                                              jac_dT_times_step_cen_y_d);
     833                                             step_cen_v_U);
    1040834      if (qtmp < q) {
    1041835        sigma = sigma_lo;
     
    1047841  }
    1048842
     843  /* AW: Tried search in the log space, but that was even worse than
     844     search in unscaled space */
     845  /*
     846  Number
     847  QualityFunctionMuOracle::PerformGoldenBisectionLog
     848  (Number sigma_up,
     849   Number sigma_lo,
     850   Number tol,
     851   const Vector& step_aff_x_L,
     852   const Vector& step_aff_x_U,
     853   const Vector& step_aff_s_L,
     854   const Vector& step_aff_s_U,
     855   const Vector& step_aff_y_c,
     856   const Vector& step_aff_y_d,
     857   const Vector& step_aff_z_L,
     858   const Vector& step_aff_z_U,
     859   const Vector& step_aff_v_L,
     860   const Vector& step_aff_v_U,
     861   const Vector& step_cen_x_L,
     862   const Vector& step_cen_x_U,
     863   const Vector& step_cen_s_L,
     864   const Vector& step_cen_s_U,
     865   const Vector& step_cen_y_c,
     866   const Vector& step_cen_y_d,
     867   const Vector& step_cen_z_L,
     868   const Vector& step_cen_z_U,
     869   const Vector& step_cen_v_L,
     870   const Vector& step_cen_v_U
     871  )
     872  {
     873    Number log_sigma;
     874    Number log_sigma_up = log(sigma_up);
     875    Number log_sigma_lo = log(sigma_lo);
     876
     877    Number log_sigma_up_in = log_sigma_up;
     878    Number log_sigma_lo_in = log_sigma_lo;
     879    Number gfac = (3.-sqrt(5.))/2.;
     880    Number log_sigma_mid1 = log_sigma_lo + gfac*(log_sigma_up-log_sigma_lo);
     881    Number log_sigma_mid2 = log_sigma_lo + (1.-gfac)*(log_sigma_up-log_sigma_lo);
     882
     883    Number qmid1 = CalculateQualityFunction(exp(log_sigma_mid1),
     884                                            step_aff_x_L,
     885                                            step_aff_x_U,
     886                                            step_aff_s_L,
     887                                            step_aff_s_U,
     888                                            step_aff_y_c,
     889                                            step_aff_y_d,
     890                                            step_aff_z_L,
     891                                            step_aff_z_U,
     892                                            step_aff_v_L,
     893                                            step_aff_v_U,
     894                                            step_cen_x_L,
     895                                            step_cen_x_U,
     896                                            step_cen_s_L,
     897                                            step_cen_s_U,
     898                                            step_cen_y_c,
     899                                            step_cen_y_d,
     900                                            step_cen_z_L,
     901                                            step_cen_z_U,
     902                                            step_cen_v_L,
     903                                            step_cen_v_U);
     904    Number qmid2 = CalculateQualityFunction(exp(log_sigma_mid2),
     905                                            step_aff_x_L,
     906                                            step_aff_x_U,