diff --git a/Solver_interface/include/CGAL/Eigen_matrix.h b/Solver_interface/include/CGAL/Eigen_matrix.h index 906229661fc..bf9162fcbdc 100644 --- a/Solver_interface/include/CGAL/Eigen_matrix.h +++ b/Solver_interface/include/CGAL/Eigen_matrix.h @@ -1,252 +1,238 @@ -// Copyright (c) 2012 INRIA Bordeaux Sud-Ouest (France), All rights reserved. -// -// This file is part of CGAL (www.cgal.org); you can redistribute it and/or -// modify it under the terms of the GNU Lesser General Public License as -// published by the Free Software Foundation; either version 3 of the License, -// or (at your option) any later version. -// -// Licensees holding a valid commercial license may use this file in -// accordance with the commercial license agreement provided with the software. -// -// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE -// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. -// -// $URL$ -// $Id$ -// -// Author(s) : Gael Guennebaud - -#ifndef CGAL_EIGEN_MATRIX_H -#define CGAL_EIGEN_MATRIX_H - -#include // include basic.h before testing #defines - -#define EIGEN_YES_I_KNOW_SPARSE_MODULE_IS_NOT_STABLE_YET -#include - -namespace CGAL { - - -/// The class Eigen_sparse_matrix -/// is a C++ wrapper around Eigen' matrix type SparseMatrix<>. -/// -/// This kind of matrix can be either symmetric or not. Symmetric -/// matrices store only the lower triangle. -/// -/// @heading Is Model for the Concepts: Model of the SparseLinearAlgebraTraits_d::Matrix concept. -/// -/// @heading Parameters: -/// @param T Number type. - -template -struct Eigen_sparse_matrix -{ -// Public types -public: - - typedef Eigen::SparseMatrix EigenType; - typedef T NT; - -// Public operations -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) - { - CGAL_precondition(dim > 0); - - m_is_symmetric = is_symmetric; - // reserve memory for a regular 3D grid - m_triplets.reserve(dim); - } - - /// Create a rectangular matrix initialized with zeros. - /// - /// @commentheading Precondition: rows == columns if is_symmetric is true. - 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) - { - CGAL_precondition(rows > 0); - CGAL_precondition(columns > 0); - if (m_is_symmetric) { - CGAL_precondition(rows == columns); - } - - m_is_symmetric = is_symmetric; - // reserve memory for a regular 3D grid - m_triplets.reserve(rows); - } - - /// Delete this object and the wrapped TAUCS matrix. - ~Eigen_sparse_matrix() - { - } - - /// Return the matrix number of rows - int row_dimension() const { return m_matrix.rows(); } - /// Return the matrix number of columns - int column_dimension() const { return m_matrix.cols(); } - - - /// Write access to a matrix coefficient: a_ij <- val. - /// - /// Optimizations: - /// - For symmetric matrices, Eigen_sparse_matrix stores only the lower triangle - /// set_coef() does nothing if (i, j) belongs to the upper triangle. - /// - Caller can optimize this call by setting 'new_coef' to true - /// if the coefficient does not already exist in the matrix. - /// - /// @commentheading Preconditions: - /// - 0 <= i < row_dimension(). - /// - 0 <= j < column_dimension(). - void set_coef(int i, int j, T val, bool new_coef = false) - { - CGAL_precondition(i < row_dimension()); - CGAL_precondition(j < column_dimension()); - - 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)); - } - } - - /// Write access to a matrix coefficient: a_ij <- a_ij+val. - /// - /// Optimizations: - /// - For symmetric matrices, Eigen_sparse_matrix stores only the lower triangle - /// add_coef() does nothing if (i, j) belongs to the upper triangle. - /// - /// @commentheading Preconditions: - /// - 0 <= i < row_dimension(). - /// - 0 <= j < column_dimension(). - void add_coef(int i, int j, T val) - { - CGAL_precondition(i < row_dimension()); - CGAL_precondition(j < column_dimension()); - - 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)); - } - - 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(); - - // turns the matrix into compressed mode: - // -> release some memory - // -> required for some external solvers - m_matrix.makeCompressed(); - return m_matrix; - } - -private: - - - /// Eigen_sparse_matrix cannot be copied (yet) - Eigen_sparse_matrix(const Eigen_sparse_matrix& rhs); - Eigen_sparse_matrix& operator=(const Eigen_sparse_matrix& rhs); - -// Fields -private: - - mutable bool m_is_already_built; - typedef Eigen::Triplet Triplet; - mutable std::vector m_triplets; - - mutable EigenType m_matrix; - - // Symmetric/hermitian? - bool m_is_symmetric; - -}; // Eigen_sparse_matrix - - - -/// The class Eigen_sparse_symmetric_matrix is a C++ wrapper -/// around a Eigen sparse matrix (type Eigen::SparseMatrix). -/// -/// Symmetric matrices store only the lower triangle. -/// -/// @heading Is Model for the Concepts: Model of the SparseLinearAlgebraTraits_d::Matrix concept. -/// -/// @heading Parameters: -/// @param T Number type. - -template -struct Eigen_sparse_symmetric_matrix - : public Eigen_sparse_matrix -{ -// Public types - typedef T NT; - -// Public operations - - /// Create a square *symmetric* matrix initialized with zeros. - Eigen_sparse_symmetric_matrix(int dim) ///< Matrix dimension. - : Eigen_sparse_matrix(dim, true /* symmetric */) - { - } - - /// Create a square *symmetric* matrix initialized with zeros. - /// - /// @commentheading Precondition: rows == columns. - Eigen_sparse_symmetric_matrix(int rows, ///< Number of rows. - int columns) ///< Number of columns. - : Eigen_sparse_matrix(rows, columns, true /* symmetric */) - { - } -}; - -template -struct Eigen_matrix : public ::Eigen::Matrix -{ - typedef ::Eigen::Matrix EigenType; - - Eigen_matrix( std::size_t n1, std::size_t n2):EigenType(n1,n2){} - - std::size_t number_of_rows () const {return this->rows();} - - std::size_t number_of_columns () const {return this->cols();} - - FT operator()( std::size_t i , std::size_t j ) const {return this->operator()(i,j);} - - void set( std::size_t i, std::size_t j,FT value){ - this->coeffRef(i,j)=value; - } - - const EigenType& eigen_object() const{ - return static_cast(*this); - } - -}; - -} //namespace CGAL - -#endif // CGAL_EIGEN_MATRIX_H +// Copyright (c) 2012 INRIA Bordeaux Sud-Ouest (France), All rights reserved. +// +// This file is part of CGAL (www.cgal.org); you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public License as +// published by the Free Software Foundation; either version 3 of the License, +// or (at your option) any later version. +// +// Licensees holding a valid commercial license may use this file in +// accordance with the commercial license agreement provided with the software. +// +// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE +// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. +// +// $URL$ +// $Id$ +// +// Author(s) : Gael Guennebaud + +#ifndef CGAL_EIGEN_MATRIX_H +#define CGAL_EIGEN_MATRIX_H + +#include // include basic.h before testing #defines + +#define EIGEN_YES_I_KNOW_SPARSE_MODULE_IS_NOT_STABLE_YET +#include + +namespace CGAL { + + + /// The class Eigen_sparse_matrix + /// is a C++ wrapper around Eigen' matrix type SparseMatrix<>. + /// + /// This kind of matrix can be either symmetric or not. Symmetric + /// matrices store only the lower triangle. + /// + /// @heading Is Model for the Concepts: Model of the SparseLinearAlgebraTraits_d::Matrix concept. + /// + /// @heading Parameters: + /// @param T Number type. + + template + struct Eigen_sparse_matrix + { + // Public types + public: + + typedef Eigen::SparseMatrix EigenType; + typedef T NT; + + // Public operations + public: + + /// Create a square matrix initialized with zeros. + Eigen_sparse_matrix(int dim, ///< Matrix dimension. + bool is_symmetric = false) ///< Symmetric/hermitian? + : m_is_uptodate(false), m_matrix(dim,dim) + { + CGAL_precondition(dim > 0); + + m_is_symmetric = is_symmetric; + // reserve memory for a regular 3D grid + m_triplets.reserve(dim); + } + + /// Create a rectangular matrix initialized with zeros. + /// + /// @commentheading Precondition: rows == columns if is_symmetric is true. + Eigen_sparse_matrix(int rows, ///< Number of rows. + int columns, ///< Number of columns. + bool is_symmetric = false) ///< Symmetric/hermitian? + : m_is_uptodate(false), m_matrix(rows,columns) + { + CGAL_precondition(rows > 0); + CGAL_precondition(columns > 0); + if (m_is_symmetric) { + CGAL_precondition(rows == columns); + } + + m_is_symmetric = is_symmetric; + // reserve memory for a regular 3D grid + m_triplets.reserve(rows); + } + + /// Delete this object and the wrapped TAUCS matrix. + ~Eigen_sparse_matrix() + { + } + + /// Return the matrix number of rows + int row_dimension() const { return m_matrix.rows(); } + /// Return the matrix number of columns + int column_dimension() const { return m_matrix.cols(); } + + + /// Write access to a matrix coefficient: a_ij <- val. + /// + /// Optimizations: + /// - For symmetric matrices, Eigen_sparse_matrix stores only the lower triangle + /// set_coef() does nothing if (i, j) belongs to the upper triangle. + /// - Caller can optimize this call by setting 'new_coef' to true + /// if the coefficient does not already exist in the matrix. + /// + /// @commentheading Preconditions: + /// - 0 <= i < row_dimension(). + /// - 0 <= j < column_dimension(). + void set_coef(int i, int j, T val, bool /* new_coef */ = false) + { + CGAL_precondition(i < row_dimension()); + CGAL_precondition(j < column_dimension()); + + if (m_is_symmetric && (j > i)) + return; + + m_triplets.push_back(Triplet(i,j,val)); + m_is_uptodate = false; + } + + /// Write access to a matrix coefficient: a_ij <- a_ij+val. + /// + /// Optimizations: + /// - For symmetric matrices, Eigen_sparse_matrix stores only the lower triangle + /// add_coef() does nothing if (i, j) belongs to the upper triangle. + /// + /// @commentheading Preconditions: + /// - 0 <= i < row_dimension(). + /// - 0 <= j < column_dimension(). + void add_coef(int i, int j, T val) + { + CGAL_precondition(i < row_dimension()); + CGAL_precondition(j < column_dimension()); + + if (m_is_symmetric && (j > i)) + return; + + m_triplets.push_back(Triplet(i,j,val)); + m_is_uptodate = false; + } + + + + const EigenType& eigen_object() const + { + 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 + m_matrix.makeCompressed(); + return m_matrix; + } + + private: + + + /// Eigen_sparse_matrix cannot be copied (yet) + Eigen_sparse_matrix(const Eigen_sparse_matrix& rhs); + Eigen_sparse_matrix& operator=(const Eigen_sparse_matrix& rhs); + + // Fields + private: + + mutable bool m_is_uptodate; + typedef Eigen::Triplet Triplet; + mutable std::vector m_triplets; + + mutable EigenType m_matrix; + + // Symmetric/hermitian? + bool m_is_symmetric; + + }; // Eigen_sparse_matrix + + + + /// The class Eigen_sparse_symmetric_matrix is a C++ wrapper + /// around a Eigen sparse matrix (type Eigen::SparseMatrix). + /// + /// Symmetric matrices store only the lower triangle. + /// + /// @heading Is Model for the Concepts: Model of the SparseLinearAlgebraTraits_d::Matrix concept. + /// + /// @heading Parameters: + /// @param T Number type. + + template + struct Eigen_sparse_symmetric_matrix + : public Eigen_sparse_matrix + { + // Public types + typedef T NT; + + // Public operations + + /// Create a square *symmetric* matrix initialized with zeros. + Eigen_sparse_symmetric_matrix(int dim) ///< Matrix dimension. + : Eigen_sparse_matrix(dim, true /* symmetric */) + { + } + + /// Create a square *symmetric* matrix initialized with zeros. + /// + /// @commentheading Precondition: rows == columns. + Eigen_sparse_symmetric_matrix(int rows, ///< Number of rows. + int columns) ///< Number of columns. + : Eigen_sparse_matrix(rows, columns, true /* symmetric */) + { + } + }; + + template + struct Eigen_matrix : public ::Eigen::Matrix + { + typedef ::Eigen::Matrix EigenType; + + Eigen_matrix( std::size_t n1, std::size_t n2):EigenType(n1,n2){} + + std::size_t number_of_rows () const {return this->rows();} + + std::size_t number_of_columns () const {return this->cols();} + + FT operator()( std::size_t i , std::size_t j ) const {return this->operator()(i,j);} + + void set( std::size_t i, std::size_t j,FT value){ + this->coeffRef(i,j)=value; + } + + const EigenType& eigen_object() const{ + return static_cast(*this); + } + + }; + +} //namespace CGAL + +#endif // CGAL_EIGEN_MATRIX_H diff --git a/Solver_interface/include/CGAL/Eigen_solver_traits.h b/Solver_interface/include/CGAL/Eigen_solver_traits.h index fa8d782d9bc..822b3aa455d 100644 --- a/Solver_interface/include/CGAL/Eigen_solver_traits.h +++ b/Solver_interface/include/CGAL/Eigen_solver_traits.h @@ -1,151 +1,174 @@ -// Copyright (c) 2012 INRIA Bordeaux Sud-Ouest (France), All rights reserved. -// -// This file is part of CGAL (www.cgal.org); you can redistribute it and/or -// modify it under the terms of the GNU Lesser General Public License as -// published by the Free Software Foundation; either version 3 of the License, -// or (at your option) any later version. -// -// Licensees holding a valid commercial license may use this file in -// accordance with the commercial license agreement provided with the software. -// -// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE -// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. -// -// $URL$ -// $Id$ -// -// Author(s) : Gael Guennebaud - -#ifndef CGAL_EIGEN_SOLVER_TRAITS_H -#define CGAL_EIGEN_SOLVER_TRAITS_H - -#include // include basic.h before testing #defines - -#include -#include -#include -#include - -namespace CGAL { - - -namespace internal { - template - struct Get_eigen_matrix{ - typedef Eigen_sparse_matrix type; - }; - - template - struct Get_eigen_matrix< ::Eigen::ConjugateGradient,FT>{ - typedef Eigen_sparse_symmetric_matrix type; - }; - - template - struct Get_eigen_matrix< ::Eigen::SimplicialCholesky,FT>{ - typedef Eigen_sparse_symmetric_matrix type; - }; -} //internal - -/// The class Eigen_solver_traits -/// is a generic traits class for solving asymmetric or symmetric positive definite (SPD) -/// sparse linear systems using one of the Eigen solvers. -/// The default solver is the iterative bi-congugate gradient stabilized solver -/// Eigen::BiCGSTAB for double. -/// -/// @heading Is Model for the Concepts: Model of the SparseLinearAlgebraTraits_d concept. - -template::EigenType> > -class Eigen_solver_traits -{ - typedef typename EigenSolverT::Scalar Scalar; -// Public types -public: - typedef Scalar NT; - typedef typename internal::Get_eigen_matrix::type Matrix; - typedef Eigen_vector Vector; - - -// Public operations -public: - - Eigen_solver_traits(): m_solver_sptr(new EigenSolverT) - { - } - - EigenSolverT& solver() { return *m_solver_sptr; } - - /// Solve the sparse linear system "A*X = B". - /// Return true on success. The solution is then (1/D) * X. - /// - /// @commentheading Preconditions: - /// - A.row_dimension() == B.dimension(). - /// - A.column_dimension() == X.dimension(). - bool linear_solver(const Matrix& A, const Vector& B, Vector& X, NT& D) - { - D = 1; // Eigen does not support homogeneous coordinates - - m_solver_sptr->compute(A.eigen_object()); - - if(m_solver_sptr->info() != Eigen::Success) - return false; - - X = m_solver_sptr->solve(B); - - return m_solver_sptr->info() == Eigen::Success; - } -protected: - boost::shared_ptr m_solver_sptr; - -}; - -//specilization of the solver for BiCGSTAB as for surface parameterization, the -//intializer should be a vector of one's (this was the case in 3.1-alpha but not in the official 3.1). -template<> -class Eigen_solver_traits< Eigen::BiCGSTAB::EigenType> > -{ - typedef Eigen::BiCGSTAB::EigenType> EigenSolverT; - typedef EigenSolverT::Scalar Scalar; -// Public types -public: - typedef Scalar NT; - typedef internal::Get_eigen_matrix::type Matrix; - typedef Eigen_vector Vector; - - -// Public operations -public: - - Eigen_solver_traits(): m_solver_sptr(new EigenSolverT) - { - } - - EigenSolverT& solver() { return *m_solver_sptr; } - - /// Solve the sparse linear system "A*X = B". - /// Return true on success. The solution is then (1/D) * X. - /// - /// @commentheading Preconditions: - /// - A.row_dimension() == B.dimension(). - /// - A.column_dimension() == X.dimension(). - bool linear_solver(const Matrix& A, const Vector& B, Vector& X, NT& D) - { - D = 1; // Eigen does not support homogeneous coordinates - - m_solver_sptr->compute(A.eigen_object()); - - if(m_solver_sptr->info() != Eigen::Success) - return false; - - X.setOnes(B.rows()); - X = m_solver_sptr->solveWithGuess(B,X); - - return m_solver_sptr->info() == Eigen::Success; - } -protected: - boost::shared_ptr m_solver_sptr; - -}; - -} //namespace CGAL - -#endif // CGAL_EIGEN_SOLVER_TRAITS_H +// Copyright (c) 2012 INRIA Bordeaux Sud-Ouest (France), All rights reserved. +// +// This file is part of CGAL (www.cgal.org); you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public License as +// published by the Free Software Foundation; either version 3 of the License, +// or (at your option) any later version. +// +// Licensees holding a valid commercial license may use this file in +// accordance with the commercial license agreement provided with the software. +// +// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE +// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. +// +// $URL$ +// $Id$ +// +// Author(s) : Gael Guennebaud + +#ifndef CGAL_EIGEN_SOLVER_TRAITS_H +#define CGAL_EIGEN_SOLVER_TRAITS_H + +#include // include basic.h before testing #defines + +#include +#include +#include +#include +#include + +namespace CGAL { + + + namespace internal { + template + struct Get_eigen_matrix{ + typedef Eigen_sparse_matrix type; + }; + + template + struct Get_eigen_matrix< ::Eigen::ConjugateGradient,FT>{ + typedef Eigen_sparse_symmetric_matrix type; + }; + + template + struct Get_eigen_matrix< ::Eigen::SimplicialCholesky,FT>{ + typedef Eigen_sparse_symmetric_matrix type; + }; + + template + struct Get_eigen_matrix< ::Eigen::SparseLU, FT> { + typedef Eigen_sparse_matrix type; + }; + } //internal + + /// The class Eigen_solver_traits + /// is a generic traits class for solving asymmetric or symmetric positive definite (SPD) + /// sparse linear systems using one of the Eigen solvers. + /// The default solver is the iterative bi-congugate gradient stabilized solver + /// Eigen::BiCGSTAB for double. + /// + /// @heading Is Model for the Concepts: Model of the SparseLinearAlgebraTraits_d concept. + + template::EigenType> > + class Eigen_solver_traits + { + typedef typename EigenSolverT::Scalar Scalar; + // Public types + public: + typedef Scalar NT; + typedef typename internal::Get_eigen_matrix::type Matrix; + typedef Eigen_vector Vector; + + + // Public operations + public: + + Eigen_solver_traits():m_mat(NULL), m_solver_sptr(new EigenSolverT) + { + } + + EigenSolverT& solver() { return *m_solver_sptr; } + + /// Solve the sparse linear system "A*X = B". + /// Return true on success. The solution is then (1/D) * X. + /// + /// @commentheading Preconditions: + /// - A.row_dimension() == B.dimension(). + /// - A.column_dimension() == X.dimension(). + bool linear_solver(const Matrix& A, const Vector& B, Vector& X, NT& D) + { + D = 1; // Eigen does not support homogeneous coordinates + + m_solver_sptr->compute(A.eigen_object()); + + if(m_solver_sptr->info() != Eigen::Success) + return false; + + X = m_solver_sptr->solve(B); + + 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 m_solver_sptr; + + }; + + //specilization of the solver for BiCGSTAB as for surface parameterization, the + //intializer should be a vector of one's (this was the case in 3.1-alpha but not in the official 3.1). + template<> + class Eigen_solver_traits< Eigen::BiCGSTAB::EigenType> > + { + typedef Eigen::BiCGSTAB::EigenType> EigenSolverT; + typedef EigenSolverT::Scalar Scalar; + // Public types + public: + typedef Scalar NT; + typedef internal::Get_eigen_matrix::type Matrix; + typedef Eigen_vector Vector; + + + // Public operations + public: + + Eigen_solver_traits(): m_solver_sptr(new EigenSolverT) + { + } + + EigenSolverT& solver() { return *m_solver_sptr; } + + /// Solve the sparse linear system "A*X = B". + /// Return true on success. The solution is then (1/D) * X. + /// + /// @commentheading Preconditions: + /// - A.row_dimension() == B.dimension(). + /// - A.column_dimension() == X.dimension(). + bool linear_solver(const Matrix& A, const Vector& B, Vector& X, NT& D) + { + D = 1; // Eigen does not support homogeneous coordinates + + m_solver_sptr->compute(A.eigen_object()); + + if(m_solver_sptr->info() != Eigen::Success) + return false; + + X.setOnes(B.rows()); + X = m_solver_sptr->solveWithGuess(B,X); + + return m_solver_sptr->info() == Eigen::Success; + } + protected: + boost::shared_ptr m_solver_sptr; + + }; + +} //namespace CGAL + +#endif // CGAL_EIGEN_SOLVER_TRAITS_H