sparse_cholesky.h 6.0 KB

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