Added class DenseMatrix that will be used LBFGSHessian

Author: cvanaret

CMakeLists.txt

@@ -78,6 +78,7 @@ file(GLOB TESTS_UNO_SOURCE_FILES
    unotest/unit_tests/ConcatenationTests.cpp
    unotest/unit_tests/COOSparseStorageTests.cpp
    unotest/unit_tests/CSCSparseStorageTests.cpp
+   unotest/unit_tests/DenseMatrixTests.cpp
    unotest/unit_tests/MatrixVectorProductTests.cpp
    unotest/unit_tests/RangeTests.cpp
    unotest/unit_tests/ScalarMultipleTests.cpp

uno/ingredients/hessian_models/HessianModelFactory.cpp

@@ -6,9 +6,12 @@
 #include "HessianModel.hpp"
 #include "ConvexifiedHessian.hpp"
 #include "ExactHessian.hpp"
-#include "LBFGSHessian.hpp"
 #include "ZeroHessian.hpp"
 #include "ingredients/subproblem_solvers/DirectSymmetricIndefiniteLinearSolver.hpp"
+#ifdef HAS_LAPACK
+#include "LBFGSHessian.hpp"
+#include "options/Options.hpp"
+#endif
 
 namespace uno {
    std::unique_ptr<HessianModel> HessianModelFactory::create(const std::string& hessian_model, size_t dimension, size_t number_nonzeros,
@@ -22,8 +25,8 @@ namespace uno {
          }
       }
 #ifdef HAS_LAPACK
-      else if (true || hessian_model == "L-BFGS") {
-         return std::make_unique<LBFGSHessian>(dimension);
+      else if (hessian_model == "L-BFGS") {
+         return std::make_unique<LBFGSHessian>(dimension, options.get_unsigned_int("quasi_newton_memory_size"));
       }
 #endif
       else if (hessian_model == "zero") {

uno/ingredients/hessian_models/LBFGSHessian.cpp

@@ -13,7 +13,13 @@ extern "C" {
 #endif
 
 namespace uno {
-   LBFGSHessian::LBFGSHessian(size_t dimension): HessianModel(), dimension(dimension) {
+   LBFGSHessian::LBFGSHessian(size_t dimension, size_t memory_size):
+         HessianModel(),
+         dimension(dimension),
+         memory_size(memory_size),
+         S_matrix(dimension, memory_size),
+         Y_matrix(dimension, memory_size),
+         M_matrix(memory_size, memory_size) {
    }
 
    void LBFGSHessian::initialize_statistics(Statistics& /*statistics*/, const Options& /*options*/) const {

uno/ingredients/hessian_models/LBFGSHessian.hpp

@@ -1,15 +1,16 @@
 // Copyright (c) 2025 Charlie Vanaret
 // Licensed under the MIT license. See LICENSE file in the project directory for details.
 
-#include "HessianModel.hpp"
-
 #ifndef UNO_LBFGSHESSIAN_H
 #define UNO_LBFGSHESSIAN_H
 
+#include "HessianModel.hpp"
+#include "linear_algebra/DenseMatrix.hpp"
+
 namespace uno {
    class LBFGSHessian : public HessianModel {
    public:
-      explicit LBFGSHessian(size_t dimension);
+      LBFGSHessian(size_t dimension, size_t memory_size);
       ~LBFGSHessian() override = default;
 
       void initialize_statistics(Statistics& statistics, const Options& options) const override;
@@ -21,6 +22,10 @@ namespace uno {
 
    protected:
       const size_t dimension;
+      const size_t memory_size;
+      DenseMatrix<double> S_matrix;
+      DenseMatrix<double> Y_matrix;
+      DenseMatrix<double> M_matrix;
    };
 } // namespace
 

uno/linear_algebra/CSCSparseStorage.hpp

@@ -8,7 +8,6 @@
 #include <vector>
 #include "SparseStorage.hpp"
 #include "linear_algebra/Vector.hpp"
-#include "tools/Infinity.hpp"
 #include "symbolic/VectorView.hpp"
 
 namespace uno {

uno/linear_algebra/DenseMatrix.hpp

@@ -0,0 +1,42 @@
+// Copyright (c) 2025 Charlie Vanaret
+// Licensed under the MIT license. See LICENSE file in the project directory for details.
+
+#ifndef UNO_DENSEMATRIX_H
+#define UNO_DENSEMATRIX_H
+
+#include <cstddef>
+#include <vector>
+
+namespace uno {
+   template <typename ElementType>
+   class DenseMatrix {
+   public:
+      DenseMatrix(size_t number_rows, size_t number_columns);
+      ~DenseMatrix() = default;
+
+      ElementType& get(size_t row_index, size_t column_index);
+      ElementType* data() const;
+
+   protected:
+      const size_t number_rows, number_columns;
+      std::vector<ElementType> matrix; // column-major ordering
+   };
+
+   template <typename ElementType>
+   DenseMatrix<ElementType>::DenseMatrix(size_t number_rows, size_t number_columns):
+         number_rows(number_rows), number_columns(number_columns),
+         matrix(number_rows * number_columns, ElementType(0)) {
+   }
+
+   template <typename ElementType>
+   ElementType& DenseMatrix<ElementType>::get(size_t row_index, size_t column_index) {
+      return this->matrix[column_index * this->number_rows + row_index];
+   }
+
+   template <typename ElementType>
+   ElementType* DenseMatrix<ElementType>::data() const {
+      return this->matrix.data();
+   }
+} // namespace
+
+#endif // UNO_DENSEMATRIX_H

unotest/unit_tests/DenseMatrixTests.cpp

@@ -0,0 +1,12 @@
+// Copyright (c) 2025 Charlie Vanaret
+// Licensed under the MIT license. See LICENSE file in the project directory for details.
+
+#include <gtest/gtest.h>
+#include "linear_algebra/DenseMatrix.hpp"
+
+using namespace uno;
+
+TEST(DenseMatrix, Access) {
+   DenseMatrix<double> M(5, 2);
+   ASSERT_EQ(M.get(1, 1), 0.);
+}