Adding Sparse matrix with prefactor related changes

Solver_interface/include/CGAL/Eigen_matrix.h

@@ -38,13 +38,13 @@ namespace CGAL {
   /// @heading Parameters:
   /// @param T Number type.
 
-template<class T>
+  template<class T, int Options = Eigen::RowMajor>
   struct Eigen_sparse_matrix
   {
     // Public types
   public:
 
-  typedef Eigen::SparseMatrix<T> EigenType;
+    typedef Eigen::SparseMatrix<T, Options> EigenType;
     typedef T NT;
 
     // Public operations
@@ -53,7 +53,7 @@ public:
     /// Create a square matrix initialized with zeros.
     Eigen_sparse_matrix(int  dim,                   ///< Matrix dimension.
       bool is_symmetric = false)  ///< Symmetric/hermitian?
-    : m_is_already_built(false), m_matrix(dim,dim)
+      : m_is_uptodate(false), m_matrix(dim,dim)
     {
       CGAL_precondition(dim > 0);
 
@@ -68,7 +68,7 @@ public:
     Eigen_sparse_matrix(int  rows,                 ///< Number of rows.
       int  columns,              ///< Number of columns.
       bool is_symmetric = false) ///< Symmetric/hermitian?
-    : m_is_already_built(false), m_matrix(rows,columns)
+      : m_is_uptodate(false), m_matrix(rows,columns)
     {
       CGAL_precondition(rows > 0);
       CGAL_precondition(columns > 0);
@@ -103,7 +103,7 @@ public:
     /// @commentheading Preconditions:
     /// - 0 <= i < row_dimension().
     /// - 0 <= j < column_dimension().
-  void set_coef(int i, int j, T  val, bool new_coef = false)
+    void set_coef(int i, int j, T  val, bool /* new_coef */ = false)
     {
       CGAL_precondition(i < row_dimension());
       CGAL_precondition(j < column_dimension());
@@ -111,18 +111,8 @@ public:
       if (m_is_symmetric && (j > i))
         return;
 
-    if (m_is_already_built)
-      m_matrix.coeffRef(i,j)=val;
-    else
-    {
-      if ( new_coef == false )
-      {
-        assemble_matrix();
-        m_matrix.coeffRef(i,j)=val;
-      }
-      else
       m_triplets.push_back(Triplet(i,j,val));
-    }
+      m_is_uptodate = false;
     }
 
     /// Write access to a matrix coefficient: a_ij <- a_ij+val.
@@ -142,23 +132,19 @@ public:
       if (m_is_symmetric && (j > i))
         return;
 
-    if (m_is_already_built)
-      m_matrix.coeffRef(i,j)+=val;
-    else
       m_triplets.push_back(Triplet(i,j,val));
+      m_is_uptodate = false;
     }  
 
-  void assemble_matrix() const
-  {
-    m_matrix.setFromTriplets(m_triplets.begin(), m_triplets.end());
-    m_is_already_built = true;
-    m_triplets.clear(); //the matrix is built and will not be rebuilt
-  }
+
 
     const EigenType& eigen_object() const
     {
-    if(!m_is_already_built) assemble_matrix();
-
+      if(!m_is_uptodate)
+      {
+        m_matrix.setFromTriplets(m_triplets.begin(), m_triplets.end());
+        m_is_uptodate = true;
+      }
       // turns the matrix into compressed mode:
       //  -> release some memory
       //  -> required for some external solvers
@@ -176,7 +162,7 @@ private:
     // Fields
   private:
 
-  mutable bool m_is_already_built;
+    mutable bool m_is_uptodate;
     typedef Eigen::Triplet<T,int> Triplet;
     mutable std::vector<Triplet> m_triplets;
 

Solver_interface/include/CGAL/Eigen_solver_traits.h

@@ -22,6 +22,7 @@
 #include <CGAL/basic.h> // include basic.h before testing #defines
 
 #include <Eigen/Sparse>
+#include <Eigen/SparseLU>
 #include <CGAL/Eigen_matrix.h>
 #include <CGAL/Eigen_vector.h>
 #include <boost/shared_ptr.hpp>
@@ -44,6 +45,11 @@ namespace internal {
     struct Get_eigen_matrix< ::Eigen::SimplicialCholesky<EigenMatrix>,FT>{
       typedef Eigen_sparse_symmetric_matrix<FT> type;
     };
+
+    template <class FT, class EigenMatrix, class EigenOrdering>
+    struct Get_eigen_matrix< ::Eigen::SparseLU<EigenMatrix, EigenOrdering >, FT> {
+      typedef Eigen_sparse_matrix<FT, ::Eigen::ColMajor> type;
+    };
   } //internal 
 
   /// The class Eigen_solver_traits
@@ -68,7 +74,7 @@ public:
     // Public operations
   public:
 
-   Eigen_solver_traits(): m_solver_sptr(new EigenSolverT)
+    Eigen_solver_traits():m_mat(NULL), m_solver_sptr(new EigenSolverT)
     {
     }
 
@@ -93,7 +99,24 @@ public:
 
       return m_solver_sptr->info() == Eigen::Success;
     }
+
+    bool pre_factor (const Matrix& A, NT& D)
+    {
+      D = 1;
+
+      m_mat = &A.eigen_object();
+      solver().compute(*m_mat);
+      return solver().info() == Eigen::Success;
+    }
+
+    bool linear_solver(const Vector& B, Vector& X)
+    {
+      CGAL_precondition(m_mat!=NULL); //pre_factor should have been called first
+      X = solver().solve(B);
+      return solver().info() == Eigen::Success;
+    }
   protected:
+    const typename Matrix::EigenType* m_mat;
     boost::shared_ptr<EigenSolverT> m_solver_sptr;
 
   };