Add refinement on the 6x6 system

Author: filikat

highs/ipm/hipo/ipm/Iterate.cpp

@@ -8,6 +8,24 @@ namespace hipo {
 NewtonDir::NewtonDir(Int m, Int n)
     : x(n, 0.0), y(m, 0.0), xl(n, 0.0), xu(n, 0.0), zl(n, 0.0), zu(n, 0.0) {}
 
+void NewtonDir::clear() {
+  std::fill(x.begin(), x.end(), 0.0);
+  std::fill(y.begin(), y.end(), 0.0);
+  std::fill(xl.begin(), xl.end(), 0.0);
+  std::fill(xu.begin(), xu.end(), 0.0);
+  std::fill(zl.begin(), zl.end(), 0.0);
+  std::fill(zu.begin(), zu.end(), 0.0);
+}
+
+void NewtonDir::add(const NewtonDir& d) {
+  vectorAdd(x, d.x);
+  vectorAdd(y, d.y);
+  vectorAdd(xl, d.xl);
+  vectorAdd(xu, d.xu);
+  vectorAdd(zl, d.zl);
+  vectorAdd(zu, d.zu);
+}
+
 Iterate::Iterate(const Model& model_input, Regularisation& r)
     : model{&model_input}, delta(model->m(), model->n()), regul{r} {
   clearIter();
@@ -257,14 +275,7 @@ void Iterate::clearRes() {
   res.r5.assign(model->n(), 0.0);
   res.r6.assign(model->n(), 0.0);
 }
-void Iterate::clearDir() {
-  delta.x.assign(model->n(), 0.0);
-  delta.xl.assign(model->n(), 0.0);
-  delta.xu.assign(model->n(), 0.0);
-  delta.y.assign(model->m(), 0.0);
-  delta.zl.assign(model->n(), 0.0);
-  delta.zu.assign(model->n(), 0.0);
-}
+void Iterate::clearDir() { delta.clear(); }
 void Iterate::clearIres() {
   ires.r1.assign(model->m(), 0.0);
   ires.r2.assign(model->n(), 0.0);

highs/ipm/hipo/ipm/Iterate.h

@@ -20,6 +20,8 @@ struct NewtonDir {
   std::vector<double> zu{};
 
   NewtonDir(Int m, Int n);
+  void clear();
+  void add(const NewtonDir& d);
 };
 
 struct Residuals {

highs/ipm/hipo/ipm/Parameters.h

@@ -37,6 +37,10 @@ const Int kMinNumberSn = 10;
 const double kDenseColThresh = 0.5;
 const Int kMinRowsForDensity = 2000;
 
+// parameters for iterative refinement
+const Int kMaxIterRefine = 3;
+const double kTolRefine = 1e-12;
+
 // static regularisation
 struct Regularisation {
   double primal = 1e-12;

highs/ipm/hipo/ipm/Refine.cpp

@@ -5,10 +5,36 @@
 namespace hipo {
 
 void Solver::refine(NewtonDir& delta) {
+  NewtonDir correction(m_, n_);
+  NewtonDir temp(m_, n_);
+
   // compute the residuals of the linear system in it_->ires
   it_->residuals6x6(delta);
 
   double omega = computeOmega(delta);
+  double old_omega{};
+
+  for (Int iter = 0; iter < kMaxIterRefine; ++iter) {
+    if (omega < kTolRefine) break;
+
+    correction.clear();
+    solveNewtonSystem(correction, it_->ires);
+
+    temp = delta;
+    temp.add(correction);
+
+    it_->residuals6x6(temp);
+
+    old_omega = omega;
+    omega = computeOmega(temp);
+
+    if (omega < old_omega) {
+      delta = temp;
+    } else {
+      omega = old_omega;
+      break;
+    }
+  }
 }
 
 static void updateOmega(double tau, double& omega1, double& omega2, double num,