Add Accelerate_sparse_matrix

Solver_interface/examples/Solver_interface/CMakeLists.txt

@@ -6,6 +6,8 @@ project(Solver_interface_Examples)
 
 find_package(CGAL REQUIRED)
 
+create_single_source_cgal_program("accelerate.cpp")
+
 # Use Eigen
 find_package(Eigen3 3.1.0 QUIET) #(requires 3.1.0 or greater)
 include(CGAL_Eigen3_support)

Solver_interface/examples/Solver_interface/accelerate.cpp

@@ -0,0 +1,34 @@
+#include <CGAL/Accelerate_sparse_matrix.h>
+
+int main()
+{
+  {
+    CGAL::Accelerate_sparse_matrix<double> sm(3);
+
+    sm.set_coef(0,0,1);
+    sm.set_coef(0,1,2);
+    sm.set_coef(0,2,3);
+    sm.set_coef(1,1,5);
+    sm.set_coef(2,2,9);
+
+    sm.assemble_matrix();
+
+    assert(sm.get_coef(1,0) == 0);
+  }
+  std::cout << std::endl;
+  {
+    CGAL::Accelerate_sparse_matrix<double> sm(3, true);// symmetric only set lower left
+
+    sm.set_coef(0,0,1);
+    sm.set_coef(1,0,2);
+    sm.set_coef(2,0,3);
+    sm.set_coef(1,1,5);
+    sm.set_coef(2,2,9);
+
+    sm.assemble_matrix();
+
+    assert(sm.get_coef(1,0) == sm.get_coef(0, 1));
+  }
+
+  return 0;
+}

Solver_interface/include/CGAL/Accelerate_sparse_matrix.h

@@ -0,0 +1,243 @@
+// Copyright (c) 2024  GeometryFactory SARL (France), All rights reserved.
+//
+// This file is part of CGAL (www.cgal.org)
+//
+// $URL$
+// $Id$
+// SPDX-License-Identifier: LGPL-3.0-or-later OR LicenseRef-Commercial
+//
+// Author(s)     : Andreas Fabri
+
+#ifndef CGAL_SOLVER_INTERFACE_ACCELERATE_SPARSE_MATRIX_H
+#define CGAL_SOLVER_INTERFACE_ACCELERATE_SPARSE_MATRIX_H
+
+#include <CGAL/basic.h> // include basic.h before testing #defines
+
+#include <map>
+
+namespace CGAL {
+
+/*!
+\ingroup PkgSolverInterfaceLS
+
+The class `Accelerate_sparse_matrix` is a wrapper around the `SparseMatrix_Double` matrix type
+<a href=https://developer.apple.com/documentation/accelerate/sparsematrix_double">`Accelerate::SparseMatrix`_Double</a>
+that represents general matrices, be they symmetric or not.
+
+\cgalModels{SparseLinearAlgebraTraits_d::Matrix}
+
+\tparam T Number type.
+
+\sa `CGAL::Accelerate_vector<T>`
+\sa `CGAL::Accelerate_matrix<T>`
+\sa `CGAL::Accelerate_sparse_symmetric_matrix<T>`
+*/
+template<class T>
+struct Accelerate_sparse_matrix
+{
+  // Public types
+public:
+  /// \name Types
+  /// @{
+
+  /// The internal matrix type from \ref thirdpartyAccelerate "Accelerate".
+  typedef int AccelerateType;
+
+  typedef T                           NT;
+  /// @}
+
+#if 0
+  Accelerate_sparse_matrix(const AccelerateType& et)
+    : m_is_already_built(true), m_has_been_changed(false), m_matrix(et), m_is_symmetric(false)
+  {}
+#endif
+
+  // Public operations
+public:
+  Accelerate_sparse_matrix()
+    : m_is_already_built(false), m_has_been_changed(false)
+  { }
+
+  /// Create a square matrix initialized with zeros.
+  Accelerate_sparse_matrix(int dim,            ///< Matrix dimension.
+                           bool is_symmetric = false)  ///< Symmetric/hermitian?
+    : m_rows(dim), m_columns(dim), m_is_already_built(false), m_has_been_changed(false),  m_is_symmetric(is_symmetric)
+  {
+    CGAL_precondition(dim > 0);
+  }
+
+
+  /// Create a rectangular matrix initialized with zeros.
+  ///
+  /// \pre rows == columns if `is_symmetric` is true.
+  Accelerate_sparse_matrix(int rows,          ///< Number of rows.
+                           int columns,       ///< Number of columns.
+                           bool is_symmetric = false) ///< Symmetric/hermitian?
+    : m_rows(rows), m_columns(columns), m_is_already_built(false), m_has_been_changed(false), m_is_symmetric(is_symmetric)
+  {
+    CGAL_precondition(rows > 0);
+    CGAL_precondition(columns > 0);
+    if(m_is_symmetric)
+    {
+      CGAL_precondition(rows == columns);
+    }
+
+    // reserve memory for a regular 3D grid
+  }
+
+  void swap(Accelerate_sparse_matrix& other)
+  {
+    std::swap(m_rows, other.m_rows);
+    std::swap(m_columns, other.m_columns);
+    std::swap(m_is_already_built, other.m_is_already_built);
+    std::swap(m_has_been_changed, other.m_has_been_changed);
+    std::swap(m_is_symmetric, other.m_is_symmetric);
+    //m_matrix.swap(other.m_matrix);
+    m_triplets.swap(other.m_triplets);
+  }
+
+  /// Delete this object and the wrapped matrix.
+  ~Accelerate_sparse_matrix() { }
+
+
+  /// Return the matrix number of rows
+  int row_dimension() const { return m_rows; }
+
+  /// Return the matrix number of columns
+  int column_dimension() const { return m_columns; }
+
+  /// Write access to a matrix coefficient: a_ij <- val.
+  ///
+  /// Users can optimize calls to this function by setting 'new_coef' to `true`
+  /// if the coefficient does not already exist in the matrix.
+  ///
+  /// \warning For symmetric matrices, `Accelerate_sparse_matrix` only stores the lower triangle
+  ///   and `set_coef()` does nothing if (i, j) belongs to the upper triangle.
+  ///
+  /// \pre 0 <= i < row_dimension().
+  /// \pre 0 <= j < column_dimension().
+  void set_coef(std::size_t i_, std::size_t j_, T val, bool new_coef = false)
+  {
+    int i = static_cast<int>(i_);
+    int j = static_cast<int>(j_);
+    CGAL_precondition(i < row_dimension());
+    CGAL_precondition(j < column_dimension());
+
+    if (m_is_symmetric && (j > i))
+      return;
+
+    m_triplets[std::make_pair(i,j)] = val;
+    if(m_is_already_built){
+      m_has_been_changed = true;
+    }
+  }
+
+  /// Write access to a matrix coefficient: a_ij <- a_ij + val.
+  ///
+  /// \warning For symmetric matrices, Accelerate_sparse_matrix only stores the lower triangle
+  /// `add_coef()` does nothing if (i, j) belongs to the upper triangle.
+  ///
+  /// \pre 0 <= i < row_dimension().
+  /// \pre 0 <= j < column_dimension().
+  void add_coef(std::size_t i_, std::size_t j_, T val)
+  {
+    int i = static_cast<int>(i_);
+    int j = static_cast<int>(j_);
+    if(m_is_symmetric && (j > i))
+      return;
+
+    m_triplets[std::make_pair(i,j)] += val;
+    if(m_is_already_built){
+      m_has_been_changed = true;
+    }
+  }
+
+  /// Read access to a matrix coefficient.
+  ///
+  /// \pre 0 <= i < row_dimension().
+  /// \pre 0 <= j < column_dimension().
+  NT get_coef (std::size_t i_, std::size_t j_) const
+  {
+    int i = static_cast<int>(i_);
+    int j = static_cast<int>(j_);
+    CGAL_precondition(i < row_dimension());
+    CGAL_precondition(j < column_dimension());
+
+    if(m_is_symmetric && j > i){
+      std::swap(i, j);
+    }
+    auto it = m_triplets.find(std::make_pair(i,j));
+    if(it != m_triplets.end()){
+      return it->second;
+    }
+    return 0;
+  }
+
+  /// \cond SKIP_IN_MANUAL
+  void assemble_matrix() const
+  {
+    std::vector<int> rows, columns;
+    std::vector<NT> values;
+    rows.reserve(m_triplets.size());
+    columns.reserve(m_triplets.size());
+    values.reserve(m_triplets.size());
+    for(const auto& triplet : m_triplets){
+      int i = triplet.first.first;
+      int j = triplet.first.second;
+      const NT& val = triplet.second;
+      rows.push_back(i);
+      columns.push_back(j);
+      values.push_back(val);
+      if(m_is_symmetric && (i != j)){
+        rows.push_back(j);
+        columns.push_back(i);
+        values.push_back(val);
+      }
+    }
+
+    for(int i = 0; i < rows.size(); ++i){
+      std::cout << rows[i] << "\t" << columns[i] << "\t" << values[i] << std::endl;
+    }
+    // todo create the swift sparse matrix from the vectors
+    m_is_already_built = true;
+    m_has_been_changed = false;
+  }
+  /// \endcond
+
+public:
+  /// \cond SKIP_IN_MANUAL
+  friend Accelerate_sparse_matrix
+  operator*(const T& c, const Accelerate_sparse_matrix& M)
+  {
+    std::cout << "todo: implement operator*(scalar, matrix)" << std::endl;
+    assert(false);
+    return M;
+  }
+
+  friend Accelerate_sparse_matrix
+  operator+(const Accelerate_sparse_matrix& M0, const Accelerate_sparse_matrix& M1)
+  {
+    std::cout << "todo: implement operator+( matrix, matrix)" << std::endl;
+    assert(false);
+    return M0;
+  }
+  /// \endcond
+
+  // Fields
+private:
+  mutable bool m_is_already_built;
+  mutable bool m_has_been_changed;
+  bool m_is_symmetric;
+  int m_rows, m_columns;
+  mutable std::map<std::pair<int,int>,NT> m_triplets;
+
+  // mutable AccelerateType m_matrix;
+
+}; // Accelerate_sparse_matrix
+
+
+
+
+} //namespace CGAL
+
+#endif // CGAL_SOLVER_INTERFACE_ACCELERATE_SPARSE_MATRIX_H