diff --git a/.gitattributes b/.gitattributes index 018131673f6..e1f146d5aa9 100644 --- a/.gitattributes +++ b/.gitattributes @@ -3436,6 +3436,14 @@ Snap_rounding_2/doc_tex/Snap_rounding_2/sr1.pdf -text svneol=unset#application/p Snap_rounding_2/test/Snap_rounding_2/cgal_test eol=lf Snap_rounding_2/test/Snap_rounding_2/cgal_test_base -text Snap_rounding_2/test/Snap_rounding_2/cgal_test_with_cmake eol=lf +Solver_interface/include/CGAL/Eigen_matrix.h -text +Solver_interface/include/CGAL/Eigen_solver_traits.h -text +Solver_interface/include/CGAL/Eigen_svd.h -text +Solver_interface/include/CGAL/Eigen_vector.h -text +Solver_interface/package_info/Solver_interface/copyright -text +Solver_interface/package_info/Solver_interface/description.txt -text +Solver_interface/package_info/Solver_interface/license.txt -text +Solver_interface/package_info/Solver_interface/maintainer -text Spatial_searching/benchmark/Spatial_searching/Compare_ANN_STANN_CGAL.cpp -text Spatial_searching/demo/Spatial_searching/Qt3/help/index.html svneol=native#text/html Spatial_searching/doc_tex/Spatial_searching/Fig1.gif -text svneol=unset#image/gif diff --git a/Solver_interface/include/CGAL/Eigen_matrix.h b/Solver_interface/include/CGAL/Eigen_matrix.h new file mode 100644 index 00000000000..ac04da4a986 --- /dev/null +++ b/Solver_interface/include/CGAL/Eigen_matrix.h @@ -0,0 +1,230 @@ +// 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_matrix(dim,dim) + { + CGAL_precondition(dim > 0); + + m_is_symmetric = is_symmetric; + + // reserve memory for a regular 3D grid + m_matrix.reserve(Eigen::VectorXi::Constant(m_matrix.outerSize(), 27)); + } + + /// 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_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_matrix.reserve(Eigen::VectorXi::Constant(m_matrix.outerSize(), 27)); + } + + /// 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(new_coef) m_matrix.insert(i,j) = val; + else m_matrix.coeffRef(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; + + m_matrix.coeffRef(i,j) += val; + } + + + + const EigenType& eigen_object() const + { + // 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 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 new file mode 100644 index 00000000000..bb96dade078 --- /dev/null +++ b/Solver_interface/include/CGAL/Eigen_solver_traits.h @@ -0,0 +1,102 @@ +// 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 + +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() + { + } + + EigenSolverT& solver() { return m_solver; } + + /// 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.compute(A.eigen_object()); + + if(m_solver.info() != Eigen::Success) + return false; + + X = m_solver.solve(B); + + return m_solver.info() == Eigen::Success; + } +protected: + EigenSolverT m_solver; + +}; + +} //namespace CGAL + +#endif // CGAL_EIGEN_SOLVER_TRAITS_H diff --git a/Solver_interface/include/CGAL/Eigen_svd.h b/Solver_interface/include/CGAL/Eigen_svd.h new file mode 100644 index 00000000000..e35f8e8d7e0 --- /dev/null +++ b/Solver_interface/include/CGAL/Eigen_svd.h @@ -0,0 +1,46 @@ +// 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_SVD_H +#define CGAL_EIGEN_SVD_H + +#include +#include +#include + + + +namespace CGAL { + +class Eigen_svd{ +public: + typedef double FT; + typedef Eigen_vector Vector; + typedef Eigen_matrix Matrix; + //solve MX=B using SVD and return the condition number of M + //The solution is stored in B + static FT solve(const Matrix& M, Vector& B){ + Eigen::JacobiSVD jacobiSvd(M.eigen_object(),::Eigen::ComputeThinU | ::Eigen::ComputeThinV); + B.eigen_object()=jacobiSvd.solve(Vector::EigenType(B.eigen_object())); + return jacobiSvd.singularValues().array().abs().maxCoeff()/ + jacobiSvd.singularValues().array().abs().minCoeff(); + } +}; + +}//namespace CGAL +#endif // CGAL_EIGEN_SVD_H diff --git a/Solver_interface/include/CGAL/Eigen_vector.h b/Solver_interface/include/CGAL/Eigen_vector.h new file mode 100644 index 00000000000..dd5c1ad9b6d --- /dev/null +++ b/Solver_interface/include/CGAL/Eigen_vector.h @@ -0,0 +1,96 @@ +// 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_VECTOR_H +#define CGAL_EIGEN_VECTOR_H + +#include // include basic.h before testing #defines + +#include + +namespace CGAL { + +/// The class Eigen_vector +/// is a C++ wrapper around Eigen' dense vector type Matrix. +/// +/// @heading Is Model for the Concepts: Model of the SparseLinearAlgebraTraits_d::Vector concept. + +template +class Eigen_vector : public Eigen::Matrix +{ +// Public types +public: + + typedef T NT; + typedef Eigen::Matrix EigenType; + +// Public operations +public: + + Eigen_vector& operator=(const Eigen_vector& other){ + return static_cast(*this) = other.eigen_object(); + } + + Eigen_vector& operator=(const EigenType& other){ + return + static_cast&>( + static_cast(*this) = other + ); + } + + /// Create a vector initialized with zeros. + Eigen_vector(int dimension) + : EigenType(dimension) + { + this->setZero(); + } + + /// Copy constructor. + Eigen_vector(const Eigen_vector& toCopy) + : EigenType(toCopy) + { + } + + ~Eigen_vector() + { + } + + /// Return the vector's number of coefficients. + int dimension() const { + return this->size(); + } + + /// Get TAUCS vector wrapped by this object. + const EigenType& eigen_object() const { + return *this; + } + EigenType& eigen_object() { + return *this; + } + + void set(int i,NT value) { + this->operator[](i)=value; + } + + NT* vector() {return this->data();} + +}; + +} //namespace CGAL + +#endif // CGAL_EIGEN_VECTOR_H diff --git a/Solver_interface/package_info/Solver_interface/copyright b/Solver_interface/package_info/Solver_interface/copyright new file mode 100644 index 00000000000..a510ca9b8ec --- /dev/null +++ b/Solver_interface/package_info/Solver_interface/copyright @@ -0,0 +1 @@ +INRIA Bordeaux Sud-Ouest (France) diff --git a/Solver_interface/package_info/Solver_interface/description.txt b/Solver_interface/package_info/Solver_interface/description.txt new file mode 100644 index 00000000000..b934d65f214 --- /dev/null +++ b/Solver_interface/package_info/Solver_interface/description.txt @@ -0,0 +1,2 @@ +Solver_interface: contains classes to interface third party linear algebra solvers +to CGAL (such as Eigen for example) diff --git a/Solver_interface/package_info/Solver_interface/license.txt b/Solver_interface/package_info/Solver_interface/license.txt new file mode 100644 index 00000000000..0b91802ccd4 --- /dev/null +++ b/Solver_interface/package_info/Solver_interface/license.txt @@ -0,0 +1 @@ +LGPL (v3 or later) diff --git a/Solver_interface/package_info/Solver_interface/maintainer b/Solver_interface/package_info/Solver_interface/maintainer new file mode 100644 index 00000000000..5f89db128de --- /dev/null +++ b/Solver_interface/package_info/Solver_interface/maintainer @@ -0,0 +1,2 @@ +Gael.Guennebaud +Sebastien Loriot