Skip to content

Commit 84689ae

Browse files
filikatfilikat
committed
Split omega into various components
1 parent b7b7059 commit 84689ae

File tree

5 files changed

+72
-20
lines changed

5 files changed

+72
-20
lines changed

‎highs/ipm/hipo/ipm/IpmData.h‎

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -21,6 +21,9 @@ struct IpmIterData {
2121
Int num_small_prod = 0;
2222
Int num_large_prod = 0;
2323
double omega = 0.0;
24+
double omega_null = 0.0;
25+
double omega_image = 0.0;
26+
double omega_rest = 0.0;
2427
double nw_back_err = 0.0;
2528
double cw_back_err = 0.0;
2629
Int large_components_cw = 0;

‎highs/ipm/hipo/ipm/Iterate.cpp‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -658,7 +658,7 @@ void Iterate::residuals6x6(const NewtonDir& d) {
658658
ires.r1 = res.r1;
659659
model->A().alphaProductPlusY(-1.0, dx, ires.r1);
660660
for (Int i = 0; i < m; ++i) {
661-
ires.r1[i] -= Rd[i] * dy[i];
661+
ires.r1[i] -= std::abs(Rd[i]) * dy[i];
662662
}
663663

664664
// ires2 = res2 - dx + dxl
@@ -684,7 +684,7 @@ void Iterate::residuals6x6(const NewtonDir& d) {
684684
model->A().alphaProductPlusY(-1.0, dy, ires.r4, true);
685685
for (Int i = 0; i < n; ++i) {
686686
double reg_p = Rp ? Rp[i] : regul.primal;
687-
ires.r4[i] += reg_p * dx[i];
687+
ires.r4[i] += std::abs(reg_p) * dx[i];
688688
}
689689

690690
// ires5 = res5 - zl * dxl - xl * dzl

‎highs/ipm/hipo/ipm/Refine.cpp‎

Lines changed: 49 additions & 12 deletions
Original file line numberDiff line numberDiff line change
@@ -16,10 +16,12 @@ void Solver::refine(NewtonDir& delta) {
1616
info_.residual_time += clock.stop();
1717

1818
clock.start();
19-
double omega = computeOmega(delta);
19+
OmegaValues omegavals = computeOmega(delta);
20+
double omega = omegavals.compute();
2021
info_.omega_time += clock.stop();
2122

2223
double old_omega{};
24+
OmegaValues old_values;
2325

2426
for (Int iter = 0; iter < kMaxIterRefine; ++iter) {
2527
if (omega < kTolRefine) break;
@@ -35,21 +37,30 @@ void Solver::refine(NewtonDir& delta) {
3537
info_.residual_time += clock.stop();
3638

3739
old_omega = omega;
40+
old_values = omegavals;
3841

3942
clock.start();
40-
omega = computeOmega(temp);
43+
omegavals = computeOmega(temp);
44+
omega = omegavals.compute();
4145
info_.omega_time += clock.stop();
4246

4347
if (omega < old_omega) {
4448
delta = temp;
4549
} else {
4650
omega = old_omega;
51+
omegavals = old_values;
4752
break;
4853
}
4954
}
5055

5156
// save worst residual seen in this ipm iteration
5257
it_->data.back().omega = std::max(it_->data.back().omega, omega);
58+
59+
double null, image, rest;
60+
omegavals.extract(null, image, rest);
61+
it_->data.back().omega_null = std::max(it_->data.back().omega_null, null);
62+
it_->data.back().omega_image = std::max(it_->data.back().omega_image, image);
63+
it_->data.back().omega_rest = std::max(it_->data.back().omega_rest, rest);
5364
}
5465

5566
static void updateOmega(double tau, double& omega1, double& omega2, double num,
@@ -63,7 +74,33 @@ static void updateOmega(double tau, double& omega1, double& omega2, double num,
6374
}
6475
}
6576

66-
double Solver::computeOmega(const NewtonDir& delta) const {
77+
double Solver::OmegaValues::compute() const {
78+
// Given the omega values of each block of equations, compute the omega value
79+
// of the whole system.
80+
double v1{}, v2{};
81+
for (Int i = 0; i < 6; ++i) {
82+
v1 = std::max(v1, omega_1[i]);
83+
v2 = std::max(v2, omega_2[i]);
84+
}
85+
return v1 + v2;
86+
}
87+
88+
void Solver::OmegaValues::extract(double& null, double& image,
89+
double& rest) const {
90+
// Given the omega values of each block of equations, compute the omega values
91+
// corresponding to the equations involging A ("null space error"), A^T
92+
// ("image space error") and the rest.
93+
94+
null = omega_1[0] + omega_2[0];
95+
image = omega_1[3] + omega_2[3];
96+
97+
rest = omega_1[1] + omega_2[1];
98+
rest = std::max(rest, omega_1[2] + omega_2[2]);
99+
rest = std::max(rest, omega_1[4] + omega_2[4]);
100+
rest = std::max(rest, omega_1[5] + omega_2[5]);
101+
}
102+
103+
Solver::OmegaValues Solver::computeOmega(const NewtonDir& delta) const {
67104
// Evaluate iterative refinement progress. Based on "Solving sparse linear
68105
// systems with sparse backward error", Arioli, Demmel, Duff.
69106

@@ -91,12 +128,12 @@ double Solver::computeOmega(const NewtonDir& delta) const {
91128
// rhs is in it_->res
92129
// residual of linear system is in it_->ires
93130

131+
OmegaValues w;
132+
94133
// tau_i =
95134
// tau_const * (inf_norm_rows(big 6x6 matrix)_i * inf_norm_delta + |rhs_i|)
96135

97136
const double tau_const = 1000.0 * (5 * n_ + m_) * 1e-16;
98-
double omega_1{}, omega_2{};
99-
100137
assert(it_->Rd);
101138

102139
// First block
@@ -106,7 +143,7 @@ double Solver::computeOmega(const NewtonDir& delta) const {
106143
inf_norm_delta +
107144
std::abs(it_->res.r1[i]));
108145
updateOmega(
109-
tau, omega_1, omega_2,
146+
tau, w.omega_1[0], w.omega_2[0],
110147
// numerator
111148
std::abs(it_->ires.r1[i]),
112149
// denominators
@@ -120,7 +157,7 @@ double Solver::computeOmega(const NewtonDir& delta) const {
120157
if (model_.hasLb(i)) {
121158
const double tau =
122159
tau_const * (1.0 * inf_norm_delta + std::abs(it_->res.r2[i]));
123-
updateOmega(tau, omega_1, omega_2,
160+
updateOmega(tau, w.omega_1[1], w.omega_2[1],
124161
// numerator
125162
std::abs(it_->ires.r2[i]),
126163
// denominators
@@ -136,7 +173,7 @@ double Solver::computeOmega(const NewtonDir& delta) const {
136173
if (model_.hasUb(i)) {
137174
const double tau =
138175
tau_const * (1.0 * inf_norm_delta + std::abs(it_->res.r3[i]));
139-
updateOmega(tau, omega_1, omega_2,
176+
updateOmega(tau, w.omega_1[2], w.omega_2[2],
140177
// numerator
141178
std::abs(it_->ires.r3[i]),
142179
// denominators
@@ -160,7 +197,7 @@ double Solver::computeOmega(const NewtonDir& delta) const {
160197
if (model_.hasUb(i)) denom1 += std::abs(delta.zu[i]);
161198
denom1 += std::abs(reg_p) * std::abs(delta.x[i]);
162199

163-
updateOmega(tau, omega_1, omega_2,
200+
updateOmega(tau, w.omega_1[3], w.omega_2[3],
164201
// numerator
165202
std::abs(it_->ires.r4[i]),
166203
// denominators
@@ -177,7 +214,7 @@ double Solver::computeOmega(const NewtonDir& delta) const {
177214
std::abs(it_->res.r5[i]));
178215

179216
updateOmega(
180-
tau, omega_1, omega_2,
217+
tau, w.omega_1[4], w.omega_2[4],
181218
// numerator
182219
std::abs(it_->ires.r5[i]),
183220
// denominators
@@ -197,7 +234,7 @@ double Solver::computeOmega(const NewtonDir& delta) const {
197234
std::abs(it_->res.r6[i]));
198235

199236
updateOmega(
200-
tau, omega_1, omega_2,
237+
tau, w.omega_1[5], w.omega_2[5],
201238
// numerator
202239
std::abs(it_->ires.r6[i]),
203240
// denominators
@@ -208,7 +245,7 @@ double Solver::computeOmega(const NewtonDir& delta) const {
208245
}
209246
}
210247

211-
return omega_1 + omega_2;
248+
return w;
212249
}
213250

214251
} // namespace hipo

‎highs/ipm/hipo/ipm/Solver.cpp‎

Lines changed: 8 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -385,7 +385,7 @@ void Solver::recoverDirection(NewtonDir& delta, const Residuals& rhs) const {
385385
}
386386

387387
// not sure if this has any effect, but IPX uses it
388-
std::vector<double> Atdy(n_);
388+
/*std::vector<double> Atdy(n_);
389389
model_.A().alphaProductPlusY(1.0, delta.y, Atdy, true);
390390
for (Int i = 0; i < n_; ++i) {
391391
if (model_.hasLb(i) || model_.hasUb(i)) {
@@ -400,7 +400,7 @@ void Solver::recoverDirection(NewtonDir& delta, const Residuals& rhs) const {
400400
delta.zu[i] = -res4[i] + delta.zl[i] + Atdy[i];
401401
}
402402
}
403-
}
403+
}*/
404404
}
405405

406406
double Solver::stepToBoundary(const std::vector<double>& x,
@@ -1111,9 +1111,10 @@ void Solver::printHeader() const {
11111111
if (!options_.timeless_log) log_stream << " time";
11121112

11131113
if (logH_.debug(1)) {
1114-
log_stream << " alpha p/d sigma af/co cor solv"
1115-
<< " static reg p/d minT maxT"
1116-
<< " (xj * zj / mu)_range_&_num max_res";
1114+
log_stream
1115+
<< " alpha p/d sigma af/co cor solv static reg p/d "
1116+
"minT maxT (xj * zj / mu)_range_&_num max_res null "
1117+
"image";
11171118
}
11181119

11191120
log_stream << "\n";
@@ -1152,6 +1153,8 @@ void Solver::printOutput() const {
11521153
log_stream << " " << integer(data.num_small_prod, 4);
11531154
log_stream << " " << integer(data.num_large_prod, 4);
11541155
log_stream << " " << sci(data.omega, 9, 1);
1156+
log_stream << " " << sci(data.omega_null, 8, 1);
1157+
log_stream << " " << sci(data.omega_image, 8, 1);
11551158
}
11561159

11571160
log_stream << "\n";

‎highs/ipm/hipo/ipm/Solver.h‎

Lines changed: 10 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -202,8 +202,17 @@ class Solver {
202202
// ===================================================================================
203203
// Functions for iterative refinement on the large 6x6 system
204204
// ===================================================================================
205+
struct OmegaValues {
206+
// store the omega values for the 6 blocks of equations. See "Solving sparse
207+
// linear systems with sparse backward error", Arioli, Demmel, Duff, for an
208+
// explanation of what omega is.
209+
double omega_1[6] = {};
210+
double omega_2[6] = {};
211+
double compute() const;
212+
void extract(double& null, double& image, double& rest) const;
213+
};
205214
void refine(NewtonDir& delta);
206-
double computeOmega(const NewtonDir& delta) const;
215+
OmegaValues computeOmega(const NewtonDir& delta) const;
207216

208217
// ===================================================================================
209218
// Steps to boundary are computed so that

0 commit comments

Comments
 (0)