Changeset 529


Ignore:
Timestamp:
Sep 29, 2005 5:12:38 PM (15 years ago)
Author:
andreasw
Message:
  • in Vector class, copy cached values in Copy and update them in Scal
  • perform iterative refinement only once for adaptive strategy
  • fix bugs in PDPerturbationHandler
  • avoid some overhead in CalculateQualityFunction?
  • minor changes in timing
Location:
branches/dev
Files:
17 edited

Legend:

Unmodified
Added
Removed
  • branches/dev/Algorithm/IpAdaptiveMuUpdate.cpp

    r526 r529  
    5252      "mu_max",
    5353      "Maximum value for barrier parameter.",
    54       0.0, true, 1e10,
     54      0.0, true, 1e5,
    5555      "This option specifies an upper bound on the barrier parameter in the "
    5656      "adaptive mu selection mode.");
  • branches/dev/Algorithm/IpFilterLineSearch.cpp

    r523 r529  
    12891289    }
    12901290
    1291     IpData().TimingStats().TryCorrector.Start();
     1291    IpData().TimingStats().TryCorrector().Start();
    12921292
    12931293    bool accept = false;
     
    14411441      Jnlst().Printf(J_DETAILED, J_LINE_SEARCH,
    14421442                     "Rejecting corrector step, because trial complementarity is too large.\n" );
    1443       IpData().TimingStats().TryCorrector.End();
     1443      IpData().TimingStats().TryCorrector().End();
    14441444      return false;
    14451445    }
     
    14741474    }
    14751475
    1476     IpData().TimingStats().TryCorrector.End();
     1476    IpData().TimingStats().TryCorrector().End();
    14771477    return accept;
    14781478  }
  • branches/dev/Algorithm/IpIpoptAlg.cpp

    r526 r529  
    145145
    146146    // Start measuring CPU time
    147     IpData().TimingStats().OverallAlgorithm.Start();
     147    IpData().TimingStats().OverallAlgorithm().Start();
    148148
    149149    if (!message_printed) {
     
    151151    }
    152152
    153     IpData().TimingStats().InitializeIterates.Start();
    154     // Initialize the iterates
    155     InitializeIterates();
    156     IpData().TimingStats().InitializeIterates.End();
    157 
    158     if (!skip_print_problem_stats_) {
    159       IpData().TimingStats().PrintProblemStatistics.Start();
    160       PrintProblemStatistics();
    161       IpData().TimingStats().PrintProblemStatistics.End();
    162     }
    163 
    164     IpData().TimingStats().CheckConvergence.Start();
    165     ConvergenceCheck::ConvergenceStatus conv_status
    166     = conv_check_->CheckConvergence();
    167     IpData().TimingStats().CheckConvergence.End();
    168 
    169153    try {
     154      IpData().TimingStats().InitializeIterates().Start();
     155      // Initialize the iterates
     156      InitializeIterates();
     157      IpData().TimingStats().InitializeIterates().End();
     158
     159      if (!skip_print_problem_stats_) {
     160        IpData().TimingStats().PrintProblemStatistics().Start();
     161        PrintProblemStatistics();
     162        IpData().TimingStats().PrintProblemStatistics().End();
     163      }
     164
     165      IpData().TimingStats().CheckConvergence().Start();
     166      ConvergenceCheck::ConvergenceStatus conv_status
     167      = conv_check_->CheckConvergence();
     168      IpData().TimingStats().CheckConvergence().End();
     169
    170170      // main loop
    171171      while (conv_status == ConvergenceCheck::CONTINUE) {
    172172        // Set the Hessian Matrix
    173         IpData().TimingStats().ActualizeHessian.Start();
     173        IpData().TimingStats().ActualizeHessian().Start();
    174174        ActualizeHessian();
    175         IpData().TimingStats().ActualizeHessian.End();
     175        IpData().TimingStats().ActualizeHessian().End();
    176176
    177177        // do all the output for this iteration
    178         IpData().TimingStats().OutputIteration.Start();
     178        IpData().TimingStats().OutputIteration().Start();
    179179        OutputIteration();
    180180        IpData().ResetInfo();
    181         IpData().TimingStats().OutputIteration.End();
     181        IpData().TimingStats().OutputIteration().End();
    182182
    183183        // update the barrier parameter if necessary
    184         IpData().TimingStats().UpdateBarrierParameter.Start();
     184        IpData().TimingStats().UpdateBarrierParameter().Start();
    185185        UpdateBarrierParameter();
    186         IpData().TimingStats().UpdateBarrierParameter.End();
     186        IpData().TimingStats().UpdateBarrierParameter().End();
    187187
    188188        // solve the primal-dual system to get the full step
    189         IpData().TimingStats().ComputeSearchDirection.Start();
     189        IpData().TimingStats().ComputeSearchDirection().Start();
    190190        ComputeSearchDirection();
    191         IpData().TimingStats().ComputeSearchDirection.End();
     191        IpData().TimingStats().ComputeSearchDirection().End();
    192192
    193193        // Compute the new iterate
    194         IpData().TimingStats().ComputeAcceptableTrialPoint.Start();
     194        IpData().TimingStats().ComputeAcceptableTrialPoint().Start();
    195195        ComputeAcceptableTrialPoint();
    196         IpData().TimingStats().ComputeAcceptableTrialPoint.End();
     196        IpData().TimingStats().ComputeAcceptableTrialPoint().End();
    197197
    198198        // Accept the new iterate
    199         IpData().TimingStats().AcceptTrialPoint.Start();
     199        IpData().TimingStats().AcceptTrialPoint().Start();
    200200        AcceptTrialPoint();
    201         IpData().TimingStats().AcceptTrialPoint.End();
     201        IpData().TimingStats().AcceptTrialPoint().End();
    202202
    203203        IpData().Set_iter_count(IpData().iter_count()+1);
    204204
    205         IpData().TimingStats().CheckConvergence.Start();
     205        IpData().TimingStats().CheckConvergence().Start();
    206206        conv_status  = conv_check_->CheckConvergence();
    207         IpData().TimingStats().CheckConvergence.End();
    208       }
    209 
    210       IpData().TimingStats().OutputIteration.Start();
     207        IpData().TimingStats().CheckConvergence().End();
     208      }
     209
     210      IpData().TimingStats().OutputIteration().Start();
    211211      OutputIteration();
    212       IpData().TimingStats().OutputIteration.End();
    213 
    214       IpData().TimingStats().OverallAlgorithm.End();
     212      IpData().TimingStats().OutputIteration().End();
     213
     214      IpData().TimingStats().OverallAlgorithm().End();
    215215
    216216      if (conv_status == ConvergenceCheck::CONVERGED) {
     
    226226    catch(TINY_STEP_DETECTED& exc) {
    227227      exc.ReportException(Jnlst(), J_DETAILED);
    228       IpData().TimingStats().UpdateBarrierParameter.End();
    229       IpData().TimingStats().OverallAlgorithm.End();
     228      IpData().TimingStats().UpdateBarrierParameter().End();
     229      IpData().TimingStats().OverallAlgorithm().End();
    230230      return STOP_AT_TINY_STEP;
    231231    }
    232232    catch(ACCEPTABLE_POINT_REACHED& exc) {
    233233      exc.ReportException(Jnlst(), J_DETAILED);
    234       IpData().TimingStats().ComputeAcceptableTrialPoint.End();
    235       IpData().TimingStats().OverallAlgorithm.End();
     234      IpData().TimingStats().ComputeAcceptableTrialPoint().End();
     235      IpData().TimingStats().OverallAlgorithm().End();
    236236      return STOP_AT_ACCEPTABLE_POINT;
    237237    }
    238238    catch(LOCALLY_INFEASIBLE& exc) {
    239239      exc.ReportException(Jnlst(), J_DETAILED);
    240       IpData().TimingStats().ComputeAcceptableTrialPoint.EndIfStarted();
    241       IpData().TimingStats().CheckConvergence.EndIfStarted();
    242       IpData().TimingStats().OverallAlgorithm.End();
     240      IpData().TimingStats().ComputeAcceptableTrialPoint().EndIfStarted();
     241      IpData().TimingStats().CheckConvergence().EndIfStarted();
     242      IpData().TimingStats().OverallAlgorithm().End();
    243243      return LOCAL_INFEASIBILITY;
    244244    }
    245245    catch(RESTORATION_FAILED& exc) {
    246246      exc.ReportException(Jnlst(), J_DETAILED);
    247       IpData().TimingStats().ComputeAcceptableTrialPoint.End();
    248       IpData().TimingStats().OverallAlgorithm.End();
     247      IpData().TimingStats().ComputeAcceptableTrialPoint().End();
     248      IpData().TimingStats().OverallAlgorithm().End();
    249249      return RESTORATION_FAILURE;
    250250    }
    251251    catch(FEASIBILITY_PROBLEM_SOLVED& exc) {
    252252      exc.ReportException(Jnlst(), J_DETAILED);
    253       IpData().TimingStats().ComputeAcceptableTrialPoint.End();
    254       IpData().TimingStats().OverallAlgorithm.End();
     253      IpData().TimingStats().ComputeAcceptableTrialPoint().End();
     254      IpData().TimingStats().OverallAlgorithm().End();
    255255      return SUCCESS;
    256256    }
    257257    catch(INTERNAL_ABORT& exc) {
    258258      exc.ReportException(Jnlst());
    259       IpData().TimingStats().OverallAlgorithm.End();
     259      IpData().TimingStats().OverallAlgorithm().End();
    260260      return INTERNAL_ERROR;
    261261    }
     
    263263    DBG_ASSERT(false && "Unknown return code in the algorithm");
    264264
    265     IpData().TimingStats().OverallAlgorithm.End();
     265    IpData().TimingStats().OverallAlgorithm().End();
    266266    return INTERNAL_ERROR;
    267267  }
     
    297297                   "\n**************************************************\n\n");
    298298
     299    bool improve_solution = false;
    299300    if (IpData().HaveDeltas()) {
     301      /*
    300302      Jnlst().Printf(J_DETAILED, J_MAIN,
    301303                     "No need to compute search direction - it has already been computed.\n");
    302304      return;
     305      */
     306      improve_solution = true;
    303307    }
    304308
     
    316320
    317321    // Get space for the search direction
    318     SmartPtr<IteratesVector> delta = IpData().curr()->MakeNewIteratesVector(true);
    319 
    320     pd_solver_->Solve(-1.0, 0.0, *rhs, *delta);
     322    SmartPtr<IteratesVector> delta =
     323      IpData().curr()->MakeNewIteratesVector(true);
     324
     325    if (improve_solution) {
     326      // We can probably avoid copying and scaling...
     327      delta->AddOneVector(-1., *IpData().delta(), 0.);
     328    }
     329
     330    bool allow_inexact = false;
     331    pd_solver_->Solve(-1.0, 0.0, *rhs, *delta, allow_inexact,
     332                      improve_solution);
    321333
    322334    // Store the search directions in the IpData object
  • branches/dev/Algorithm/IpOrigIpoptNLP.cpp

    r526 r529  
    649649
    650650    jnlst.Printf(level, category,
    651                  "Function Evaluations................: %10.2f\n",
     651                 "Function Evaluations................: %10.3f\n",
    652652                 f_eval_time_.TotalTime()+
    653653                 c_eval_time_.TotalTime()+
     
    657657                 h_eval_time_.TotalTime());
    658658    jnlst.Printf(level, category,
    659                  " Objective function.................: %10.2f\n",
     659                 " Objective function.................: %10.3f\n",
    660660                 f_eval_time_.TotalTime());
    661661    jnlst.Printf(level, category,
    662                  " Equality constraints...............: %10.2f\n",
     662                 " Equality constraints...............: %10.3f\n",
    663663                 c_eval_time_.TotalTime());
    664664    jnlst.Printf(level, category,
    665                  " Inequality constraints.............: %10.2f\n",
     665                 " Inequality constraints.............: %10.3f\n",
    666666                 d_eval_time_.TotalTime());
    667667    jnlst.Printf(level, category,
    668                  " Equality constraint Jacobian.......: %10.2f\n",
     668                 " Equality constraint Jacobian.......: %10.3f\n",
    669669                 jac_c_eval_time_.TotalTime());
    670670    jnlst.Printf(level, category,
    671                  " Inequality constraint Jacobian.....: %10.2f\n",
     671                 " Inequality constraint Jacobian.....: %10.3f\n",
    672672                 jac_d_eval_time_.TotalTime());
    673673    jnlst.Printf(level, category,
    674                  " Lagrangian Hessian.................: %10.2f\n",
     674                 " Lagrangian Hessian.................: %10.3f\n",
    675675                 h_eval_time_.TotalTime());
    676676  }
  • branches/dev/Algorithm/IpPDFullSpaceSolver.cpp

    r524 r529  
    106106                                const IteratesVector& rhs,
    107107                                IteratesVector& res,
    108                                 bool allow_inexact)
     108                                bool allow_inexact,
     109                                bool improve_solution /* = false */)
    109110  {
    110111    DBG_START_METH("PDFullSpaceSolver::Solve",dbg_verbosity);
     112    DBG_ASSERT(!allow_inexact || !improve_solution);
     113    DBG_ASSERT(!improve_solution || beta==0.);
    111114
    112115    // Timing of PDSystem solver starts here
    113     IpData().TimingStats().PDSystemSolverTotal.Start();
     116    IpData().TimingStats().PDSystemSolverTotal().Start();
    114117
    115118    DBG_PRINT_VECTOR(2, "rhs_x", *rhs.x());
     
    121124    DBG_PRINT_VECTOR(2, "rhs_vL", *rhs.v_L());
    122125    DBG_PRINT_VECTOR(2, "rhs_vU", *rhs.v_U());
     126    DBG_PRINT_VECTOR(2, "res_x in", *res.x());
     127    DBG_PRINT_VECTOR(2, "res_s in", *res.s());
     128    DBG_PRINT_VECTOR(2, "res_c in", *res.y_c());
     129    DBG_PRINT_VECTOR(2, "res_d in", *res.y_d());
     130    DBG_PRINT_VECTOR(2, "res_zL in", *res.z_L());
     131    DBG_PRINT_VECTOR(2, "res_zU in", *res.z_U());
     132    DBG_PRINT_VECTOR(2, "res_vL in", *res.v_L());
     133    DBG_PRINT_VECTOR(2, "res_vU in", *res.v_U());
    123134
    124135    // if beta is nonzero, keep a copy of the incoming values in res_ */
     
    167178    while (!done) {
    168179
    169       bool solve_retval =
    170         SolveOnce(resolve_with_better_quality, pretend_singular,
    171                   *W, *J_c, *J_d, *Px_L, *Px_U, *Pd_L, *Pd_U, *z_L, *z_U,
    172                   *v_L, *v_U, *slack_x_L, *slack_x_U, *slack_s_L, *slack_s_U,
    173                   *sigma_x, *sigma_s, 1., 0., rhs, res);
    174       resolve_with_better_quality = false;
    175       pretend_singular = false;
     180      // if improve_solution is true, we are given already a solution
     181      // from the calling function, so we can skip the first solve
     182      bool solve_retval = true;
     183      if (!improve_solution) {
     184        solve_retval =
     185          SolveOnce(resolve_with_better_quality, pretend_singular,
     186                    *W, *J_c, *J_d, *Px_L, *Px_U, *Pd_L, *Pd_U, *z_L, *z_U,
     187                    *v_L, *v_U, *slack_x_L, *slack_x_U, *slack_s_L, *slack_s_U,
     188                    *sigma_x, *sigma_s, 1., 0., rhs, res);
     189        resolve_with_better_quality = false;
     190        pretend_singular = false;
     191      }
     192      improve_solution = false;
    176193
    177194      if (!solve_retval) {
     
    181198        // might want to use a more explicit cue later.
    182199        res.Set(0.0);
    183         IpData().TimingStats().PDSystemSolverTotal.End();
     200        IpData().TimingStats().PDSystemSolverTotal().End();
    184201        return;
    185202      }
    186203
     204      if (allow_inexact) {
     205        // no safety checks required
     206        break;
     207      }
     208
    187209      // Get space for the residual
    188210      SmartPtr<IteratesVector> resid = res.MakeNewIteratesVector(true);
    189211
     212      // ToDo don't to that after max refinement?
    190213      ComputeResiduals(*W, *J_c, *J_d, *Px_L, *Px_U, *Pd_L, *Pd_U,
    191214                       *z_L, *z_U, *v_L, *v_U, *slack_x_L, *slack_x_U,
     
    212235                    *v_L, *v_U, *slack_x_L, *slack_x_U, *slack_s_L, *slack_s_U,
    213236                    *sigma_x, *sigma_s, -1., 1., *resid, res);
    214         DBG_ASSERT(solve_retval);
     237        ASSERT_EXCEPTION(solve_retval, INTERNAL_ABORT,
     238                         "SolveOnce returns false.");
    215239
    216240        ComputeResiduals(*W, *J_c, *J_d, *Px_L, *Px_U, *Pd_L, *Pd_U,
     
    312336    DBG_PRINT_VECTOR(2, "res_vU", *res.v_U());
    313337
    314     IpData().TimingStats().PDSystemSolverTotal.End();
     338    IpData().TimingStats().PDSystemSolverTotal().End();
    315339  }
    316340
     
    351375    DBG_START_METH("PDFullSpaceSolver::SolveOnce",dbg_verbosity);
    352376
    353     IpData().TimingStats().PDSystemSolverSolveOnce.Start();
     377    IpData().TimingStats().PDSystemSolverSolveOnce().Start();
    354378
    355379    // Compute the right hand side for the augmented system formulation
     
    416440                                    false, 0);
    417441      if (retval!=SYMSOLVER_SUCCESS) {
    418         IpData().TimingStats().PDSystemSolverSolveOnce.End();
     442        IpData().TimingStats().PDSystemSolverSolveOnce().End();
    419443        return false;
    420444      }
     
    462486          // Get new perturbation factors from the perturbation
    463487          // handlers for the singular case
    464           perturbHandler_->PerturbForSingularity(delta_x, delta_s,
    465                                                  delta_c, delta_d);
     488          bool pert_return =
     489                             perturbHandler_->PerturbForSingularity(delta_x, delta_s,
     490                                                                    delta_c, delta_d);
     491          if (!pert_return) {
     492            Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA,
     493                           "PerturbForSingularity can't be done\n");
     494            IpData().TimingStats().PDSystemSolverSolveOnce().End();
     495            return false;
     496          }
    466497        }
    467498        else if (retval==SYMSOLVER_WRONG_INERTIA &&
     
    489520          }
    490521          if (assume_singular) {
    491             perturbHandler_->PerturbForSingularity(delta_x, delta_s,
    492                                                    delta_c, delta_d);
     522            bool pert_return =
     523                               perturbHandler_->PerturbForSingularity(delta_x, delta_s,
     524                                                                      delta_c, delta_d);
     525            if (!pert_return) {
     526              Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA,
     527                             "PerturbForSingularity can't be done for assume singular.\n");
     528              IpData().TimingStats().PDSystemSolverSolveOnce().End();
     529              return false;
     530            }
    493531            IpData().Append_info_string("a");
    494532          }
     
    498536          // Get new perturbation factors from the perturbation
    499537          // handlers for the case of wrong inertia
    500           perturbHandler_->PerturbForWrongInertia(delta_x, delta_s,
    501                                                   delta_c, delta_d);
     538          bool pert_return =
     539                             perturbHandler_->PerturbForWrongInertia(delta_x, delta_s,
     540                                                                     delta_c, delta_d);
     541          if (!pert_return) {
     542            Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA,
     543                           "PerturbForWrongInertia can't be done for assume singular.\n");
     544            IpData().TimingStats().PDSystemSolverSolveOnce().End();
     545            return false;
     546          }
    502547        }
    503548      } // while (retval!=SYMSOLVER_SUCCESS && !fail) {
     
    521566    res.AddOneVector(alpha, *sol, beta);
    522567
    523     IpData().TimingStats().PDSystemSolverSolveOnce.End();
     568    IpData().TimingStats().PDSystemSolverSolveOnce().End();
    524569
    525570    return true;
     
    552597    DBG_START_METH("PDFullSpaceSolver::ComputeResiduals", dbg_verbosity);
    553598
     599    DBG_PRINT_VECTOR(2, "res", res);
     600    IpData().TimingStats().ComputeResiduals().Start();
     601
    554602    // Get the current sizes of the perturbation factors
    555603    Number delta_x;
     
    619667    tmp->ElementWiseMultiply(v_U);
    620668    resid.v_U_NonConst()->AddTwoVectors(-1., *tmp, -1., *rhs.v_U(), 1.);
     669
     670    DBG_PRINT_VECTOR(2, "resid", resid);
    621671
    622672    if (Jnlst().ProduceOutput(J_MOREVECTOR, J_LINEAR_ALGEBRA)) {
     
    642692                     "max-norm resid_vU %e\n", resid.v_U()->Amax());
    643693    }
     694    IpData().TimingStats().ComputeResiduals().End();
    644695  }
    645696
  • branches/dev/Algorithm/IpPDFullSpaceSolver.hpp

    r501 r529  
    5555                       const IteratesVector& rhs,
    5656                       IteratesVector& res,
    57                        bool allow_inexact=false);
     57                       bool allow_inexact=false,
     58                       bool improve_solution=false);
    5859
    5960    /** Methods for IpoptType */
  • branches/dev/Algorithm/IpPDPerturbationHandler.cpp

    r517 r529  
    3434      "Maximum value of regularization parameter for handling negative curvature.",
    3535      0, true,
    36       1e40,
     36      1e40, // ToDo make this 1e20 ?
    3737      "In order to guarantee that the search directions are indeed proper "
    3838      "descent directions, Ipopt requires that the inertia of the "
     
    176176
    177177    if (hess_degenerate_ == DEGENERATE) {
     178      delta_x_curr_ = 0.;
     179      delta_s_curr_ = 0.;
    178180      bool retval = get_deltas_for_wrong_inertia(delta_x, delta_s,
    179181                    delta_c, delta_d);
    180       Jnlst().Printf(J_ERROR, J_MAIN, "Cannot determine appropriate modification of Hessian matrix.\nThis can happen if there are NaN or Inf numbers in the user-provided Hessian.\n");
    181182      ASSERT_EXCEPTION(retval, INTERNAL_ABORT,
    182                        "get_deltas_for_wrong_inertia returns false.");
     183                       "get_deltas_for_wrong_inertia returns false for degenerate Hessian before any further modification.");
    183184    }
    184185    else {
     
    264265        retval = get_deltas_for_wrong_inertia(delta_x, delta_s,
    265266                                              delta_c, delta_d);
    266         ASSERT_EXCEPTION(retval, INTERNAL_ABORT,
    267                          "get_deltas_for_wrong_inertia returns false.");
     267        if (!retval) {
     268          Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA,
     269                         "Can't get_deltas_for_wrong_inertia for delta_x_curr_ = %e and delta_c_curr_ = %e\n",
     270                         delta_x_curr_, delta_c_curr_);
     271          return false;
     272        }
    268273      }
    269274      else {
  • branches/dev/Algorithm/IpPDSystemSolver.hpp

    r501 r529  
    9999     *  system to best accuracy; for example, we don't want iterative
    100100     *  refinement during the computation of the second order
    101      *  correction.
     101     *  correction.  On the other hand, if improve_solution is true,
     102     *  the solution given in res should be improved (here beta has to
     103     *  be zero, and res is assume to be the solution for the system
     104     *  using rhs, without the factor alpha...).
    102105     */
    103106    virtual void Solve(Number alpha,
     
    105108                       const IteratesVector& rhs,
    106109                       IteratesVector& res,
    107                        bool allow_inexact=false) =0;
     110                       bool allow_inexact=false,
     111                       bool improve_solution=false) =0;
    108112
    109113  private:
  • branches/dev/Algorithm/IpProbingMuOracle.cpp

    r526 r529  
    7575    SmartPtr<IteratesVector> step = rhs->MakeNewIteratesVector(true);
    7676
    77     // Now solve the primal-dual system to get the step
     77    // Now solve the primal-dual system to get the affine step.  We
     78    // allow a somewhat inexact solution here
     79    bool allow_inexact = true;
    7880    pd_solver_->Solve(-1.0, 0.0,
    7981                      *rhs,
    80                       *step
    81                       //                      true           // don't need high accuracy
     82                      *step,
     83                      allow_inexact
    8284                     );
    8385
  • branches/dev/Algorithm/IpQualityFunctionMuOracle.cpp

    r527 r529  
    6262      "Maximum value of the centering parameter.",
    6363      0.0, true, 1e2);
     64    roptions->AddLowerBoundedNumberOption(
     65      "sigma_min",
     66      "Minimum value of the centering parameter.",
     67      0.0, false, 1e-6);
    6468    roptions->AddStringOption4(
    6569      "quality_function_norm_type",
     
    122126
    123127    options.GetNumericValue("sigma_max", sigma_max_, prefix);
     128    options.GetNumericValue("sigma_min", sigma_min_, prefix);
    124129
    125130    options.GetEnumValue("quality_function_norm_type", enum_int, prefix);
     
    136141                            quality_function_section_qf_tol_, prefix);
    137142
     143    initialized_ = false;
     144
    138145    return true;
    139146  }
     
    143150    DBG_START_METH("QualityFunctionMuOracle::CalculateMu",
    144151                   dbg_verbosity);
    145 
    146     count_qf_evals_ = 0;
    147152
    148153    ///////////////////////////////////////////////////////////////////////////
     
    188193    SmartPtr<IteratesVector> step_aff = IpData().curr()->MakeNewIteratesVector(true);
    189194
    190     // Now solve the primal-dual system to get the step
     195    // Now solve the primal-dual system to get the step.  We allow a
     196    // somewhat inexact solution, iterative refinement will be done
     197    // after mu is known
     198    bool allow_inexact = true;
    191199    pd_solver_->Solve(-1.0, 0.0,
    192200                      *rhs_aff,
    193201                      *step_aff,
    194                       false           // want accurate solution here
    195                       // because we can use it to
    196                       // compute the overall search
    197                       // direction
     202                      allow_inexact
    198203                     );
    199204
     
    228233
    229234    // Now solve the primal-dual system to get the step
     235    allow_inexact = true;
    230236    pd_solver_->Solve(1.0, 0.0,
    231237                      *rhs_cen,
    232238                      *step_cen,
    233                       false           // want accurate solution here
    234                       // because we can use it to
    235                       // compute the overall search
    236                       // direction
     239                      allow_inexact
    237240                     );
    238241
     
    240243
    241244    // Start the timing for the quality function search here
    242     IpData().TimingStats().QualityFunctionSearch.Start();
     245    IpData().TimingStats().QualityFunctionSearch().Start();
     246
     247    // Some initializations
     248    if (!initialized_) {
     249      n_dual_ = IpData().curr()->x()->Dim() + IpData().curr()->s()->Dim();
     250      n_pri_ = IpData().curr()->y_c()->Dim() + IpData().curr()->y_d()->Dim();
     251      n_comp_ = IpData().curr()->z_L()->Dim() + IpData().curr()->z_U()->Dim() +
     252                IpData().curr()->v_L()->Dim() + IpData().curr()->v_U()->Dim();
     253      initialized_ = true;
     254    }
     255
     256    count_qf_evals_ = 0;
     257
     258    // Compute some quantities used for the quality function evaluations
     259    // (This way we try to avoid retrieving numbers from cache...
     260
     261    curr_slack_x_L_ = IpCq().curr_slack_x_L();
     262    curr_slack_x_U_ = IpCq().curr_slack_x_U();
     263    curr_slack_s_L_ = IpCq().curr_slack_s_L();
     264    curr_slack_s_U_ = IpCq().curr_slack_s_U();
     265
     266    curr_z_L_ = IpData().curr()->z_L();
     267    curr_z_U_ = IpData().curr()->z_U();
     268    curr_v_L_ = IpData().curr()->v_L();
     269    curr_v_U_ = IpData().curr()->v_U();
     270
     271    IpData().TimingStats().Task5().Start();
     272    switch (quality_function_norm_) {
     273      case NM_NORM_1:
     274      curr_grad_lag_x_asum_ = IpCq().curr_grad_lag_x()->Asum();
     275      curr_grad_lag_s_asum_ = IpCq().curr_grad_lag_s()->Asum();
     276      curr_c_asum_ = IpCq().curr_c()->Asum();
     277      curr_d_minus_s_asum_ = IpCq().curr_d_minus_s()->Asum();
     278      break;
     279      case NM_NORM_2_SQUARED:
     280      case NM_NORM_2:
     281      curr_grad_lag_x_nrm2_ = IpCq().curr_grad_lag_x()->Nrm2();
     282      curr_grad_lag_s_nrm2_ = IpCq().curr_grad_lag_s()->Nrm2();
     283      curr_c_nrm2_ = IpCq().curr_c()->Nrm2();
     284      curr_d_minus_s_nrm2_ = IpCq().curr_d_minus_s()->Nrm2();
     285      break;
     286      case NM_NORM_MAX:
     287      curr_grad_lag_x_amax_ = IpCq().curr_grad_lag_x()->Amax();
     288      curr_grad_lag_s_amax_ = IpCq().curr_grad_lag_s()->Amax();
     289      curr_c_amax_ = IpCq().curr_c()->Amax();
     290      curr_d_minus_s_amax_ = IpCq().curr_d_minus_s()->Amax();
     291      break;
     292      default:
     293      DBG_ASSERT(false && "Unknown value for quality_function_norm_");
     294    }
     295    IpData().TimingStats().Task5().End();
     296
     297    // Some initializations
     298    if (!initialized_) {
     299      n_dual_ = IpData().curr()->x()->Dim() + IpData().curr()->s()->Dim();
     300      n_pri_ = IpData().curr()->y_c()->Dim() + IpData().curr()->y_d()->Dim();
     301      n_comp_ = IpData().curr()->z_L()->Dim() + IpData().curr()->z_U()->Dim() +
     302                IpData().curr()->v_L()->Dim() + IpData().curr()->v_U()->Dim();
     303      initialized_ = true;
     304    }
     305
     306    count_qf_evals_ = 0;
     307
     308    // Compute some quantities used for the quality function evaluations
     309    // (This way we try to avoid retrieving numbers from cache...
     310
     311    curr_slack_x_L_ = IpCq().curr_slack_x_L();
     312    curr_slack_x_U_ = IpCq().curr_slack_x_U();
     313    curr_slack_s_L_ = IpCq().curr_slack_s_L();
     314    curr_slack_s_U_ = IpCq().curr_slack_s_U();
     315
     316    curr_z_L_ = IpData().curr()->z_L();
     317    curr_z_U_ = IpData().curr()->z_U();
     318    curr_v_L_ = IpData().curr()->v_L();
     319    curr_v_U_ = IpData().curr()->v_U();
     320
     321    IpData().TimingStats().Task5().Start();
     322    switch (quality_function_norm_) {
     323      case NM_NORM_1:
     324      curr_grad_lag_x_asum_ = IpCq().curr_grad_lag_x()->Asum();
     325      curr_grad_lag_s_asum_ = IpCq().curr_grad_lag_s()->Asum();
     326      curr_c_asum_ = IpCq().curr_c()->Asum();
     327      curr_d_minus_s_asum_ = IpCq().curr_d_minus_s()->Asum();
     328      break;
     329      case NM_NORM_2_SQUARED:
     330      case NM_NORM_2:
     331      curr_grad_lag_x_nrm2_ = IpCq().curr_grad_lag_x()->Nrm2();
     332      curr_grad_lag_s_nrm2_ = IpCq().curr_grad_lag_s()->Nrm2();
     333      curr_c_nrm2_ = IpCq().curr_c()->Nrm2();
     334      curr_d_minus_s_nrm2_ = IpCq().curr_d_minus_s()->Nrm2();
     335      break;
     336      case NM_NORM_MAX:
     337      curr_grad_lag_x_amax_ = IpCq().curr_grad_lag_x()->Amax();
     338      curr_grad_lag_s_amax_ = IpCq().curr_grad_lag_s()->Amax();
     339      curr_c_amax_ = IpCq().curr_c()->Amax();
     340      curr_d_minus_s_amax_ = IpCq().curr_d_minus_s()->Amax();
     341      break;
     342      default:
     343      DBG_ASSERT(false && "Unknown value for quality_function_norm_");
     344    }
     345    IpData().TimingStats().Task5().End();
    243346
    244347    // We now compute the step for the slack variables.  This safes
     
    289392                                           *step_cen->v_U());
    290393
    291     Number sigma_1minus = 1.-quality_function_section_sigma_tol_;
     394    Number sigma_1minus = 1.-Max(1e-4, quality_function_section_sigma_tol_);
    292395    Number qf_1minus = CalculateQualityFunction(sigma_1minus,
    293396                       *step_aff_x_L,
     
    318421      Number sigma_up = Min(sigma_max_, mu_max/avrg_compl);
    319422      Number sigma_lo = 1.;
    320       // ToDo maybe we should use different tolerances for sigma>1
    321       sigma = PerformGoldenSection(sigma_up, -100., sigma_lo, qf_1,
    322                                    quality_function_section_sigma_tol_,
    323                                    quality_function_section_qf_tol_,
    324                                    *step_aff_x_L,
    325                                    *step_aff_x_U,
    326                                    *step_aff_s_L,
    327                                    *step_aff_s_U,
    328                                    *step_aff->y_c(),
    329                                    *step_aff->y_d(),
    330                                    *step_aff->z_L(),
    331                                    *step_aff->z_U(),
    332                                    *step_aff->v_L(),
    333                                    *step_aff->v_U(),
    334                                    *step_cen_x_L,
    335                                    *step_cen_x_U,
    336                                    *step_cen_s_L,
    337                                    *step_cen_s_U,
    338                                    *step_cen->y_c(),
    339                                    *step_cen->y_d(),
    340                                    *step_cen->z_L(),
    341                                    *step_cen->z_U(),
    342                                    *step_cen->v_L(),
    343                                    *step_cen->v_U());
     423      if (sigma_lo >= sigma_up) {
     424        sigma = sigma_up;
     425      }
     426      else {
     427        // ToDo maybe we should use different tolerances for sigma>1
     428        sigma = PerformGoldenSection(sigma_up, -100., sigma_lo, qf_1,
     429                                     quality_function_section_sigma_tol_,
     430                                     quality_function_section_qf_tol_,
     431                                     *step_aff_x_L,
     432                                     *step_aff_x_U,
     433                                     *step_aff_s_L,
     434                                     *step_aff_s_U,
     435                                     *step_aff->y_c(),
     436                                     *step_aff->y_d(),
     437                                     *step_aff->z_L(),
     438                                     *step_aff->z_U(),
     439                                     *step_aff->v_L(),
     440                                     *step_aff->v_U(),
     441                                     *step_cen_x_L,
     442                                     *step_cen_x_U,
     443                                     *step_cen_s_L,
     444                                     *step_cen_s_U,
     445                                     *step_cen->y_c(),
     446                                     *step_cen->y_d(),
     447                                     *step_cen->z_L(),
     448                                     *step_cen->z_U(),
     449                                     *step_cen->v_L(),
     450                                     *step_cen->v_U());
     451      }
    344452    }
    345453    else {
    346454      // Search for sigma less than 1
    347455
    348       Number sigma_lo = mu_min/avrg_compl;
    349       Number sigma_up = Max(sigma_lo, sigma_1minus);
    350       // Todo: What if lower bound is less than sigma_lo? Skip search?
    351       sigma = PerformGoldenSection(sigma_up, qf_1minus, sigma_lo, -100.,
    352                                    quality_function_section_sigma_tol_,
    353                                    quality_function_section_qf_tol_,
    354                                    *step_aff_x_L,
    355                                    *step_aff_x_U,
    356                                    *step_aff_s_L,
    357                                    *step_aff_s_U,
    358                                    *step_aff->y_c(),
    359                                    *step_aff->y_d(),
    360                                    *step_aff->z_L(),
    361                                    *step_aff->z_U(),
    362                                    *step_aff->v_L(),
    363                                    *step_aff->v_U(),
    364                                    *step_cen_x_L,
    365                                    *step_cen_x_U,
    366                                    *step_cen_s_L,
    367                                    *step_cen_s_U,
    368                                    *step_cen->y_c(),
    369                                    *step_cen->y_d(),
    370                                    *step_cen->z_L(),
    371                                    *step_cen->z_U(),
    372                                    *step_cen->v_L(),
    373                                    *step_cen->v_U());
     456      Number sigma_lo = Max(sigma_min_, mu_min/avrg_compl);
     457      Number sigma_up = Min(Max(sigma_lo, sigma_1minus), mu_max/avrg_compl);
     458      if (sigma_lo >= sigma_up) {
     459        // Skip the search, we are already at the minimum
     460        sigma = sigma_lo;
     461      }
     462      else {
     463        sigma = PerformGoldenSection(sigma_up, qf_1minus, sigma_lo, -100.,
     464                                     quality_function_section_sigma_tol_,
     465                                     quality_function_section_qf_tol_,
     466                                     *step_aff_x_L,
     467                                     *step_aff_x_U,
     468                                     *step_aff_s_L,
     469                                     *step_aff_s_U,
     470                                     *step_aff->y_c(),
     471                                     *step_aff->y_d(),
     472                                     *step_aff->z_L(),
     473                                     *step_aff->z_U(),
     474                                     *step_aff->v_L(),
     475                                     *step_aff->v_U(),
     476                                     *step_cen_x_L,
     477                                     *step_cen_x_U,
     478                                     *step_cen_s_L,
     479                                     *step_cen_s_U,
     480                                     *step_cen->y_c(),
     481                                     *step_cen->y_d(),
     482                                     *step_cen->z_L(),
     483                                     *step_cen->z_U(),
     484                                     *step_cen->v_L(),
     485                                     *step_cen->v_U());
     486      }
    374487    }
    375488
     
    412525
    413526    // End timing of quality function search
    414     IpData().TimingStats().QualityFunctionSearch.End();
     527    IpData().TimingStats().QualityFunctionSearch().End();
    415528
    416529    Jnlst().Printf(J_DETAILED, J_BARRIER_UPDATE,
     
    427540    step->AddTwoVectors(sigma, *step_cen, 1.0, *IpData().delta_aff(), 0.0);
    428541
     542    DBG_PRINT_VECTOR(2, "step", *step);
    429543    IpData().set_delta(step);
    430544    IpData().SetHaveDeltas(true);
     
    451565    tmp_v_L_ = NULL;
    452566    tmp_v_U_ = NULL;
     567
     568    curr_slack_x_L_ = NULL;
     569    curr_slack_x_U_ = NULL;
     570    curr_slack_s_L_ = NULL;
     571    curr_slack_s_U_ = NULL;
    453572
    454573    // DELETEME
     
    497616    count_qf_evals_++;
    498617
    499     Index n_dual = IpData().curr()->x()->Dim() + IpData().curr()->s()->Dim();
    500     Index n_pri = IpData().curr()->y_c()->Dim() + IpData().curr()->y_d()->Dim();
    501     Index n_comp = IpData().curr()->z_L()->Dim() + IpData().curr()->z_U()->Dim() +
    502                    IpData().curr()->v_L()->Dim() + IpData().curr()->v_U()->Dim();
    503 
    504     IpData().TimingStats().Task1.Start();
     618    IpData().TimingStats().Task1().Start();
    505619    tmp_step_x_L_->AddTwoVectors(1., step_aff_x_L, sigma, step_cen_x_L, 0.);
    506620    tmp_step_x_U_->AddTwoVectors(1., step_aff_x_U, sigma, step_cen_x_U, 0.);
     
    511625    tmp_step_v_L_->AddTwoVectors(1., step_aff_v_L, sigma, step_cen_v_L, 0.);
    512626    tmp_step_v_U_->AddTwoVectors(1., step_aff_v_U, sigma, step_cen_v_U, 0.);
    513     IpData().TimingStats().Task1.End();
     627    IpData().TimingStats().Task1().End();
    514628
    515629    // Compute the fraction-to-the-boundary step sizes
    516     IpData().TimingStats().Task2.Start();
     630    IpData().TimingStats().Task2().Start();
    517631    Number tau = IpData().curr_tau();
    518632    Number alpha_primal = IpCq().uncached_slack_frac_to_the_bound(tau,
     
    527641                        *tmp_step_v_L_,
    528642                        *tmp_step_v_U_);
    529     IpData().TimingStats().Task2.End();
     643    IpData().TimingStats().Task2().End();
    530644
    531645    Number xi = 0.; // centrality measure
    532646
    533     IpData().TimingStats().Task1.Start();
    534     tmp_slack_x_L_->AddTwoVectors(1., *IpCq().curr_slack_x_L(),
     647    IpData().TimingStats().Task1().Start();
     648    tmp_slack_x_L_->AddTwoVectors(1., *curr_slack_x_L_,
    535649                                  alpha_primal, *tmp_step_x_L_, 0.);
    536     tmp_slack_x_U_->AddTwoVectors(1., *IpCq().curr_slack_x_U(),
     650    tmp_slack_x_U_->AddTwoVectors(1., *curr_slack_x_U_,
    537651                                  alpha_primal, *tmp_step_x_U_, 0.);
    538     tmp_slack_s_L_->AddTwoVectors(1., *IpCq().curr_slack_s_L(),
     652    tmp_slack_s_L_->AddTwoVectors(1., *curr_slack_s_L_,
    539653                                  alpha_primal, *tmp_step_s_L_, 0.);
    540     tmp_slack_s_U_->AddTwoVectors(1., *IpCq().curr_slack_s_U(),
     654    tmp_slack_s_U_->AddTwoVectors(1., *curr_slack_s_U_,
    541655                                  alpha_primal, *tmp_step_s_U_, 0.);
    542656
    543     tmp_z_L_->AddTwoVectors(1., *IpData().curr()->z_L(),
     657    tmp_z_L_->AddTwoVectors(1., *curr_z_L_,
    544658                            alpha_dual, *tmp_step_z_L_, 0.);
    545     tmp_z_U_->AddTwoVectors(1., *IpData().curr()->z_U(),
     659    tmp_z_U_->AddTwoVectors(1., *curr_z_U_,
    546660                            alpha_dual, *tmp_step_z_U_, 0.);
    547     tmp_v_L_->AddTwoVectors(1., *IpData().curr()->v_L(),
     661    tmp_v_L_->AddTwoVectors(1., *curr_v_L_,
    548662                            alpha_dual, *tmp_step_v_L_, 0.);
    549     tmp_v_U_->AddTwoVectors(1., *IpData().curr()->v_U(),
     663    tmp_v_U_->AddTwoVectors(1., *curr_v_U_,
    550664                            alpha_dual, *tmp_step_v_U_, 0.);
    551     IpData().TimingStats().Task1.End();
    552 
    553     IpData().TimingStats().Task3.Start();
     665    IpData().TimingStats().Task1().End();
     666
     667    IpData().TimingStats().Task3().Start();
    554668    tmp_slack_x_L_->ElementWiseMultiply(*tmp_z_L_);
    555669    tmp_slack_x_U_->ElementWiseMultiply(*tmp_z_U_);
    556670    tmp_slack_s_L_->ElementWiseMultiply(*tmp_v_L_);
    557671    tmp_slack_s_U_->ElementWiseMultiply(*tmp_v_U_);
    558     IpData().TimingStats().Task3.End();
     672    IpData().TimingStats().Task3().End();
    559673
    560674    DBG_PRINT_VECTOR(2, "compl_x_L", *tmp_slack_x_L_);
     
    567681    Number compl_inf=-1.;
    568682
    569     IpData().TimingStats().Task5.Start();
     683    IpData().TimingStats().Task5().Start();
    570684    switch (quality_function_norm_) {
    571685      case NM_NORM_1:
    572       dual_inf = (1.-alpha_dual)*(IpCq().curr_grad_lag_x()->Asum() +
    573                                   IpCq().curr_grad_lag_s()->Asum());
    574 
    575       primal_inf = (1.-alpha_primal)*(IpCq().curr_c()->Asum() +
    576                                       IpCq().curr_d_minus_s()->Asum());
     686      dual_inf = (1.-alpha_dual)*(curr_grad_lag_x_asum_ +
     687                                  curr_grad_lag_s_asum_);
     688
     689      primal_inf = (1.-alpha_primal)*(curr_c_asum_ +
     690                                      curr_d_minus_s_asum_);
    577691
    578692      compl_inf = tmp_slack_x_L_->Asum() + tmp_slack_x_U_->Asum() +
    579693                  tmp_slack_s_L_->Asum() + tmp_slack_s_U_->Asum();
    580694
    581       dual_inf /= n_dual;
    582       if (n_pri>0) {
    583         primal_inf /= n_pri;
    584       }
    585       DBG_ASSERT(n_comp>0);
    586       compl_inf /= n_comp;
     695      dual_inf /= n_dual_;
     696      if (n_pri_>0) {
     697        primal_inf /= n_pri_;
     698      }
     699      DBG_ASSERT(n_comp_>0);
     700      compl_inf /= n_comp_;
    587701      break;
    588702      case NM_NORM_2_SQUARED:
    589703      dual_inf =
    590         pow(1.-alpha_dual, 2)*(pow(IpCq().curr_grad_lag_x()->Nrm2(), 2) +
    591                                pow(IpCq().curr_grad_lag_s()->Nrm2(), 2));
     704        pow(1.-alpha_dual, 2)*(pow(curr_grad_lag_x_nrm2_, 2) +
     705                               pow(curr_grad_lag_s_nrm2_, 2));
    592706      primal_inf =
    593         pow(1.-alpha_primal, 2)*(pow(IpCq().curr_c()->Nrm2(), 2) +
    594                                  pow(IpCq().curr_d_minus_s()->Nrm2(), 2));
     707        pow(1.-alpha_primal, 2)*(pow(curr_c_nrm2_, 2) +
     708                                 pow(curr_d_minus_s_nrm2_, 2));
    595709      compl_inf =
    596710        pow(tmp_slack_x_L_->Nrm2(), 2) + pow(tmp_slack_x_U_->Nrm2(), 2) +
    597711        pow(tmp_slack_s_L_->Nrm2(), 2) + pow(tmp_slack_s_U_->Nrm2(), 2);
    598712
    599       dual_inf /= n_dual;
    600       if (n_pri>0) {
    601         primal_inf /= n_pri;
    602       }
    603       DBG_ASSERT(n_comp>0);
    604       compl_inf /= n_comp;
     713      dual_inf /= n_dual_;
     714      if (n_pri_>0) {
     715        primal_inf /= n_pri_;
     716      }
     717      DBG_ASSERT(n_comp_>0);
     718      compl_inf /= n_comp_;
    605719      break;
    606720      case NM_NORM_MAX:
    607721      dual_inf =
    608         (1.-alpha_dual)*Max(IpCq().curr_grad_lag_x()->Amax(),
    609                             IpCq().curr_grad_lag_s()->Amax());
     722        (1.-alpha_dual)*Max(curr_grad_lag_x_amax_,
     723                            curr_grad_lag_s_amax_);
    610724      primal_inf =
    611         (1.-alpha_primal)*Max(IpCq().curr_c()->Amax(),
    612                               IpCq().curr_d_minus_s()->Amax());
     725        (1.-alpha_primal)*Max(curr_c_amax_,
     726                              curr_d_minus_s_amax_);
    613727      compl_inf =
    614728        Max(tmp_slack_x_L_->Amax(), tmp_slack_x_U_->Amax(),
     
    617731      case NM_NORM_2:
    618732      dual_inf =
    619         (1.-alpha_dual)*sqrt(pow(IpCq().curr_grad_lag_x()->Nrm2(), 2) +
    620                              pow(IpCq().curr_grad_lag_s()->Nrm2(), 2));
     733        (1.-alpha_dual)*sqrt(pow(curr_grad_lag_x_nrm2_, 2) +
     734                             pow(curr_grad_lag_s_nrm2_, 2));
    621735      primal_inf =
    622         (1.-alpha_primal)*sqrt(pow(IpCq().curr_c()->Nrm2(), 2) +
    623                                pow(IpCq().curr_d_minus_s()->Nrm2(), 2));
     736        (1.-alpha_primal)*sqrt(pow(curr_c_nrm2_, 2) +
     737                               pow(curr_d_minus_s_nrm2_, 2));
    624738      compl_inf =
    625739        sqrt(pow(tmp_slack_x_L_->Nrm2(), 2) + pow(tmp_slack_x_U_->Nrm2(), 2) +
    626740             pow(tmp_slack_s_L_->Nrm2(), 2) + pow(tmp_slack_s_U_->Nrm2(), 2));
    627741
    628       dual_inf /= sqrt((Number)n_dual);
    629       if (n_pri>0) {
    630         primal_inf /= sqrt((Number)n_pri);
    631       }
    632       DBG_ASSERT(n_comp>0);
    633       compl_inf /= sqrt((Number)n_comp);
     742      dual_inf /= sqrt((Number)n_dual_);
     743      if (n_pri_>0) {
     744        primal_inf /= sqrt((Number)n_pri_);
     745      }
     746      DBG_ASSERT(n_comp_>0);
     747      compl_inf /= sqrt((Number)n_comp_);
    634748      break;
    635749      default:
    636750      DBG_ASSERT(false && "Unknown value for quality_function_norm_");
    637751    }
    638     IpData().TimingStats().Task5.End();
     752    IpData().TimingStats().Task5().End();
    639753
    640754    Number quality_function = dual_inf + primal_inf + compl_inf;
    641755
    642756    if (quality_function_centrality_!=CEN_NONE) {
    643       IpData().TimingStats().Task4.Start();
     757      IpData().TimingStats().Task4().Start();
    644758      xi = IpCq().CalcCentralityMeasure(*tmp_slack_x_L_, *tmp_slack_x_U_,
    645759                                        *tmp_slack_s_L_, *tmp_slack_s_U_);
    646       IpData().TimingStats().Task4.End();
     760      IpData().TimingStats().Task4().End();
    647761    }
    648762    switch (quality_function_centrality_) {
  • branches/dev/Algorithm/IpQualityFunctionMuOracle.hpp

    r527 r529  
    189189    /** Upper bound on centering parameter sigma */
    190190    Number sigma_max_;
     191    /** Lower bound on centering parameter sigma */
     192    Number sigma_min_;
    191193    /** Norm to be used for the quality function. */
    192194    NormEnum quality_function_norm_;
     
    231233    //@}
    232234
     235    /* Counter for the qualify function evaluations */
    233236    Index count_qf_evals_;
     237
     238    /**@name Quantities used many times in CalculateQualityFunction,
     239     * which we store here instead of retrieving them from cache every
     240     * time.  I (AW) don't know if that really makes a difference, but
     241     * some of those things showed up in gprof. */
     242    //@{
     243    bool initialized_;
     244    Index n_dual_;
     245    Index n_pri_;
     246    Index n_comp_;
     247
     248    SmartPtr<const Vector> curr_slack_x_L_;
     249    SmartPtr<const Vector> curr_slack_x_U_;
     250    SmartPtr<const Vector> curr_slack_s_L_;
     251    SmartPtr<const Vector> curr_slack_s_U_;
     252
     253    SmartPtr<const Vector> curr_z_L_;
     254    SmartPtr<const Vector> curr_z_U_;
     255    SmartPtr<const Vector> curr_v_L_;
     256    SmartPtr<const Vector> curr_v_U_;
     257
     258    Number curr_grad_lag_x_asum_;
     259    Number curr_grad_lag_s_asum_;
     260    Number curr_c_asum_;
     261    Number curr_d_minus_s_asum_;
     262
     263    Number curr_grad_lag_x_nrm2_;
     264    Number curr_grad_lag_s_nrm2_;
     265    Number curr_c_nrm2_;
     266    Number curr_d_minus_s_nrm2_;
     267
     268    Number curr_grad_lag_x_amax_;
     269    Number curr_grad_lag_s_amax_;
     270    Number curr_c_amax_;
     271    Number curr_d_minus_s_amax_;
     272    //@}
    234273  };
    235274
  • branches/dev/Algorithm/IpTimingStatistics.cpp

    r526 r529  
    2121
    2222    jnlst.Printf(level, category,
    23                  "OverallAlgorithm....................: %10.2f\n",
    24                  OverallAlgorithm.TotalTime());
     23                 "OverallAlgorithm....................: %10.3f\n",
     24                 OverallAlgorithm_.TotalTime());
    2525    jnlst.Printf(level, category,
    26                  " PrintProblemStatistics.............: %10.2f\n",
    27                  PrintProblemStatistics.TotalTime());
     26                 " PrintProblemStatistics.............: %10.3f\n",
     27                 PrintProblemStatistics_.TotalTime());
    2828    jnlst.Printf(level, category,
    29                  " InitializeIterates.................: %10.2f\n",
    30                  InitializeIterates.TotalTime());
     29                 " InitializeIterates.................: %10.3f\n",
     30                 InitializeIterates_.TotalTime());
    3131    jnlst.Printf(level, category,
    32                  " ActualizeHessian...................: %10.2f\n",
    33                  ActualizeHessian.TotalTime());
     32                 " ActualizeHessian...................: %10.3f\n",
     33                 ActualizeHessian_.TotalTime());
    3434    jnlst.Printf(level, category,
    35                  " OutputIteration....................: %10.2f\n",
    36                  OutputIteration.TotalTime());
     35                 " OutputIteration....................: %10.3f\n",
     36                 OutputIteration_.TotalTime());
    3737    jnlst.Printf(level, category,
    38                  " UpdateBarrierParameter.............: %10.2f\n",
    39                  UpdateBarrierParameter.TotalTime());
     38                 " UpdateBarrierParameter.............: %10.3f\n",
     39                 UpdateBarrierParameter_.TotalTime());
    4040    jnlst.Printf(level, category,
    41                  " ComputeSearchDirection.............: %10.2f\n",
    42                  ComputeSearchDirection.TotalTime());
     41                 " ComputeSearchDirection.............: %10.3f\n",
     42                 ComputeSearchDirection_.TotalTime());
    4343    jnlst.Printf(level, category,
    44                  " ComputeAcceptableTrialPoint........: %10.2f\n",
    45                  ComputeAcceptableTrialPoint.TotalTime());
     44                 " ComputeAcceptableTrialPoint........: %10.3f\n",
     45                 ComputeAcceptableTrialPoint_.TotalTime());
    4646    jnlst.Printf(level, category,
    47                  " AcceptTrialPoint...................: %10.2f\n",
    48                  AcceptTrialPoint.TotalTime());
     47                 " AcceptTrialPoint...................: %10.3f\n",
     48                 AcceptTrialPoint_.TotalTime());
    4949    jnlst.Printf(level, category,
    50                  " CheckConvergence...................: %10.2f\n",
    51                  CheckConvergence.TotalTime());
     50                 " CheckConvergence...................: %10.3f\n",
     51                 CheckConvergence_.TotalTime());
    5252
    5353    jnlst.Printf(level, category,
    54                  "PDSystemSolverTotal.................: %10.2f\n",
    55                  PDSystemSolverTotal.TotalTime());
     54                 "PDSystemSolverTotal.................: %10.3f\n",
     55                 PDSystemSolverTotal_.TotalTime());
    5656    jnlst.Printf(level, category,
    57                  " PDSystemSolverSolveOnce............: %10.2f\n",
    58                  PDSystemSolverSolveOnce.TotalTime());
     57                 " PDSystemSolverSolveOnce............: %10.3f\n",
     58                 PDSystemSolverSolveOnce_.TotalTime());
    5959    jnlst.Printf(level, category,
    60                  " LinearSystemScaling................: %10.2f\n",
    61                  LinearSystemScaling.TotalTime());
     60                 " ComputeResiduals...................: %10.3f\n",
     61                 ComputeResiduals_.TotalTime());
    6262    jnlst.Printf(level, category,
    63                  " LinearSystemSymbolicFactorization..: %10.2f\n",
    64                  LinearSystemSymbolicFactorization.TotalTime());
     63                 " LinearSystemScaling................: %10.3f\n",
     64                 LinearSystemScaling_.TotalTime());
    6565    jnlst.Printf(level, category,
    66                  " LinearSystemFactorization..........: %10.2f\n",
    67                  LinearSystemFactorization.TotalTime());
     66                 " LinearSystemSymbolicFactorization..: %10.3f\n",
     67                 LinearSystemSymbolicFactorization_.TotalTime());
    6868    jnlst.Printf(level, category,
    69                  " LinearSystemBackSolve..............: %10.2f\n",
    70                  LinearSystemBackSolve.TotalTime());
     69                 " LinearSystemFactorization..........: %10.3f\n",
     70                 LinearSystemFactorization_.TotalTime());
    7171    jnlst.Printf(level, category,
    72                  "QualityFunctionSearch...............: %10.2f\n",
    73                  QualityFunctionSearch.TotalTime());
     72                 " LinearSystemBackSolve..............: %10.3f\n",
     73                 LinearSystemBackSolve_.TotalTime());
    7474    jnlst.Printf(level, category,
    75                  "TryCorrector........................: %10.2f\n",
    76                  TryCorrector.TotalTime());
     75                 "QualityFunctionSearch...............: %10.3f\n",
     76                 QualityFunctionSearch_.TotalTime());
    7777    jnlst.Printf(level, category,
    78                  "Task1...............................: %10.2f\n",
    79                  Task1.TotalTime());
     78                 "TryCorrector........................: %10.3f\n",
     79                 TryCorrector_.TotalTime());
    8080    jnlst.Printf(level, category,
    81                  "Task2...............................: %10.2f\n",
    82                  Task2.TotalTime());
     81                 "Task1...............................: %10.3f\n",
     82                 Task1_.TotalTime());
    8383    jnlst.Printf(level, category,
    84                  "Task3...............................: %10.2f\n",
    85                  Task3.TotalTime());
     84                 "Task2...............................: %10.3f\n",
     85                 Task2_.TotalTime());
    8686    jnlst.Printf(level, category,
    87                  "Task4...............................: %10.2f\n",
    88                  Task4.TotalTime());
     87                 "Task3...............................: %10.3f\n",
     88                 Task3_.TotalTime());
    8989    jnlst.Printf(level, category,
    90                  "Task5...............................: %10.2f\n",
    91                  Task5.TotalTime());
     90                 "Task4...............................: %10.3f\n",
     91                 Task4_.TotalTime());
     92    jnlst.Printf(level, category,
     93                 "Task5...............................: %10.3f\n",
     94                 Task5_.TotalTime());
    9295  }
    9396} // namespace Ipopt
  • branches/dev/Algorithm/IpTimingStatistics.hpp

    r528 r529  
    172172                                  EJournalCategory category) const;
    173173
    174     /**@name All timed tasks. */
    175     //@{
    176     TimedTask OverallAlgorithm;
    177     TimedTask PrintProblemStatistics;
    178     TimedTask InitializeIterates;
    179     TimedTask ActualizeHessian;
    180     TimedTask OutputIteration;
    181     TimedTask UpdateBarrierParameter;
    182     TimedTask ComputeSearchDirection;
    183     TimedTask ComputeAcceptableTrialPoint;
    184     TimedTask AcceptTrialPoint;
    185     TimedTask CheckConvergence;
    186 
    187     TimedTask PDSystemSolverTotal;
    188     TimedTask PDSystemSolverSolveOnce;
    189     TimedTask LinearSystemScaling;
    190     TimedTask LinearSystemSymbolicFactorization;
    191     TimedTask LinearSystemFactorization;
    192     TimedTask LinearSystemBackSolve;
    193     TimedTask QualityFunctionSearch;
    194     TimedTask TryCorrector;
    195 
    196     TimedTask Task1;
    197     TimedTask Task2;
    198     TimedTask Task3;
    199     TimedTask Task4;
    200     TimedTask Task5;
    201     TimedTask Task6;
     174    /**@name Accessor methods to all timed tasks. */
     175    //@{
     176    TimedTask& OverallAlgorithm()
     177    {
     178      return OverallAlgorithm_;
     179    }
     180    TimedTask& PrintProblemStatistics()
     181    {
     182      return PrintProblemStatistics_;
     183    }
     184    TimedTask& InitializeIterates()
     185    {
     186      return InitializeIterates_;
     187    }
     188    TimedTask& ActualizeHessian()
     189    {
     190      return ActualizeHessian_;
     191    }
     192    TimedTask& OutputIteration()
     193    {
     194      return OutputIteration_;
     195    }
     196    TimedTask& UpdateBarrierParameter()
     197    {
     198      return UpdateBarrierParameter_;
     199    }
     200    TimedTask& ComputeSearchDirection()
     201    {
     202      return ComputeSearchDirection_;
     203    }
     204    TimedTask& ComputeAcceptableTrialPoint()
     205    {
     206      return ComputeAcceptableTrialPoint_;
     207    }
     208    TimedTask& AcceptTrialPoint()
     209    {
     210      return AcceptTrialPoint_;
     211    }
     212    TimedTask& CheckConvergence()
     213    {
     214      return CheckConvergence_;
     215    }
     216
     217    TimedTask& PDSystemSolverTotal()
     218    {
     219      return PDSystemSolverTotal_;
     220    }
     221    TimedTask& PDSystemSolverSolveOnce()
     222    {
     223      return PDSystemSolverSolveOnce_;
     224    }
     225    TimedTask& ComputeResiduals()
     226    {
     227      return ComputeResiduals_;
     228    }
     229    TimedTask& LinearSystemScaling()
     230    {
     231      return LinearSystemScaling_;
     232    }
     233    TimedTask& LinearSystemSymbolicFactorization()
     234    {
     235      return LinearSystemSymbolicFactorization_;
     236    }
     237    TimedTask& LinearSystemFactorization()
     238    {
     239      return LinearSystemFactorization_;
     240    }
     241    TimedTask& LinearSystemBackSolve()
     242    {
     243      return LinearSystemBackSolve_;
     244    }
     245    TimedTask& QualityFunctionSearch()
     246    {
     247      return QualityFunctionSearch_;
     248    }
     249    TimedTask& TryCorrector()
     250    {
     251      return TryCorrector_;
     252    }
     253
     254    TimedTask& Task1()
     255    {
     256      return Task1_;
     257    }
     258    TimedTask& Task2()
     259    {
     260      return Task2_;
     261    }
     262    TimedTask& Task3()
     263    {
     264      return Task3_;
     265    }
     266    TimedTask& Task4()
     267    {
     268      return Task4_;
     269    }
     270    TimedTask& Task5()
     271    {
     272      return Task5_;
     273    }
     274    TimedTask& Task6()
     275    {
     276      return Task6_;
     277    }
    202278    //@}
    203279
     
    218294    //@}
    219295
     296    /**@name All timed tasks. */
     297    //@{
     298    TimedTask OverallAlgorithm_;
     299    TimedTask PrintProblemStatistics_;
     300    TimedTask InitializeIterates_;
     301    TimedTask ActualizeHessian_;
     302    TimedTask OutputIteration_;
     303    TimedTask UpdateBarrierParameter_;
     304    TimedTask ComputeSearchDirection_;
     305    TimedTask ComputeAcceptableTrialPoint_;
     306    TimedTask AcceptTrialPoint_;
     307    TimedTask CheckConvergence_;
     308
     309    TimedTask PDSystemSolverTotal_;
     310    TimedTask PDSystemSolverSolveOnce_;
     311    TimedTask ComputeResiduals_;
     312    TimedTask LinearSystemScaling_;
     313    TimedTask LinearSystemSymbolicFactorization_;
     314    TimedTask LinearSystemFactorization_;
     315    TimedTask LinearSystemBackSolve_;
     316    TimedTask QualityFunctionSearch_;
     317    TimedTask TryCorrector_;
     318
     319    TimedTask Task1_;
     320    TimedTask Task2_;
     321    TimedTask Task3_;
     322    TimedTask Task4_;
     323    TimedTask Task5_;
     324    TimedTask Task6_;
     325    //@}
    220326  };
    221327
  • branches/dev/Algorithm/LinearSolvers/IpMa27TSolverInterface.cpp

    r524 r529  
    251251    DBG_START_METH("Ma27TSolverInterface::SymbolicFactorization",dbg_verbosity);
    252252
    253     IpData().TimingStats().LinearSystemSymbolicFactorization.Start();
     253    IpData().TimingStats().LinearSystemSymbolicFactorization().Start();
    254254
    255255    // Get memory for the IW workspace
     
    308308    a_ = new double[la_];
    309309
    310     IpData().TimingStats().LinearSystemSymbolicFactorization.End();
     310    IpData().TimingStats().LinearSystemSymbolicFactorization().End();
    311311
    312312    return SYMSOLVER_SUCCESS;
     
    321321    DBG_START_METH("Ma27TSolverInterface::Factorization",dbg_verbosity);
    322322    // Check if la should be increased
    323     IpData().TimingStats().LinearSystemFactorization.Start();
     323    IpData().TimingStats().LinearSystemFactorization().Start();
    324324    if (la_increase_) {
    325325      double* a_old = a_;
     
    397397                     "MA27BD returned iflag=%d.\n Increase liw from %d to %d and la from %d to %d and factorize again.\n",
    398398                     iflag, liw_old, liw_, la_old, la_);
    399       IpData().TimingStats().LinearSystemFactorization.End();
     399      IpData().TimingStats().LinearSystemFactorization().End();
    400400      return SYMSOLVER_CALL_AGAIN;
    401401    }
     
    403403    // Check if the system is singular, and if some other error occurred
    404404    if (iflag==-5 || iflag==3) {
    405       IpData().TimingStats().LinearSystemFactorization.End();
     405      IpData().TimingStats().LinearSystemFactorization().End();
    406406      return SYMSOLVER_SINGULAR;
    407407    }
    408408    else if (iflag != 0) {
    409409      // There is some error
    410       IpData().TimingStats().LinearSystemFactorization.End();
     410      IpData().TimingStats().LinearSystemFactorization().End();
    411411      return SYMSOLVER_FATAL_ERROR;
    412412    }
     
    433433                     "In Ma27TSolverInterface::Factorization: negevals_ = %d, but numberOfNegEVals = %d\n",
    434434                     negevals_, numberOfNegEVals);
    435       IpData().TimingStats().LinearSystemFactorization.End();
     435      IpData().TimingStats().LinearSystemFactorization().End();
    436436      return SYMSOLVER_WRONG_INERTIA;
    437437    }
    438438
    439     IpData().TimingStats().LinearSystemFactorization.End();
     439    IpData().TimingStats().LinearSystemFactorization().End();
    440440    return SYMSOLVER_SUCCESS;
    441441  }
     
    445445  {
    446446    DBG_START_METH("Ma27TSolverInterface::Backsolve",dbg_verbosity);
    447     IpData().TimingStats().LinearSystemBackSolve.Start();
     447    IpData().TimingStats().LinearSystemBackSolve().Start();
    448448
    449449    ipfint N=dim_;
     
    472472    delete [] IW1;
    473473
    474     IpData().TimingStats().LinearSystemBackSolve.End();
     474    IpData().TimingStats().LinearSystemBackSolve().End();
    475475    return SYMSOLVER_SUCCESS;
    476476  }
  • branches/dev/Algorithm/LinearSolvers/IpTSymLinearSolver.cpp

    r523 r529  
    104104    bool retval = true;
    105105    if (IsValid(scaling_method_)) {
    106       IpData().TimingStats().LinearSystemScaling.Start();
     106      IpData().TimingStats().LinearSystemScaling().Start();
    107107      retval = scaling_method_->Initialize(Jnlst(), IpNLP(), IpData(), IpCq(),
    108108                                           options, prefix);
    109       IpData().TimingStats().LinearSystemScaling.End();
     109      IpData().TimingStats().LinearSystemScaling().End();
    110110    }
    111111    return retval;
     
    154154                                          &rhs_vals[irhs*(dim_)]);
    155155      if (IsValid(scaling_method_)) {
    156         IpData().TimingStats().LinearSystemScaling.Start();
     156        IpData().TimingStats().LinearSystemScaling().Start();
    157157        for (Index i=0; i<dim_; i++) {
    158158          rhs_vals[irhs*(dim_)+i] *= scaling_factors_[i];
    159159        }
    160         IpData().TimingStats().LinearSystemScaling.End();
     160        IpData().TimingStats().LinearSystemScaling().End();
    161161      }
    162162    }
     
    198198      for (Index irhs=0; irhs<nrhs; irhs++) {
    199199        if (IsValid(scaling_method_)) {
    200           IpData().TimingStats().LinearSystemScaling.Start();
     200          IpData().TimingStats().LinearSystemScaling().Start();
    201201          for (Index i=0; i<dim_; i++) {
    202202            rhs_vals[irhs*(dim_)+i] *= scaling_factors_[i];
    203203          }
    204           IpData().TimingStats().LinearSystemScaling.End();
     204          IpData().TimingStats().LinearSystemScaling().End();
    205205        }
    206206        TripletHelper::PutValuesInVector(dim_, &rhs_vals[irhs*(dim_)],
     
    266266      delete [] scaling_factors_;
    267267      if (IsValid(scaling_method_)) {
    268         IpData().TimingStats().LinearSystemScaling.Start();
     268        IpData().TimingStats().LinearSystemScaling().Start();
    269269        scaling_factors_ = new double[dim_];
    270         IpData().TimingStats().LinearSystemScaling.End();
     270        IpData().TimingStats().LinearSystemScaling().End();
    271271      }
    272272
     
    342342
    343343    if (IsValid(scaling_method_)) {
    344       IpData().TimingStats().LinearSystemScaling.Start();
     344      IpData().TimingStats().LinearSystemScaling().Start();
    345345      DBG_ASSERT(scaling_factors_);
    346346      if (new_matrix) {
     
    371371        }
    372372      }
    373       IpData().TimingStats().LinearSystemScaling.End();
     373      IpData().TimingStats().LinearSystemScaling().End();
    374374    }
    375375
  • branches/dev/Apps/CUTErInterface/CUTErInterface.f

    r527 r529  
    6464C
    6565      logical equatn(CUTE_MMAX), linear(CUTE_MMAX)
    66       integer cnr_input
    67       logical efirst, lfirst, nvfrst
     66      integer i, cnr_input
     67      logical efirst, lfirst, nvfrst, ex
     68      double precision init_val
    6869C
    6970C     Initialize the CUTEr interface and get the initial point
     
    8081     2     efirst, lfirst, nvfrst)
    8182      close(cnr_input)
     83C
     84C     See if we want to set a different initial point
     85C
     86      inquire(file='INITPOINT.VAL', exist=ex)
     87      if (ex) then
     88         open(70, file='INITPOINT.VAL', status='old')
     89         read(70,'(d25.16)') init_val
     90         do i = 1, N
     91            X(i) = init_val
     92         enddo
     93         close(70)
     94      endif
    8295C
    8396C     Obtain the number of nonzeros in Jacobian and Hessian
  • branches/dev/LinAlg/IpVector.hpp

    r501 r529  
    1717
    1818#include <vector>
     19
     20#ifdef HAVE_CMATH
     21# include <cmath>
     22#else
     23# ifdef HAVE_MATH_H
     24#  include <math.h>
     25# else
     26#  error "don't have header file for math"
     27# endif
     28#endif
    1929
    2030namespace Ipopt
     
    326336    mutable Number cached_sumlogs_;
    327337
    328     //     /** Cache for 2-norm */
    329     //     mutable CachedResults<Number> nrm2_cache_;
    330     //     /** Cache for Asum */
    331     //     mutable CachedResults<Number> asum_cache_;
    332     //     /** Cache for Amax */
    333     //     mutable CachedResults<Number> amax_cache_;
    334     //     /** Cache for Max */
    335     //     mutable CachedResults<Number> max_cache_;
    336     //     /** Cache for Min */
    337     //     mutable CachedResults<Number> min_cache_;
    338     //     /** Cache for Sum */
    339     //     mutable CachedResults<Number> sum_cache_;
    340     //     /** Cache for SumLogs */
    341     //     mutable CachedResults<Number> sumlogs_cache_;
    342     /** Cache for FracToBound */
    343     mutable CachedResults<Number> frac_to_bound_cache_;
     338    //     AW: I removed this cache since it gets in the way for the
     339    //         quality function search
     340    //     /** Cache for FracToBound */
     341    //     mutable CachedResults<Number> frac_to_bound_cache_;
    344342    //@}
    345343
     
    416414      min_cache_tag_(0),
    417415      sum_cache_tag_(0),
    418       sumlogs_cache_tag_(0),
    419       frac_to_bound_cache_(4)
     416      sumlogs_cache_tag_(0)
    420417  {
    421418    DBG_ASSERT(IsValid(owner_space_));
     
    431428  Vector* Vector::MakeNewCopy() const
    432429  {
     430    // ToDo: We can probably copy also the cached values for Norms etc here
    433431    Vector* copy = MakeNew();
    434432    copy->Copy(*this);
     
    441439    CopyImpl(x);
    442440    ObjectChanged();
     441    // Also copy any cached scalar values from the original vector
     442    // ToDo: Check if that is too much overhead
     443    TaggedObject::Tag x_tag = x.GetTag();
     444    if (x_tag == x.nrm2_cache_tag_) {
     445      nrm2_cache_tag_ = GetTag();
     446      cached_nrm2_ = x.cached_nrm2_;
     447    }
     448    if (x_tag == x.asum_cache_tag_) {
     449      asum_cache_tag_ = GetTag();
     450      cached_asum_ = x.cached_asum_;
     451    }
     452    if (x_tag == x.amax_cache_tag_) {
     453      amax_cache_tag_ = GetTag();
     454      cached_amax_ = x.cached_amax_;
     455    }
     456    if (x_tag == x.max_cache_tag_) {
     457      max_cache_tag_ = GetTag();
     458      cached_max_ = x.cached_max_;
     459    }
     460    if (x_tag == x.min_cache_tag_) {
     461      min_cache_tag_ = GetTag();
     462      cached_min_ = x.cached_min_;
     463    }
     464    if (x_tag == x.sum_cache_tag_) {
     465      sum_cache_tag_ = GetTag();
     466      cached_sum_ = x.cached_sum_;
     467    }
     468    if (x_tag == x.sumlogs_cache_tag_) {
     469      sumlogs_cache_tag_ = GetTag();
     470      cached_sumlogs_ = x.cached_sumlogs_;
     471    }
    443472  }
    444473
     
    447476  {
    448477    if (alpha!=1.) {
     478      TaggedObject::Tag old_tag = GetTag();
    449479      ScalImpl(alpha);
    450480      ObjectChanged();
     481      if (old_tag == nrm2_cache_tag_) {
     482        nrm2_cache_tag_ = GetTag();
     483        cached_nrm2_ *= fabs(alpha);
     484      }
     485      if (old_tag == asum_cache_tag_) {
     486        asum_cache_tag_ = GetTag();
     487        cached_asum_ *= fabs(alpha);
     488      }
     489      if (old_tag == amax_cache_tag_) {
     490        amax_cache_tag_ = GetTag();
     491        cached_amax_ *= fabs(alpha);
     492      }
     493      if (old_tag == max_cache_tag_) {
     494        if (alpha>=0.) {
     495          max_cache_tag_ = GetTag();
     496          cached_max_ *= alpha;
     497        }
     498        else if (alpha<0.) {
     499          min_cache_tag_ = GetTag();
     500          cached_min_ = cached_max_ * alpha;
     501        }
     502      }
     503      if (old_tag == min_cache_tag_) {
     504        if (alpha>=0.) {
     505          min_cache_tag_ = GetTag();
     506          cached_min_ *= alpha;
     507        }
     508        else if (alpha<0.) {
     509          max_cache_tag_ = GetTag();
     510          cached_max_ = cached_min_ * alpha;
     511        }
     512      }
     513      if (old_tag == sum_cache_tag_) {
     514        sum_cache_tag_ = GetTag();
     515        cached_sum_ *= alpha;
     516      }
     517      if (old_tag == sumlogs_cache_tag_) {
     518        sumlogs_cache_tag_ = GetTag();
     519        cached_sumlogs_ += ((Number)Dim())*log(alpha);
     520      }
    451521    }
    452522  }
     
    530600  void Vector::Set(Number alpha)
    531601  {
     602    // Could initialize caches here
    532603    SetImpl(alpha);
    533604    ObjectChanged();
     
    558629  void Vector::ElementWiseMax(const Vector& x)
    559630  {
     631    // Could initialize some caches here
    560632    ElementWiseMaxImpl(x);
    561633    ObjectChanged();
     
    565637  void Vector::ElementWiseMin(const Vector& x)
    566638  {
     639    // Could initialize some caches here
    567640    ElementWiseMinImpl(x);
    568641    ObjectChanged();
     
    572645  void Vector::ElementWiseAbs()
    573646  {
     647    // Could initialize some caches here
    574648    ElementWiseAbsImpl();
    575649    ObjectChanged();
     
    586660  void Vector::AddScalar(Number scalar)
    587661  {
     662    // Could initialize some caches here
    588663    AddScalarImpl(scalar);
    589664    ObjectChanged();
     
    627702  Number Vector::FracToBound(const Vector& delta, Number tau) const
    628703  {
     704    /* AW: I avoid the caching here, since it leads to overhead in the
     705       quality function search.  Caches for this are in
     706       CalculatedQuantities.
    629707    Number retValue;
    630708    std::vector<const TaggedObject*> tdeps(1);
     
    637715    }
    638716    return retValue;
     717    */
     718    return FracToBoundImpl(delta, tau);
    639719  }
    640720
Note: See TracChangeset for help on using the changeset viewer.