Skip to content

Commit 2c24c64

Browse files
filikatfilikat
committed
Remove previous refinement
1 parent 63074e4 commit 2c24c64

File tree

13 files changed

+50
-273
lines changed

13 files changed

+50
-273
lines changed

‎highs/ipm/hipo/factorhighs/DataCollector.cpp‎

Lines changed: 0 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -203,14 +203,6 @@ void DataCollector::printTimes(const Log& log) const {
203203
<< fix(times[kTimeSolveSolve_sparse], 8, 4) << "\n";
204204
log_stream << "\t\tswap: " << fix(times[kTimeSolveSolve_swap], 8, 4)
205205
<< "\n";
206-
log_stream << "\tResidual: "
207-
<< fix(times[kTimeSolveResidual], 8, 4) << " ("
208-
<< fix(times[kTimeSolveResidual] / times[kTimeSolve] * 100, 4, 1)
209-
<< "%)\n";
210-
log_stream << "\tOmega: "
211-
<< fix(times[kTimeSolveOmega], 8, 4) << " ("
212-
<< fix(times[kTimeSolveOmega] / times[kTimeSolve] * 100, 4, 1)
213-
<< "%)\n";
214206
#endif
215207
log_stream << "----------------------------------------------------\n";
216208

‎highs/ipm/hipo/factorhighs/FactorHiGHS.cpp‎

Lines changed: 1 addition & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -48,8 +48,6 @@ Int FHsolver::factorise(const Symbolic& S, const std::vector<Int>& rows,
4848
return fact_obj.run(N_);
4949
}
5050

51-
Int FHsolver::solve(std::vector<double>& x, Int* solve_count, double* omega) {
52-
return N_.solve(x, solve_count, omega);
53-
}
51+
Int FHsolver::solve(std::vector<double>& x) { return N_.solve(x); }
5452

5553
} // namespace hipo

‎highs/ipm/hipo/factorhighs/FactorHiGHS.h‎

Lines changed: 2 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -91,13 +91,8 @@ class FHsolver {
9191
const std::vector<Int>& ptr, const std::vector<double>& vals);
9292

9393
// Perform solve phase with rhs given by x, which is overwritten with the
94-
// solution. solve_count returns the number of solves performed during the
95-
// phase (including refinement). omega returns the final residual after
96-
// refinement.
97-
// For now refinement is performed automatically as part of solve,
98-
// this will change in the future.
99-
Int solve(std::vector<double>& x, Int* solve_count = nullptr,
100-
double* omega = nullptr);
94+
// solution.
95+
Int solve(std::vector<double>& x);
10196

10297
// If multiple factorisation are performed, call newIter() before each
10398
// factorisation. This is used only to collect data for debugging, if

‎highs/ipm/hipo/factorhighs/FactorHiGHSSettings.h‎

Lines changed: 0 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -47,10 +47,6 @@ const Int kBlockParallelThreshold = 5;
4747
// regularisation
4848
const double kDynamicDiagCoeff = 1e-24;
4949

50-
// refinement
51-
const Int kMaxRefinementIter = 3;
52-
const double kRefinementTolerance = 1e-12;
53-
5450
// metis
5551
const Int kMetisSeed = 42;
5652

‎highs/ipm/hipo/factorhighs/Factorise.cpp‎

Lines changed: 4 additions & 21 deletions
Original file line numberDiff line numberDiff line change
@@ -74,29 +74,17 @@ Factorise::Factorise(const Symbolic& S, const std::vector<Int>& rowsA,
7474
min_diag_ = std::min(min_diag_, val);
7575
}
7676

77-
// infinity norm of columns of A
78-
inf_norm_cols_.assign(n_, 0.0);
79-
for (Int col = 0; col < n_; ++col) {
80-
for (Int el = ptrA_[col]; el < ptrA_[col + 1]; ++el) {
81-
Int row = rowsA_[el];
82-
double val = valA_[el];
83-
inf_norm_cols_[col] = std::max(inf_norm_cols_[col], std::abs(val));
84-
if (row != col)
85-
inf_norm_cols_[row] = std::max(inf_norm_cols_[row], std::abs(val));
86-
}
87-
}
88-
8977
// one norm of columns of A
90-
one_norm_cols_.assign(n_, 0.0);
78+
std::vector<double> one_norm_cols(n_, 0.0);
9179
for (Int col = 0; col < n_; ++col) {
9280
for (Int el = ptrA_[col]; el < ptrA_[col + 1]; ++el) {
9381
Int row = rowsA_[el];
9482
double val = valA_[el];
95-
one_norm_cols_[col] += std::abs(val);
96-
if (row != col) one_norm_cols_[row] += std::abs(val);
83+
one_norm_cols[col] += std::abs(val);
84+
if (row != col) one_norm_cols[row] += std::abs(val);
9785
}
9886
}
99-
A_norm1_ = *std::max_element(one_norm_cols_.begin(), one_norm_cols_.end());
87+
A_norm1_ = *std::max_element(one_norm_cols.begin(), one_norm_cols.end());
10088

10189
data_.setNorms(A_norm1_, max_diag_);
10290
}
@@ -421,11 +409,6 @@ bool Factorise::run(Numeric& num) {
421409
num.total_reg_ = std::move(total_reg_);
422410
num.swaps_ = std::move(swaps_);
423411
num.pivot_2x2_ = std::move(pivot_2x2_);
424-
num.ptrA_ = std::move(ptrA_);
425-
num.rowsA_ = std::move(rowsA_);
426-
num.valA_ = std::move(valA_);
427-
num.one_norm_cols_ = std::move(one_norm_cols_);
428-
num.inf_norm_cols_ = std::move(inf_norm_cols_);
429412
num.data_ = &data_;
430413

431414
#if HIPO_TIMING_LEVEL >= 1

‎highs/ipm/hipo/factorhighs/Factorise.h‎

Lines changed: 0 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -51,8 +51,6 @@ class Factorise {
5151
double max_diag_{};
5252
double min_diag_{};
5353
double A_norm1_{};
54-
std::vector<double> one_norm_cols_{};
55-
std::vector<double> inf_norm_cols_{};
5654

5755
// regularisation
5856
std::vector<double> total_reg_{};

‎highs/ipm/hipo/factorhighs/Numeric.cpp‎

Lines changed: 1 addition & 173 deletions
Original file line numberDiff line numberDiff line change
@@ -13,8 +13,7 @@
1313

1414
namespace hipo {
1515

16-
Int Numeric::solve(std::vector<double>& x, Int* solve_count,
17-
double* omega) const {
16+
Int Numeric::solve(std::vector<double>& x) const {
1817
// Return the number of solves performed
1918

2019
if (!sn_columns_ || !S_) return kRetInvalidPointer;
@@ -31,13 +30,8 @@ Int Numeric::solve(std::vector<double>& x, Int* solve_count,
3130
#if HIPO_TIMING_LEVEL >= 2
3231
Clock clock_fine{};
3332
#endif
34-
3533
// permute rhs
3634
permuteVectorInverse(x, S_->iperm());
37-
38-
// make a copy of permuted rhs, for refinement
39-
const std::vector<double> rhs(x);
40-
4135
#if HIPO_TIMING_LEVEL >= 2
4236
if (data_) data_->sumTime(kTimeSolvePrepare, clock_fine.stop());
4337
clock_fine.start();
@@ -52,16 +46,11 @@ Int Numeric::solve(std::vector<double>& x, Int* solve_count,
5246
if (data_) data_->sumTime(kTimeSolveSolve, clock_fine.stop());
5347
#endif
5448

55-
// iterative refinement
56-
auto refine_data = refine(rhs, x);
57-
5849
#if HIPO_TIMING_LEVEL >= 2
5950
clock_fine.start();
6051
#endif
61-
6252
// unpermute solution
6353
permuteVector(x, S_->iperm());
64-
6554
#if HIPO_TIMING_LEVEL >= 2
6655
if (data_) data_->sumTime(kTimeSolvePrepare, clock_fine.stop());
6756
#endif
@@ -70,168 +59,7 @@ Int Numeric::solve(std::vector<double>& x, Int* solve_count,
7059
if (data_) data_->sumTime(kTimeSolve, clock.stop());
7160
#endif
7261

73-
if (solve_count) *solve_count = refine_data.first + 1;
74-
if (omega) *omega = refine_data.second;
75-
7662
return kRetOk;
7763
}
7864

79-
std::vector<double> Numeric::residual(const std::vector<double>& rhs,
80-
const std::vector<double>& x) const {
81-
// Compute the residual rhs - A * x - Reg * x
82-
std::vector<double> res(rhs);
83-
symProduct(ptrA_, rowsA_, valA_, x, res, -1.0);
84-
for (Int i = 0; i < x.size(); ++i) res[i] -= total_reg_[i] * x[i];
85-
86-
return res;
87-
}
88-
89-
std::vector<double> Numeric::residualQuad(const std::vector<double>& rhs,
90-
const std::vector<double>& x) const {
91-
std::vector<HighsCDouble> res(rhs.size());
92-
for (Int i = 0; i < rhs.size(); ++i) res[i] = rhs[i];
93-
94-
symProductQuad(ptrA_, rowsA_, valA_, x, res, -1.0);
95-
96-
for (Int i = 0; i < x.size(); ++i)
97-
res[i] -= (HighsCDouble)total_reg_[i] * (HighsCDouble)x[i];
98-
99-
std::vector<double> result(res.size());
100-
for (Int i = 0; i < res.size(); ++i) {
101-
result[i] = (double)res[i];
102-
}
103-
104-
return result;
105-
}
106-
107-
std::pair<Int, double> Numeric::refine(const std::vector<double>& rhs,
108-
std::vector<double>& x) const {
109-
// Return the number of solver performed
110-
111-
double old_omega{};
112-
Int solves_counter{};
113-
114-
#if HIPO_TIMING_LEVEL >= 2
115-
Clock clock{};
116-
#endif
117-
118-
// compute residual
119-
std::vector<double> res = residualQuad(rhs, x);
120-
121-
#if HIPO_TIMING_LEVEL >= 2
122-
if (data_) data_->sumTime(kTimeSolveResidual, clock.stop());
123-
clock.start();
124-
#endif
125-
126-
double omega = computeOmega(rhs, x, res);
127-
128-
#if HIPO_TIMING_LEVEL >= 2
129-
if (data_) data_->sumTime(kTimeSolveOmega, clock.stop());
130-
#endif
131-
132-
// if(log_) log_->printDevVerbose(" # start %.2e\n", omega);
133-
134-
Int iter = 0;
135-
for (; iter < kMaxRefinementIter; ++iter) {
136-
// termination criterion
137-
if (omega < kRefinementTolerance) break;
138-
139-
#if HIPO_TIMING_LEVEL >= 2
140-
clock.start();
141-
#endif
142-
143-
// compute correction
144-
SH_->forwardSolve(res);
145-
SH_->diagSolve(res);
146-
SH_->backwardSolve(res);
147-
++solves_counter;
148-
149-
#if HIPO_TIMING_LEVEL >= 2
150-
if (data_) data_->sumTime(kTimeSolveSolve, clock.stop());
151-
clock.start();
152-
#endif
153-
154-
// add correction
155-
std::vector<double> temp(x);
156-
vectorAdd(temp, res);
157-
158-
#if HIPO_TIMING_LEVEL >= 2
159-
if (data_) data_->sumTime(kTimeSolvePrepare, clock.stop());
160-
clock.start();
161-
#endif
162-
163-
// compute new residual
164-
res = residualQuad(rhs, temp);
165-
166-
#if HIPO_TIMING_LEVEL >= 2
167-
if (data_) data_->sumTime(kTimeSolveResidual, clock.stop());
168-
clock.start();
169-
#endif
170-
171-
old_omega = omega;
172-
omega = computeOmega(rhs, temp, res);
173-
174-
// if(log_) log_->printDevVerbose(" # refine %.2e\n", omega);
175-
176-
#if HIPO_TIMING_LEVEL >= 2
177-
if (data_) data_->sumTime(kTimeSolveOmega, clock.stop());
178-
#endif
179-
180-
if (omega < old_omega) {
181-
x = temp;
182-
} else {
183-
omega = old_omega;
184-
// if(log_) log_->printDevVerbose(" ## reject\n");
185-
break;
186-
}
187-
}
188-
189-
return {solves_counter, omega};
190-
}
191-
192-
double Numeric::computeOmega(const std::vector<double>& b,
193-
const std::vector<double>& x,
194-
const std::vector<double>& res) const {
195-
// Termination of iterative refinement based on "Solving sparse linear systems
196-
// with sparse backward error", Arioli, Demmel, Duff.
197-
198-
const Int n = x.size();
199-
200-
// infinity norm of x
201-
const double inf_norm_x = infNorm(x);
202-
203-
// |A|*|x|
204-
std::vector<double> abs_prod(n);
205-
for (Int col = 0; col < n; ++col) {
206-
for (Int el = ptrA_[col]; el < ptrA_[col + 1]; ++el) {
207-
Int row = rowsA_[el];
208-
double val = valA_[el];
209-
abs_prod[row] += std::abs(val) * std::abs(x[col]);
210-
if (row != col) abs_prod[col] += std::abs(val) * std::abs(x[row]);
211-
}
212-
}
213-
214-
double omega_1{};
215-
double omega_2{};
216-
217-
for (Int i = 0; i < n; ++i) {
218-
// threshold 1000 * n * eps * (||Ai|| * ||x|| + |bi|)
219-
double tau =
220-
1000.0 * n * 1e-16 * (inf_norm_cols_[i] * inf_norm_x + std::abs(b[i]));
221-
222-
if (abs_prod[i] + std::abs(b[i]) > tau) {
223-
// case 1, denominator is large enough
224-
double omega = std::abs(res[i]) / (abs_prod[i] + std::abs(b[i]));
225-
omega_1 = std::max(omega_1, omega);
226-
} else {
227-
// case 2, denominator would be small, change it
228-
double omega =
229-
std::abs(res[i]) / (abs_prod[i] + one_norm_cols_[i] * inf_norm_x);
230-
omega_2 = std::max(omega_2, omega);
231-
}
232-
}
233-
234-
return omega_1 + omega_2;
235-
}
236-
23765
} // namespace hipo

‎highs/ipm/hipo/factorhighs/Numeric.h‎

Lines changed: 2 additions & 23 deletions
Original file line numberDiff line numberDiff line change
@@ -31,14 +31,6 @@ class Numeric {
3131
// object to handle solve phase in different formats
3232
mutable std::unique_ptr<SolveHandler> SH_;
3333

34-
// lower triangle of original matrix, permuted
35-
std::vector<Int> rowsA_{};
36-
std::vector<Int> ptrA_{};
37-
std::vector<double> valA_{};
38-
39-
// norms of columns of matrix
40-
std::vector<double> inf_norm_cols_, one_norm_cols_;
41-
4234
DataCollector* data_ = nullptr;
4335

4436
friend class Factorise;
@@ -47,21 +39,8 @@ class Numeric {
4739
// dynamic regularisation applied to the matrix
4840
std::vector<double> total_reg_{};
4941

50-
// Full solve with refinement
51-
// Return also number of solves and final residual
52-
Int solve(std::vector<double>& x, Int* solve_count = nullptr,
53-
double* omega = nullptr) const;
54-
55-
// Iterative refinement
56-
std::pair<Int, double> refine(const std::vector<double>& rhs,
57-
std::vector<double>& x) const;
58-
std::vector<double> residual(const std::vector<double>& rhs,
59-
const std::vector<double>& x) const;
60-
std::vector<double> residualQuad(const std::vector<double>& rhs,
61-
const std::vector<double>& x) const;
62-
double computeOmega(const std::vector<double>& b,
63-
const std::vector<double>& x,
64-
const std::vector<double>& res) const;
42+
// Full solve
43+
Int solve(std::vector<double>& x) const;
6544
};
6645

6746
} // namespace hipo

‎highs/ipm/hipo/factorhighs/Timing.h‎

Lines changed: 0 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -30,8 +30,6 @@ enum TimeItems {
3030
kTimeSolveSolve_dense, // TIMING_LEVEL 2
3131
kTimeSolveSolve_sparse, // TIMING_LEVEL 2
3232
kTimeSolveSolve_swap, // TIMING_LEVEL 2
33-
kTimeSolveResidual, // TIMING_LEVEL 2
34-
kTimeSolveOmega, // TIMING_LEVEL 2
3533
kTimeBlasStart, //
3634
kTimeBlas_copy = kTimeBlasStart, // TIMING_LEVEL 3
3735
kTimeBlas_axpy, // TIMING_LEVEL 3

0 commit comments

Comments
 (0)