User Manual

Authors

Simon Giraudot, Pierre Alliez, Frédéric Cazals, Gaël Guennebaud, Bruno Lévy, Marc Pouget, and Laurent Saboret

Several CGAL packages have to solve linear systems with dense or sparse matrices. This package provides concepts and models for that purpose.

We generally provide models using the Eigen library. Wrappers for the Eigen classes Eigen_matrix and Eigen_vector are also provided when needed.

It is straightforward to develop equivalent models for other solvers, for example those found in the Intel Math Kernel Library (MKL).

Matrix Diagonalization

The concept DiagonalizeTraits<T,dim> defines an interface for the diagonalization and computation of eigenvectors and eigenvalues of a symmetric matrix. T is the number type and dim is the dimension of the matrices and vector (set to 3 by default). We provide two models for this concept:

• Eigen_diagonalize_traits<T,dim> uses the Eigen library.
• Diagonalize_traits<T,dim> is an internal implementation that does not depend on another library.

Although both models achieve the same computation, Eigen_diagonalize_traits<T,dim> is faster and should thus be used if the Eigen library is available. The eigenvalues are stored in ascending order and eigenvectors are stored in accordance.

This is an example of an eigen decomposition of a matrix using this class:

File Solver_interface/diagonalize_matrix.cpp

#include <iostream>

#ifdef CGAL_EIGEN3_ENABLED
#include <CGAL/Eigen_diagonalize_traits.h>
#else
#include <CGAL/Diagonalize_traits.h>
#endif

typedef double                                          FT;
typedef CGAL::cpp11::array<FT, 6>                       Eigen_matrix;
typedef CGAL::cpp11::array<FT, 3>                       Eigen_vector;
typedef CGAL::cpp11::array<FT, 9>                       Eigen_three_vectors;

// If Eigen is enabled, use it, otherwise fallback to the internal model
#ifdef CGAL_EIGEN3_ENABLED
typedef CGAL::Eigen_diagonalize_traits<FT, 3>           Diagonalize_traits;
#else
typedef CGAL::Diagonalize_traits<FT, 3>                 Diagonalize_traits;
#endif

int main(void)
{
  Eigen_matrix covariance = {{ 0., 0., 0., 0., 0., 0. }};

  // Fill matrix with random numbers
  for(std::size_t i=0; i<6; ++i)
    covariance[i] = rand();

  Eigen_vector eigenvalues;
  Eigen_three_vectors eigenvectors;

  if(!(Diagonalize_traits::diagonalize_selfadjoint_covariance_matrix(covariance,
                                                                     eigenvalues,
                                                                     eigenvectors)))
  {
    std::cerr << "Error: cannot diagonalize matrix" << std::endl;
    return -1;
  }

  // Print result
  for(std::size_t i=0; i<3; ++i)
  {
    std::cout << "Eigenvalue " << i+1 << " = " << eigenvalues[i] << std::endl
              << "  with eigenvector [ ";
    for(std::size_t j=0; j<3; ++j)
      std::cout << eigenvectors[3*i + j] << " ";
    std::cout << "]" << std::endl;
  }

  return 0;
}

Singular Value Decomposition

The concept SvdTraits defines an interface for solving in the least square sense a linear system with a singular value decomposition. The field type is double. We provide the model Eigen_svd that uses the Eigen library.

Here is a simple example that shows how to handle matrices, vectors and this solver:

File Solver_interface/singular_value_decomposition.cpp

#ifdef CGAL_EIGEN3_ENABLED
#include <CGAL/Eigen_matrix.h>
#include <CGAL/Eigen_vector.h>
#include <CGAL/Eigen_svd.h>
typedef CGAL::Eigen_svd Svd;
#endif

typedef Svd::FT     FT;
typedef Svd::Vector Eigen_vector;
typedef Svd::Matrix Eigen_matrix;

int main(void)
{
  std::size_t degree = 3;

  Eigen_vector B(degree);
  Eigen_matrix M(degree, degree);

  // Fill B and M with random numbers
  for(std::size_t i=0; i<degree; ++i)
  {
    B.set(i, rand());
    for(std::size_t j = 0; j < degree; ++j)
      M.set(i, j, rand());
  }

  // Solve MX=B
  std::cout << Svd::solve(M, B) << std::endl;

  // Print result
  std::cout << "Solution of SVD = [ ";
  for(std::size_t i=0; i<degree; ++i)
    std::cout << B.vector()[i] << " ";
  std::cout << "]" << std::endl;

  return 0;
}

Sparse Solvers

We define 3 concepts for sparse linear algebra:

• SparseLinearAlgebraTraits_d
• SparseLinearAlgebraWithFactorTraits_d
• NormalEquationSparseLinearAlgebraTraits_d

An interface to the sparse solvers from the Eigen library is provided as a model for these 3 concepts through the class Eigen_solver_traits<T>. This solver traits class can be used for an iterative or a direct, symmetric or general sparse solvers. The specific solver to be used must be given as template parameter.

Each CGAL package using a sparse solver specifies which type of matrix and solver is required:

typedef CGAL::Eigen_sparse_matrix<double>::EigenType EigenMatrix;

//iterative general solver 
typedef CGAL::Eigen_solver_traits< Eigen::BiCGSTAB<EigenMatrix> > Iterative_general_solver;

//iterative symmetric solver 
typedef CGAL::Eigen_solver_traits< Eigen::ConjugateGradient<EigenMatrix> > Iterative_symmetric_solver;

//direct symmetric solver 
typedef CGAL::Eigen_solver_traits< Eigen::SimplicialCholesky<EigenMatrix> > Direct_symmetric_solver;

Here is an example that shows how to fill the sparse matrix and call the solver:

File Solver_interface/sparse_solvers.cpp

#include <CGAL/Eigen_solver_traits.h>
#include <CGAL/Eigen_matrix.h>

typedef CGAL::Eigen_solver_traits<> Eigen_solver;
typedef Eigen_solver::NT            FT;
typedef Eigen_solver::Matrix        Eigen_matrix;
typedef Eigen_solver::Vector        Eigen_vector;

int main(void)
{
  srand(static_cast<unsigned int>(time (NULL)));
  std::size_t degree = 3000;
  std::size_t nb_nonzero_coef = 100;

  Eigen_vector B(degree); // Zero vector
  Eigen_matrix A(degree);

  // Randomly make some coefficients of the matrix non-zero
  for(std::size_t i=0; i<nb_nonzero_coef; ++i)
  {
    int x = rand() % degree;
    int y = rand() % degree;

    FT value = rand() / (FT)RAND_MAX;

    A.add_coef(x, y, value);
    A.add_coef(y, x, value);
  }

  Eigen_vector X(degree);
  FT d;

  Eigen_solver solver;
  if(!(solver.linear_solver(A, B, X, d)))
  {
    std::cerr << "Error: linear solver failed" << std::endl;
    return -1;
  }

  std::cerr << "Linear solve succeeded" << std::endl;
  return 0;
}

Implementation History

This package is the result of the increasing needs for linear solvers in CGAL. The first packages that introduced the solver concepts were Triangulated Surface Mesh Parameterization Reference, Poisson Surface Reconstruction Reference and Estimation of Local Differential Properties of Point-Sampled Surfaces Reference. At that time, these packages were relying on Taucs, LAPACK, BLAS and OpenNL. Gaël Guennebaud then introduced new models using the Eigen library that became the only supported models by CGAL. Later on the packages Triangulated Surface Mesh Skeletonization and Triangulated Surface Mesh Deformation extended the existing concepts.

Simon Giraudot was responsible for gathering all concepts and classes, and also wrote this user manual with the help of Andreas Fabri.