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 #ifndef CERES_INTERNAL_SPARSE_CHOLESKY_H_ 32 #define CERES_INTERNAL_SPARSE_CHOLESKY_H_ 33 34 // This include must come before any #ifndef check on Ceres compile options. 35 #include "ceres/internal/port.h" 36 37 #include <memory> 38 #include "ceres/linear_solver.h" 39 #include "glog/logging.h" 40 41 namespace ceres { 42 namespace internal { 43 44 // An interface that abstracts away the internal details of various 45 // sparse linear algebra libraries and offers a simple API for solving 46 // symmetric positive definite linear systems using a sparse Cholesky 47 // factorization. 48 // 49 // Instances of SparseCholesky are expected to cache the symbolic 50 // factorization of the linear system. They do this on the first call 51 // to Factorize or FactorAndSolve. Subsequent calls to Factorize and 52 // FactorAndSolve are expected to have the same sparsity structure. 53 // 54 // Example usage: 55 // 56 // std::unique_ptr<SparseCholesky> 57 // sparse_cholesky(SparseCholesky::Create(SUITE_SPARSE, AMD)); 58 // 59 // CompressedRowSparseMatrix lhs = ...; 60 // std::string message; 61 // CHECK_EQ(sparse_cholesky->Factorize(&lhs, &message), LINEAR_SOLVER_SUCCESS); 62 // Vector rhs = ...; 63 // Vector solution = ...; 64 // CHECK_EQ(sparse_cholesky->Solve(rhs.data(), solution.data(), &message), 65 // LINEAR_SOLVER_SUCCESS); 66 67 class SparseCholesky { 68 public: 69 static std::unique_ptr<SparseCholesky> Create( 70 const LinearSolver::Options& options); 71 72 virtual ~SparseCholesky(); 73 74 // Due to the symmetry of the linear system, sparse linear algebra 75 // libraries only use one half of the input matrix. Whether it is 76 // the upper or the lower triangular part of the matrix depends on 77 // the library and the re-ordering strategy being used. This 78 // function tells the user the storage type expected of the input 79 // matrix for the sparse linear algebra library and reordering 80 // strategy used. 81 virtual CompressedRowSparseMatrix::StorageType StorageType() const = 0; 82 83 // Computes the numeric factorization of the given matrix. If this 84 // is the first call to Factorize, first the symbolic factorization 85 // will be computed and cached and the numeric factorization will be 86 // computed based on that. 87 // 88 // Subsequent calls to Factorize will use that symbolic 89 // factorization assuming that the sparsity of the matrix has 90 // remained constant. 91 virtual LinearSolverTerminationType Factorize( 92 CompressedRowSparseMatrix* lhs, std::string* message) = 0; 93 94 // Computes the solution to the equation 95 // 96 // lhs * solution = rhs 97 virtual LinearSolverTerminationType Solve(const double* rhs, 98 double* solution, 99 std::string* message) = 0; 100 101 // Convenience method which combines a call to Factorize and 102 // Solve. Solve is only called if Factorize returns 103 // LINEAR_SOLVER_SUCCESS. 104 virtual LinearSolverTerminationType FactorAndSolve( 105 CompressedRowSparseMatrix* lhs, 106 const double* rhs, 107 double* solution, 108 std::string* message); 109 110 }; 111 112 class IterativeRefiner; 113 114 // Computes an initial solution using the given instance of 115 // SparseCholesky, and then refines it using the IterativeRefiner. 116 class RefinedSparseCholesky : public SparseCholesky { 117 public: 118 RefinedSparseCholesky(std::unique_ptr<SparseCholesky> sparse_cholesky, 119 std::unique_ptr<IterativeRefiner> iterative_refiner); 120 virtual ~RefinedSparseCholesky(); 121 122 virtual CompressedRowSparseMatrix::StorageType StorageType() const; 123 virtual LinearSolverTerminationType Factorize( 124 CompressedRowSparseMatrix* lhs, std::string* message); 125 virtual LinearSolverTerminationType Solve(const double* rhs, 126 double* solution, 127 std::string* message); 128 129 private: 130 std::unique_ptr<SparseCholesky> sparse_cholesky_; 131 std::unique_ptr<IterativeRefiner> iterative_refiner_; 132 CompressedRowSparseMatrix* lhs_ = nullptr; 133 }; 134 135 } // namespace internal 136 } // namespace ceres 137 138 #endif // CERES_INTERNAL_SPARSE_CHOLESKY_H_ 139