Partial sovles with multiple rhs

Author: filikat

highs/ipm/hipo/factorhighs/FactorHighs.cpp

@@ -38,9 +38,13 @@ Int FHsolver::factorise(const Symbolic& S, Int n, Int nz, const Int* rows,
 }
 
 Int FHsolver::solve(double* x, Int k) const { return N_.solve(x, k); }
-Int FHsolver::forwardSolve(double* x) const { return N_.forwardSolve(x); }
-Int FHsolver::diagSolve(double* x) const { return N_.diagSolve(x); }
-Int FHsolver::backwardSolve(double* x) const { return N_.backwardSolve(x); }
+Int FHsolver::forwardSolve(double* x, Int k) const {
+  return N_.forwardSolve(x, k);
+}
+Int FHsolver::diagSolve(double* x, Int k) const { return N_.diagSolve(x, k); }
+Int FHsolver::backwardSolve(double* x, Int k) const {
+  return N_.backwardSolve(x, k);
+}
 
 void FHsolver::getRegularisation(double* reg) { N_.getReg(reg); }
 

highs/ipm/hipo/factorhighs/FactorHighs.h

@@ -103,9 +103,9 @@ class FHsolver {
   Int solve(double* x, Int k = 1) const;
 
   // Perform partial solves with L, D, L^T (including permutation)
-  Int forwardSolve(double* x) const;
-  Int diagSolve(double* x) const;
-  Int backwardSolve(double* x) const;
+  Int forwardSolve(double* x, Int k = 1) const;
+  Int diagSolve(double* x, Int k = 1) const;
+  Int backwardSolve(double* x, Int k = 1) const;
 
   // If multiple factorisation are performed, call newIter() before each
   // factorisation. This is used only to collect data for debugging, if

highs/ipm/hipo/factorhighs/FactorHighs_c_api.cpp

@@ -35,16 +35,16 @@ HighsInt FactorHighs_solve(void* FH, double* x, HighsInt k) {
   return ((hipo::FHsolver*)FH)->solve(x, k);
 }
 
-HighsInt FactorHighs_forwardSolve(void* FH, double* x) {
-  return ((hipo::FHsolver*)FH)->forwardSolve(x);
+HighsInt FactorHighs_forwardSolve(void* FH, double* x, HighsInt k) {
+  return ((hipo::FHsolver*)FH)->forwardSolve(x, k);
 }
 
-HighsInt FactorHighs_diagSolve(void* FH, double* x) {
-  return ((hipo::FHsolver*)FH)->diagSolve(x);
+HighsInt FactorHighs_diagSolve(void* FH, double* x, HighsInt k) {
+  return ((hipo::FHsolver*)FH)->diagSolve(x, k);
 }
 
-HighsInt FactorHighs_backwardSolve(void* FH, double* x) {
-  return ((hipo::FHsolver*)FH)->backwardSolve(x);
+HighsInt FactorHighs_backwardSolve(void* FH, double* x, HighsInt k) {
+  return ((hipo::FHsolver*)FH)->backwardSolve(x, k);
 }
 
 void FactorHighs_setRegularisation(void* FH, double reg_p, double reg_d) {

highs/ipm/hipo/factorhighs/FactorHighs_c_api.h

@@ -31,9 +31,9 @@ HighsInt FactorHighs_factorise(void* FH, const void* S, HighsInt n, HighsInt nz,
 HighsInt FactorHighs_solve(void* FH, double* x, HighsInt k);
 
 // Partial solves
-HighsInt FactorHighs_forwardSolve(void* FH, double* x);
-HighsInt FactorHighs_diagSolve(void* FH, double* x);
-HighsInt FactorHighs_backwardSolve(void* FH, double* x);
+HighsInt FactorHighs_forwardSolve(void* FH, double* x, HighsInt k);
+HighsInt FactorHighs_diagSolve(void* FH, double* x, HighsInt k);
+HighsInt FactorHighs_backwardSolve(void* FH, double* x, HighsInt k);
 
 // Reordering
 HighsInt FactorHighs_reorderMetis(void* FH, HighsInt n, HighsInt nz,

highs/ipm/hipo/factorhighs/Numeric.cpp

@@ -44,24 +44,6 @@ Int Numeric::solve(double* x) const {
   return kRetOk;
 }
 
-Int Numeric::solve(double* x, Int k) const {
-  if (k == 1)
-    return solve(x);
-  else {
-    highs::parallel::TaskGroup tg;
-    const Int n = S_->size();
-    std::atomic<bool> fail{false};
-    for (Int i = 0; i < k; ++i) {
-      tg.spawn([=, &fail]() {
-        Int status = solve(&x[i * n]);
-        if (status) fail.store(true, std::memory_order_relaxed);
-      });
-    }
-    tg.taskWait();
-    return fail;
-  }
-}
-
 Int Numeric::forwardSolve(double* x) const {
   if (!sn_columns_ || !S_ || !data_ || !options_) return kRetInvalidPointer;
   HybridSolveHandler SH(*S_, *sn_columns_, swaps_, pivot_2x2_, *data_,
@@ -86,6 +68,32 @@ Int Numeric::backwardSolve(double* x) const {
   return kRetOk;
 }
 
+#define SOLVE_MULTIPLE(f)                                        \
+  if (k == 1)                                                    \
+    return f(x);                                                 \
+  else {                                                         \
+    highs::parallel::TaskGroup tg;                               \
+    const Int n = S_->size();                                    \
+    std::atomic<bool> fail{false};                               \
+    for (Int i = 0; i < k; ++i) {                                \
+      tg.spawn([=, &fail]() {                                    \
+        Int status = f(&x[i * n]);                               \
+        if (status) fail.store(true, std::memory_order_relaxed); \
+      });                                                        \
+    }                                                            \
+    tg.taskWait();                                               \
+    return fail;                                                 \
+  }
+
+Int Numeric::solve(double* x, Int k) const { SOLVE_MULTIPLE(solve); }
+Int Numeric::forwardSolve(double* x, Int k) const {
+  SOLVE_MULTIPLE(forwardSolve);
+}
+Int Numeric::diagSolve(double* x, Int k) const { SOLVE_MULTIPLE(diagSolve); }
+Int Numeric::backwardSolve(double* x, Int k) const {
+  SOLVE_MULTIPLE(backwardSolve);
+}
+
 void Numeric::getReg(double* reg) {
   std::memcpy(reg, total_reg_.data(), total_reg_.size() * sizeof(double));
 }

highs/ipm/hipo/factorhighs/Numeric.h

@@ -46,6 +46,10 @@ class Numeric {
   Int diagSolve(double* x) const;
   Int backwardSolve(double* x) const;
 
+  Int forwardSolve(double* x, Int k) const;
+  Int diagSolve(double* x, Int k) const;
+  Int backwardSolve(double* x, Int k) const;
+
   void getReg(double* reg);
   void inertia(Int& pos, Int& neg, Int& zero, double tol) const;
 };