internal/ceres/eigensparse.cc

   1 // Ceres Solver - A fast non-linear least squares minimizer
   2 // Copyright 2017 Google Inc. All rights reserved.
   3 // http://ceres-solver.org/
   4 //
   5 // Redistribution and use in source and binary forms, with or without
   6 // modification, are permitted provided that the following conditions are met:
   7 //
   8 // * Redistributions of source code must retain the above copyright notice,
   9 //   this list of conditions and the following disclaimer.
  10 // * Redistributions in binary form must reproduce the above copyright notice,
  11 //   this list of conditions and the following disclaimer in the documentation
  12 //   and/or other materials provided with the distribution.
  13 // * Neither the name of Google Inc. nor the names of its contributors may be
  14 //   used to endorse or promote products derived from this software without
  15 //   specific prior written permission.
  16 //
  17 // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
  18 // AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
  19 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
  20 // ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
  21 // LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
  22 // CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
  23 // SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
  24 // INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
  25 // CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
  26 // ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
  27 // POSSIBILITY OF SUCH DAMAGE.
  28 //
  29 // Author: sameeragarwal@google.com (Sameer Agarwal)
  30
  31 #include "ceres/eigensparse.h"
  32
  33 #ifdef CERES_USE_EIGEN_SPARSE
  34
  35 #include <sstream>
  36 #include "Eigen/SparseCholesky"
  37 #include "Eigen/SparseCore"
  38 #include "ceres/compressed_row_sparse_matrix.h"
  39 #include "ceres/linear_solver.h"
  40
  41 namespace ceres {
  42 namespace internal {
  43
  44 template <typename Solver>
  45 class EigenSparseCholeskyTemplate : public EigenSparseCholesky {
  46  public:
  47   EigenSparseCholeskyTemplate() : analyzed_(false) {}
  48   virtual ~EigenSparseCholeskyTemplate() {}
  49   virtual CompressedRowSparseMatrix::StorageType StorageType() const {
  50     return CompressedRowSparseMatrix::LOWER_TRIANGULAR;
  51   }
  52
  53   virtual LinearSolverTerminationType Factorize(
  54       const Eigen::SparseMatrix<double>& lhs, std::string* message) {
  55     if (!analyzed_) {
  56       solver_.analyzePattern(lhs);
  57
  58       if (VLOG_IS_ON(2)) {
  59         std::stringstream ss;
  60         solver_.dumpMemory(ss);
  61         VLOG(2) << "Symbolic Analysis\n" << ss.str();
  62       }
  63
  64       if (solver_.info() != Eigen::Success) {
  65         *message = "Eigen failure. Unable to find symbolic factorization.";
  66         return LINEAR_SOLVER_FATAL_ERROR;
  67       }
  68
  69       analyzed_ = true;
  70     }
  71
  72     solver_.factorize(lhs);
  73     if (solver_.info() != Eigen::Success) {
  74       *message = "Eigen failure. Unable to find numeric factorization.";
  75       return LINEAR_SOLVER_FAILURE;
  76     }
  77     return LINEAR_SOLVER_SUCCESS;
  78   }
  79
  80   virtual LinearSolverTerminationType Solve(const double* rhs,
  81                                             double* solution,
  82                                             std::string* message) {
  83     CHECK(analyzed_) << "Solve called without a call to Factorize first.";
  84
  85     VectorRef(solution, solver_.cols()) =
  86         solver_.solve(ConstVectorRef(rhs, solver_.cols()));
  87     if (solver_.info() != Eigen::Success) {
  88       *message = "Eigen failure. Unable to do triangular solve.";
  89       return LINEAR_SOLVER_FAILURE;
  90     }
  91     return LINEAR_SOLVER_SUCCESS;
  92   }
  93
  94   virtual LinearSolverTerminationType Factorize(CompressedRowSparseMatrix* lhs,
  95                                                 std::string* message) {
  96     CHECK_EQ(lhs->storage_type(), StorageType());
  97     Eigen::MappedSparseMatrix<double, Eigen::ColMajor> eigen_lhs(
  98         lhs->num_rows(),
  99         lhs->num_rows(),
 100         lhs->num_nonzeros(),
 101         lhs->mutable_rows(),
 102         lhs->mutable_cols(),
 103         lhs->mutable_values());
 104     return Factorize(eigen_lhs, message);
 105   }
 106
 107  private:
 108   bool analyzed_;
 109   Solver solver_;
 110 };
 111
 112 EigenSparseCholesky* EigenSparseCholesky::Create(
 113     const OrderingType ordering_type) {
 114   // The preprocessor gymnastics here are dealing with the fact that
 115   // before version 3.2.2, Eigen did not support a third template
 116   // parameter to specify the ordering and it always defaults to AMD.
 117 #if EIGEN_VERSION_AT_LEAST(3, 2, 2)
 118   typedef Eigen::SimplicialLDLT<Eigen::SparseMatrix<double>,
 119                                 Eigen::Upper,
 120                                 Eigen::AMDOrdering<int> >
 121       WithAMDOrdering;
 122   typedef Eigen::SimplicialLDLT<Eigen::SparseMatrix<double>,
 123                                 Eigen::Upper,
 124                                 Eigen::NaturalOrdering<int> >
 125       WithNaturalOrdering;
 126   if (ordering_type == AMD) {
 127     return new EigenSparseCholeskyTemplate<WithAMDOrdering>();
 128   } else {
 129     return new EigenSparseCholeskyTemplate<WithNaturalOrdering>();
 130   }
 131 #else
 132   typedef Eigen::SimplicialLDLT<Eigen::SparseMatrix<double>, Eigen::Upper>
 133       WithAMDOrdering;
 134   return new EigenSparseCholeskyTemplate<WithAMDOrdering>();
 135 #endif
 136 }
 137
 138 EigenSparseCholesky::~EigenSparseCholesky() {}
 139
 140 }  // namespace internal
 141 }  // namespace ceres
 142
 143 #endif  // CERES_USE_EIGEN_SPARSE