C++ Jnlst函数代码示例

  bool PDFullSpaceSolver::InitializeImpl(const OptionsList& options,
                                         const std::string& prefix)
  {
    // Check for the algorithm options
    options.GetIntegerValue("min_refinement_steps", min_refinement_steps_, prefix);
    options.GetIntegerValue("max_refinement_steps", max_refinement_steps_, prefix);
    ASSERT_EXCEPTION(max_refinement_steps_ >= min_refinement_steps_, OPTION_INVALID,
                     "Option \"max_refinement_steps\": This value must be larger than or equal to min_refinement_steps (default 1)");

    options.GetNumericValue("residual_ratio_max", residual_ratio_max_, prefix);
    options.GetNumericValue("residual_ratio_singular", residual_ratio_singular_, prefix);
    ASSERT_EXCEPTION(residual_ratio_singular_ >= residual_ratio_max_, OPTION_INVALID,
                     "Option \"residual_ratio_singular\": This value must be not smaller than residual_ratio_max.");
    options.GetNumericValue("residual_improvement_factor", residual_improvement_factor_, prefix);
    options.GetNumericValue("neg_curv_test_tol", neg_curv_test_tol_, prefix);
    options.GetBoolValue("neg_curv_test_reg", neg_curv_test_reg_, prefix);

    // Reset internal flags and data
    augsys_improved_ = false;

    if (!augSysSolver_->Initialize(Jnlst(), IpNLP(), IpData(), IpCq(),
                                   options, prefix)) {
      return false;
    }

    return perturbHandler_->Initialize(Jnlst(), IpNLP(), IpData(), IpCq(),
                                       options, prefix);
  }

  bool WsmpSolverInterface::IncreaseQuality()
  {
    DBG_START_METH("WsmpSolverInterface::IncreaseQuality",dbg_verbosity);

    if (factorizations_since_recomputed_ordering_ == -1 ||
        factorizations_since_recomputed_ordering_ > 2) {
      DPARM_[14] = 1.0;
      pivtol_changed_ = true;
      IpData().Append_info_string("RO ");
      factorizations_since_recomputed_ordering_ = 0;
      Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA,
                     "Triggering WSMP's recomputation of the ordering for next factorization.\n");
      return true;
    }
    if (wsmp_pivtol_ == wsmp_pivtolmax_) {
      return false;
    }
    pivtol_changed_ = true;

    Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA,
                   "Increasing pivot tolerance for WSMP from %7.2e ",
                   wsmp_pivtol_);
    wsmp_pivtol_ = Min(wsmp_pivtolmax_, pow(wsmp_pivtol_,0.75));
    Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA,
                   "to %7.2e.\n",
                   wsmp_pivtol_);
    return true;
  }

 bool CGPenaltyLSAcceptor::ArmijoHolds(Number alpha_primal_test)
 {
   DBG_START_METH("CGPenaltyLSAcceptor::ArmijoHolds",
                  dbg_verbosity);
   bool accept = false;
   Number trial_penalty_function = CGPenCq().trial_penalty_function();
   DBG_ASSERT(IsFiniteNumber(trial_penalty_function));
   Jnlst().Printf(J_DETAILED, J_LINE_SEARCH,
                  "Checking acceptability for trial step size alpha_primal_test=%13.6e:\n", alpha_primal_test);
   Jnlst().Printf(J_DETAILED, J_LINE_SEARCH,
                  " New values of penalty function     = %23.16e  (reference %23.16e):\n", trial_penalty_function, reference_penalty_function_);
   if (Jnlst().ProduceOutput(J_DETAILED, J_LINE_SEARCH)) {
     Jnlst().Printf(J_DETAILED, J_LINE_SEARCH,
                    "curr_barr  = %23.16e curr_inf  = %23.16e\n",
                    IpCq().curr_barrier_obj(),
                    IpCq().curr_constraint_violation());
     Jnlst().Printf(J_DETAILED, J_LINE_SEARCH,
                    "trial_barr = %23.16e trial_inf = %23.16e\n",
                    IpCq().trial_barrier_obj(),
                    IpCq().trial_constraint_violation());
   }
   // Now check the Armijo condition
   accept = Compare_le(trial_penalty_function-reference_penalty_function_,
                       eta_penalty_*alpha_primal_test*reference_direct_deriv_penalty_function_,
                       reference_penalty_function_);
   return accept;
 }

  char CGPenaltyLSAcceptor::UpdateForNextIteration(Number alpha_primal_test)
  {
    char info_alpha_primal_char='n';
    if (pen_curr_mu_ > IpData().curr_mu()) {
      pen_curr_mu_ = IpData().curr_mu();
      best_KKT_error_ = -1.;
    }
    // See if the current iterate has the least KKT errors so far
    // If so, store the current iterate
    if (CurrentIsBest()) {
      StoreBestPoint();
    }
    // update piecewise penalty parameters
    PiecewisePenalty_.Print( Jnlst() );
    if (!accepted_by_Armijo_) {
      PiecewisePenalty_.UpdateEntry(IpCq().trial_barrier_obj(),
                                    IpCq().trial_constraint_violation());
    }
    PiecewisePenalty_.Print( Jnlst() );

    // update regular penalty parameter
    if (CGPenData().CurrPenaltyPert() != 0) {
      info_alpha_primal_char = UpdatePenaltyParameter();
    }
    return info_alpha_primal_char;
  }

  /* The functions SchurSolve do IFT step, if S_==NULL, and DenseGenSchurDriver otherwise. */
  bool DenseGenSchurDriver::SchurSolve(SmartPtr<IteratesVector> lhs, // new left hand side will be stored here
				       SmartPtr<const IteratesVector> rhs, // rhs r_s
				       SmartPtr<Vector> delta_u,  // should be (u_p - u_0) WATCH OUT FOR THE SIGN! I like it this way, so that u_0+delta_u = u_p, but victor always used it the other way round, so be careful. At the end, delta_nu is saved in here.
				       SmartPtr<IteratesVector> sol) // the vector K^(-1)*r_s which usually should have been computed before.

  {
    DBG_START_METH("DenseGenSchurDriver::SchurSolve", dbg_verbosity);
    DBG_ASSERT(IsValid(S_));
    bool retval;

    // set up rhs of equation (3.48a)
    SmartPtr<Vector> delta_rhs = delta_u->MakeNew();
    data_B()->Multiply(*sol, *delta_rhs);
    delta_rhs->Print(Jnlst(),J_VECTOR,J_USER1,"delta_rhs");
    delta_rhs->Scal(-1.0);
    delta_rhs->Axpy(1.0, *delta_u);
    delta_rhs->Print(Jnlst(),J_VECTOR,J_USER1,"rhs 3.48a");

    // solve equation (3.48a) for delta_nu
    SmartPtr<DenseVector> delta_nu = dynamic_cast<DenseVector*>(GetRawPtr(delta_rhs))->MakeNewDenseVector();
    delta_nu->Copy(*delta_rhs);
    S_->LUSolveVector(*delta_nu); // why is LUSolveVector not bool??
    delta_nu->Print(Jnlst(),J_VECTOR,J_USER1,"delta_nu");

    // solve equation (3.48b) for lhs (=delta_s)
    SmartPtr<IteratesVector> new_rhs = lhs->MakeNewIteratesVector();
    data_A()->TransMultiply(*delta_nu, *new_rhs);
    new_rhs->Axpy(-1.0, *rhs);
    new_rhs->Scal(-1.0);
    new_rhs->Print(Jnlst(),J_VECTOR,J_USER1,"new_rhs");
    retval = backsolver_->Solve(lhs, ConstPtr(new_rhs));

    return retval;
  }

 ESymSolverStatus MumpsSolverInterface::Solve(Index nrhs, double *rhs_vals)
 {
   DBG_START_METH("MumpsSolverInterface::Solve", dbg_verbosity);
   DMUMPS_STRUC_C* mumps_data = (DMUMPS_STRUC_C*)mumps_ptr_;
   ESymSolverStatus retval = SYMSOLVER_SUCCESS;
   if (HaveIpData()) {
     IpData().TimingStats().LinearSystemBackSolve().Start();
   }
   for (Index i = 0; i < nrhs; i++) {
     Index offset = i * mumps_data->n;
     mumps_data->rhs = &(rhs_vals[offset]);
     mumps_data->job = 3;//solve
     Jnlst().Printf(J_MOREDETAILED, J_LINEAR_ALGEBRA,
                    "Calling MUMPS-3 for solve at cpu time %10.3f (wall %10.3f).\n", CpuTime(), WallclockTime());
     dmumps_c(mumps_data);
     Jnlst().Printf(J_MOREDETAILED, J_LINEAR_ALGEBRA,
                    "Done with MUMPS-3 for solve at cpu time %10.3f (wall %10.3f).\n", CpuTime(), WallclockTime());
     int error = mumps_data->info[0];
     if (error < 0) {
       Jnlst().Printf(J_ERROR, J_LINEAR_ALGEBRA,
                      "Error=%d returned from MUMPS in Solve.\n",
                      error);
       retval = SYMSOLVER_FATAL_ERROR;
     }
   }
   if (HaveIpData()) {
     IpData().TimingStats().LinearSystemBackSolve().End();
   }
   return retval;
 }

 void IpoptAlgorithm::UpdateHessian()
 {
   Jnlst().Printf(J_DETAILED, J_MAIN, "\n**************************************************\n");
   Jnlst().Printf(J_DETAILED, J_MAIN, "*** Update HessianMatrix for Iteration %d:", IpData().iter_count());
   Jnlst().Printf(J_DETAILED, J_MAIN, "\n**************************************************\n\n");
   hessian_updater_->UpdateHessian();
 }

  bool LoqoMuOracle::CalculateMu(Number mu_min, Number mu_max,
                                 Number& new_mu)
  {
    DBG_START_METH("LoqoMuOracle::CalculateMu",
                   dbg_verbosity);

    Number avrg_compl = IpCq().curr_avrg_compl();
    Jnlst().Printf(J_DETAILED, J_BARRIER_UPDATE,
                   "  Average complemantarity is %lf\n", avrg_compl);

    Number xi = IpCq().curr_centrality_measure();
    Jnlst().Printf(J_DETAILED, J_BARRIER_UPDATE,
                   "  Xi (distance from uniformity) is %lf\n", xi);

    //Number factor = 1.-tau_min_;   //This is the original values
    Number factor = 0.05;   //This is the value I used otherwise
    Number sigma = 0.1*pow(Min(factor*(1.-xi)/xi,2.),3.);

    Number mu = sigma*avrg_compl;
    Jnlst().Printf(J_DETAILED, J_BARRIER_UPDATE,
                   "  Barrier parameter proposed by LOQO rule is %lf\n", mu);

    // DELETEME
    char ssigma[40];
    sprintf(ssigma, " sigma=%8.2e", sigma);
    IpData().Append_info_string(ssigma);
    sprintf(ssigma, " xi=%8.2e ", IpCq().curr_centrality_measure());
    IpData().Append_info_string(ssigma);

    new_mu = Max(Min(mu_max, mu), mu_min);
    return true;
  }

  void
  AdaptiveMuUpdate::RememberCurrentPointAsAccepted()
  {
    switch (adaptive_mu_globalization_) {
    case KKT_ERROR : {
        Number curr_error = quality_function_pd_system();
        Index num_refs = (Index)refs_vals_.size();
        if (num_refs >= num_refs_max_) {
          refs_vals_.pop_front();
        }
        refs_vals_.push_back(curr_error);

        if (Jnlst().ProduceOutput(J_MOREDETAILED, J_BARRIER_UPDATE)) {
          Index num_refs = 0;
          std::list<Number>::iterator iter;
          for (iter = refs_vals_.begin(); iter != refs_vals_.end();
               iter++) {
            num_refs++;
            Jnlst().Printf(J_MOREDETAILED, J_BARRIER_UPDATE,
                           "pd system reference[%2d] = %.6e\n", num_refs, *iter);
          }
        }
      }
      break;
    case FILTER_OBJ_CONSTR : {
        /*
               Number theta = IpCq().curr_constraint_violation();
               filter_.AddEntry(IpCq().curr_f() - filter_margin_fact_*theta,
                                IpCq().curr_constraint_violation() - filter_margin_fact_*theta,
                                IpData().iter_count());
               filter_.Print(Jnlst());
        */
        filter_.AddEntry(IpCq().curr_f(),
                         IpCq().curr_constraint_violation(),
                         IpData().iter_count());
        filter_.Print(Jnlst());
      }
      break;
    case NEVER_MONOTONE_MODE : {
        // Nothing to be done
      }
      break;
    default:
      DBG_ASSERT(false && "Unknown corrector_type value.");
    }

    if (restore_accepted_iterate_) {
      // Keep pointers to this iterate so that it could be restored
      accepted_point_ = IpData().curr();
    }
  }

  void TSymLinearSolver::GiveMatrixToSolver(bool new_matrix,
      const SymMatrix& sym_A)
  {
    DBG_START_METH("TSymLinearSolver::GiveMatrixToSolver",dbg_verbosity);
    DBG_PRINT((1,"new_matrix = %d\n",new_matrix));

    double* pa = solver_interface_->GetValuesArrayPtr();
    double* atriplet;

    if (matrix_format_!=SparseSymLinearSolverInterface::Triplet_Format) {
      atriplet = new double[nonzeros_triplet_];
    }
    else {
      atriplet = pa;
    }

    //DBG_PRINT_MATRIX(3, "Aunscaled", sym_A);
    TripletHelper::FillValues(nonzeros_triplet_, sym_A, atriplet);
    if (DBG_VERBOSITY()>=3) {
      for (Index i=0; i<nonzeros_triplet_; i++) {
        DBG_PRINT((3, "KKTunscaled(%6d,%6d) = %24.16e\n", airn_[i], ajcn_[i], atriplet[i]));
      }
    }

    if (use_scaling_) {
      IpData().TimingStats().LinearSystemScaling().Start();
      DBG_ASSERT(scaling_factors_);
      if (new_matrix || just_switched_on_scaling_) {
        // only compute scaling factors if the matrix has not been
        // changed since the last call to this method
        bool retval =
          scaling_method_->ComputeSymTScalingFactors(dim_, nonzeros_triplet_,
              airn_, ajcn_,
              atriplet, scaling_factors_);
        if (!retval) {
          Jnlst().Printf(J_ERROR, J_LINEAR_ALGEBRA,
                         "Error during computation of scaling factors.\n");
          THROW_EXCEPTION(ERROR_IN_LINEAR_SCALING_METHOD, "scaling_method_->ComputeSymTScalingFactors returned false.")
        }
        // complain if not in debug mode
        if (Jnlst().ProduceOutput(J_MOREVECTOR, J_LINEAR_ALGEBRA)) {
          for (Index i=0; i<dim_; i++) {
            Jnlst().Printf(J_MOREVECTOR, J_LINEAR_ALGEBRA,
                           "scaling factor[%6d] = %22.17e\n",
                           i, scaling_factors_[i]);
          }
        }
        just_switched_on_scaling_ = false;
      }

  ESymSolverStatus WsmpSolverInterface::Solve(
    const Index* ia,
    const Index* ja,
    Index nrhs,
    double *rhs_vals)
  {
    DBG_START_METH("WsmpSolverInterface::Solve",dbg_verbosity);

    IpData().TimingStats().LinearSystemBackSolve().Start();

    // Call WSMP to solve for some right hand sides (including
    // iterative refinement)
    // ToDo: Make iterative refinement an option?
    ipfint N = dim_;
    ipfint LDB = dim_;
    ipfint NRHS = nrhs;
    ipfint NAUX = 0;
    IPARM_[1] = 4; // Forward and Backward Elimintation
    IPARM_[2] = 5; // Iterative refinement
    IPARM_[5] = 1;
    DPARM_[5] = 1e-12;

    ipfint idmy;
    double ddmy;
    F77_FUNC(wssmp,WSSMP)(&N, ia, ja, a_, &ddmy, PERM_, INVP_,
                          rhs_vals, &LDB, &NRHS, &ddmy, &NAUX,
                          &idmy, IPARM_, DPARM_);
    IpData().TimingStats().LinearSystemBackSolve().End();

    Index ierror = IPARM_[63];
    if (ierror!=0) {
      if (ierror==-102) {
        Jnlst().Printf(J_ERROR, J_LINEAR_ALGEBRA,
                       "Error: WSMP is not able to allocate sufficient amount of memory during ordering/symbolic factorization.\n");
      }
      else {
        Jnlst().Printf(J_ERROR, J_LINEAR_ALGEBRA,
                       "Error in WSMP during ordering/symbolic factorization phase.\n     Error code is %d.\n", ierror);
      }
      return SYMSOLVER_FATAL_ERROR;
    }
    Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA,
                   "Number of iterative refinement steps in WSSMP: %d\n",
                   IPARM_[5]);


    return SYMSOLVER_SUCCESS;
  }

  bool Ma86SolverInterface::IncreaseQuality()
  {
    if (control_.u >= umax_) {
      return false;
    }
    pivtol_changed_ = true;

    Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA,
                   "Increasing pivot tolerance for HSL_MA86 from %7.2e ",
                   control_.u);
    control_.u = Min(umax_, pow(control_.u,0.75));
    Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA,
                   "to %7.2e.\n",
                   control_.u);
    return true;
  }

  void IpoptAlgorithm::PrintProblemStatistics()
  {
    if (!Jnlst().ProduceOutput(J_SUMMARY, J_STATISTICS)) {
      // nothing to print
      return;
    }

    SmartPtr<const Vector> x = IpData().curr()->x();
    SmartPtr<const Vector> x_L = IpNLP().x_L();
    SmartPtr<const Vector> x_U = IpNLP().x_U();
    SmartPtr<const Matrix> Px_L = IpNLP().Px_L();
    SmartPtr<const Matrix> Px_U = IpNLP().Px_U();

    Index nx_tot, nx_only_lower, nx_both, nx_only_upper;
    calc_number_of_bounds(*IpData().curr()->x(), *IpNLP().x_L(), *IpNLP().x_U(),
                          *IpNLP().Px_L(), *IpNLP().Px_U(),
                          nx_tot, nx_only_lower, nx_both, nx_only_upper);

    Index ns_tot, ns_only_lower, ns_both, ns_only_upper;
    calc_number_of_bounds(*IpData().curr()->s(), *IpNLP().d_L(), *IpNLP().d_U(),
                          *IpNLP().Pd_L(), *IpNLP().Pd_U(),
                          ns_tot, ns_only_lower, ns_both, ns_only_upper);

    Jnlst().Printf(J_SUMMARY, J_STATISTICS,
                   "Total number of variables............................: %8d\n",nx_tot);
    Jnlst().Printf(J_SUMMARY, J_STATISTICS,
                   "                     variables with only lower bounds: %8d\n",
                   nx_only_lower);
    Jnlst().Printf(J_SUMMARY, J_STATISTICS,
                   "                variables with lower and upper bounds: %8d\n",nx_both);
    Jnlst().Printf(J_SUMMARY, J_STATISTICS,
                   "                     variables with only upper bounds: %8d\n",
                   nx_only_upper);
    Jnlst().Printf(J_SUMMARY, J_STATISTICS,
                   "Total number of equality constraints.................: %8d\n",
                   IpData().curr()->y_c()->Dim());
    Jnlst().Printf(J_SUMMARY, J_STATISTICS,
                   "Total number of inequality constraints...............: %8d\n",ns_tot);
    Jnlst().Printf(J_SUMMARY, J_STATISTICS,
                   "        inequality constraints with only lower bounds: %8d\n",
                   ns_only_lower);
    Jnlst().Printf(J_SUMMARY, J_STATISTICS,
                   "   inequality constraints with lower and upper bounds: %8d\n",ns_both);
    Jnlst().Printf(J_SUMMARY, J_STATISTICS,
                   "        inequality constraints with only upper bounds: %8d\n\n",
                   ns_only_upper);
  }

  bool Ma27TSolverInterface::IncreaseQuality()
  {
    DBG_START_METH("Ma27TSolverInterface::IncreaseQuality",dbg_verbosity);
    if (pivtol_ == pivtolmax_) {
      return false;
    }
    pivtol_changed_ = true;

    Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA,
                   "Indreasing pivot tolerance for MA27 from %7.2e ",
                   pivtol_);
    pivtol_ = Min(pivtolmax_, pow(pivtol_,0.75));
    Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA,
                   "to %7.2e.\n",
                   pivtol_);
    return true;
  }

  void IpoptAlgorithm::ComputeFeasibilityMultipliers()
  {
    DBG_START_METH("IpoptAlgorithm::ComputeFeasibilityMultipliers",
                   dbg_verbosity);
    DBG_ASSERT(IpCq().IsSquareProblem());

    // if we don't have an object for computing least square
    // multipliers we don't compute them
    if (IsNull(eq_multiplier_calculator_)) {
      Jnlst().Printf(J_WARNING, J_SOLUTION,
                     "This is a square problem, but multipliers cannot be recomputed at solution, since no eq_mult_calculator object is available in IpoptAlgorithm\n");
      return;
    }

    SmartPtr<IteratesVector> iterates = IpData().curr()->MakeNewContainer();
    SmartPtr<Vector> tmp = iterates->z_L()->MakeNew();
    tmp->Set(0.);
    iterates->Set_z_L(*tmp);
    tmp = iterates->z_U()->MakeNew();
    tmp->Set(0.);
    iterates->Set_z_U(*tmp);
    tmp = iterates->v_L()->MakeNew();
    tmp->Set(0.);
    iterates->Set_v_L(*tmp);
    tmp = iterates->v_U()->MakeNew();
    tmp->Set(0.);
    iterates->Set_v_U(*tmp);
    SmartPtr<Vector> y_c = iterates->y_c()->MakeNew();
    SmartPtr<Vector> y_d = iterates->y_d()->MakeNew();
    IpData().set_trial(iterates);
    IpData().AcceptTrialPoint();
    bool retval = eq_multiplier_calculator_->CalculateMultipliers(*y_c, *y_d);
    if (retval) {
      //Check if following line is really necessary
      iterates = IpData().curr()->MakeNewContainer();
      iterates->Set_y_c(*y_c);
      iterates->Set_y_d(*y_d);
      IpData().set_trial(iterates);
      IpData().AcceptTrialPoint();
    }
    else {
      Jnlst().Printf(J_WARNING, J_SOLUTION,
                     "Cannot recompute multipliers for feasibility problem.  Error in eq_mult_calculator\n");
    }
  }