Number alpha_dual_aff = IpCq().dual_frac_to_the_bound(1.0,
*step->z_L(),
*step->z_U(),
*step->v_L(),
*step->v_U());
Jnlst().Printf(J_DETAILED, J_BARRIER_UPDATE,
" The affine maximal step sizes are\n"
" alpha_primal_aff = %23.16e\n"
" alpha_dual_aff = %23.16e\n",
alpha_primal_aff,
alpha_dual_aff);
// now compute the average complementarity at the affine step
// ToDo shoot for mu_min instead of 0?
Number mu_aff = CalculateAffineMu(alpha_primal_aff, alpha_dual_aff, *step);
Jnlst().Printf(J_DETAILED, J_BARRIER_UPDATE,
" The average complementariy at the affine step is %23.16e\n",
mu_aff);
// get the current average complementarity
Number mu_curr = IpCq().curr_avrg_compl();
Jnlst().Printf(J_DETAILED, J_BARRIER_UPDATE,
" The average complementariy at the current point is %23.16e\n",
mu_curr);
DBG_ASSERT(mu_curr>0.);
// Apply Mehrotra's rule
Number sigma = pow((mu_aff/mu_curr),3);
// Make sure, sigma is not too large
sigma = Min(sigma, sigma_max_);
Number mu = sigma*mu_curr;
// Store the affine search direction (in case it is needed in the
// line search for a corrector step)
IpData().set_delta_aff(step);
IpData().SetHaveAffineDeltas(true);
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, mu_max), mu_min);
return true;
}
Number ProbingMuOracle::CalculateAffineMu
(
Number alpha_primal,
Number alpha_dual,
const IteratesVector& step)
{
// Get the current values of the slack variables and bound multipliers
SmartPtr<const Vector> slack_x_L = IpCq().curr_slack_x_L();
SmartPtr<const Vector> slack_x_U = IpCq().curr_slack_x_U();
SmartPtr<const Vector> slack_s_L = IpCq().curr_slack_s_L();
SmartPtr<const Vector> slack_s_U = IpCq().curr_slack_s_U();
SmartPtr<const Vector> z_L = IpData().curr()->z_L();
SmartPtr<const Vector> z_U = IpData().curr()->z_U();
SmartPtr<const Vector> v_L = IpData().curr()->v_L();
SmartPtr<const Vector> v_U = IpData().curr()->v_U();
SmartPtr<Vector> tmp_slack;
SmartPtr<Vector> tmp_mult;
SmartPtr<const Matrix> P;
Index ncomp = 0;
Number sum =0.;
// For each combination of slack and multiplier, compute the new
// values and their dot products.
// slack_x_L
if (slack_x_L->Dim()>0) {
ncomp += slack_x_L->Dim();
P = IpNLP().Px_L();
tmp_slack = slack_x_L->MakeNew();
tmp_slack->Copy(*slack_x_L);
P->TransMultVector(alpha_primal, *step.x(), 1.0, *tmp_slack);
tmp_mult = z_L->MakeNew();
tmp_mult->Copy(*z_L);