Adding option to pass factorized quadratic matrix

include/jrl-qp/SolverOptions.h

@@ -16,6 +16,10 @@ struct JRLQP_DLLAPI SolverOptions
   int maxIter_ = 500;
   double bigBnd_ = 1e100;
   bool warmStart_ = false;
+  /** If true, the input matrix G must contain the lower triangular matrix L
+   * such that G = L L^T. The upper part of the matrix is ignored.
+   */
+  bool factorizedG_ = false;
   std::uint32_t logFlags_ = 0;
   std::ostream * logStream_ = &defaultStream_;
 
@@ -76,6 +80,16 @@ struct JRLQP_DLLAPI SolverOptions
     return *this;
   }
 
+  bool factorizedG() const
+  {
+    return factorizedG_;
+  }
+  SolverOptions & factorizedG(bool fact)
+  {
+    factorizedG_ = fact;
+    return *this;
+  }
+
   std::ostream & logStream() const
   {
     return *logStream_;

src/GoldfarbIdnaniSolver.cpp

@@ -55,7 +55,7 @@ TerminationStatus GoldfarbIdnaniSolver::solve(MatrixRef G,
 
 internal::InitTermination GoldfarbIdnaniSolver::init_()
 {
-  auto ret = Eigen::internal::llt_inplace<double, Eigen::Lower>::blocked(pb_.G);
+  auto ret = (!options_.factorizedG_) ? Eigen::internal::llt_inplace<double, Eigen::Lower>::blocked(pb_.G) : -1;
   auto L = pb_.G.template triangularView<Eigen::Lower>();
 
   if(ret >= 0) return TerminationStatus::NON_POS_HESSIAN;

tests/CMakeLists.txt

@@ -57,6 +57,7 @@ add_dependencies(BlockGISolverTest extra-multiik-archive)
 addUnitTest(BoxAndSingleConstraintSolverTest)
 addUnitTest(GoldfarbIdnaniSolverTest QPSReader.cpp)
 addUnitTest(InternalTest)
+addUnitTest(OptionsTest)
 addUnitTest(RandomProblemsTest)
 addUnitTest(StructuredTest)
 addUnitTest(triBlockDiagLLTTest)

tests/OptionsTest.cpp

@@ -0,0 +1,37 @@
+/* Copyright 2024 CNRS-AIST JRL, Inria */
+
+#define DOCTEST_CONFIG_IMPLEMENT_WITH_MAIN
+#define DOCTEST_CONFIG_SUPER_FAST_ASSERTS
+#include "doctest/doctest.h"
+
+#include <jrl-qp/GoldfarbIdnaniSolver.h>
+#include <jrl-qp/test/randomProblems.h>
+
+using namespace jrl::qp::test;
+
+TEST_CASE("Test FactorizedG")
+{
+  jrl::qp::GoldfarbIdnaniSolver solver(7, 11, false);
+  jrl::qp::SolverOptions options;
+
+  for(int i = 0; i < 10; ++i)
+  {
+    auto pb = QPProblem(randomProblem(ProblemCharacteristics(7, 7, 3, 8)));
+    pb.C.transposeInPlace();
+
+    auto llt = pb.G.llt();
+    Eigen::MatrixXd L = llt.matrixL();
+
+    options.factorizedG(true);
+    solver.options(options);
+    solver.solve(L, pb.a, pb.C, pb.l, pb.u, pb.xl, pb.xu);
+    Eigen::VectorXd x1 = solver.solution();
+
+    options.factorizedG(false);
+    solver.options(options);
+    solver.solve(pb.G, pb.a, pb.C, pb.l, pb.u, pb.xl, pb.xu);
+    Eigen::VectorXd x2 = solver.solution();
+
+    FAST_CHECK_UNARY(x1.isApprox(x2));
+  }
+}