sparse_normal_cholesky_solver.cc 4.2 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112
  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. #include "ceres/sparse_normal_cholesky_solver.h"
  31. #include <algorithm>
  32. #include <cstring>
  33. #include <ctime>
  34. #include <memory>
  35. #include "ceres/block_sparse_matrix.h"
  36. #include "ceres/inner_product_computer.h"
  37. #include "ceres/internal/eigen.h"
  38. #include "ceres/iterative_refiner.h"
  39. #include "ceres/linear_solver.h"
  40. #include "ceres/sparse_cholesky.h"
  41. #include "ceres/triplet_sparse_matrix.h"
  42. #include "ceres/types.h"
  43. #include "ceres/wall_time.h"
  44. namespace ceres::internal {
  45. SparseNormalCholeskySolver::SparseNormalCholeskySolver(
  46. const LinearSolver::Options& options)
  47. : options_(options) {
  48. sparse_cholesky_ = SparseCholesky::Create(options);
  49. }
  50. SparseNormalCholeskySolver::~SparseNormalCholeskySolver() = default;
  51. LinearSolver::Summary SparseNormalCholeskySolver::SolveImpl(
  52. BlockSparseMatrix* A,
  53. const double* b,
  54. const LinearSolver::PerSolveOptions& per_solve_options,
  55. double* x) {
  56. EventLogger event_logger("SparseNormalCholeskySolver::Solve");
  57. LinearSolver::Summary summary;
  58. summary.num_iterations = 1;
  59. summary.termination_type = LinearSolverTerminationType::SUCCESS;
  60. summary.message = "Success.";
  61. const int num_cols = A->num_cols();
  62. VectorRef xref(x, num_cols);
  63. xref.setZero();
  64. rhs_.resize(num_cols);
  65. rhs_.setZero();
  66. A->LeftMultiplyAndAccumulate(b, rhs_.data());
  67. event_logger.AddEvent("Compute RHS");
  68. if (per_solve_options.D != nullptr) {
  69. // Temporarily append a diagonal block to the A matrix, but undo
  70. // it before returning the matrix to the user.
  71. std::unique_ptr<BlockSparseMatrix> regularizer =
  72. BlockSparseMatrix::CreateDiagonalMatrix(per_solve_options.D,
  73. A->block_structure()->cols);
  74. event_logger.AddEvent("Diagonal");
  75. A->AppendRows(*regularizer);
  76. event_logger.AddEvent("Append");
  77. }
  78. event_logger.AddEvent("Append Rows");
  79. if (inner_product_computer_.get() == nullptr) {
  80. inner_product_computer_ =
  81. InnerProductComputer::Create(*A, sparse_cholesky_->StorageType());
  82. event_logger.AddEvent("InnerProductComputer::Create");
  83. }
  84. inner_product_computer_->Compute();
  85. event_logger.AddEvent("InnerProductComputer::Compute");
  86. if (per_solve_options.D != nullptr) {
  87. A->DeleteRowBlocks(A->block_structure()->cols.size());
  88. }
  89. summary.termination_type = sparse_cholesky_->FactorAndSolve(
  90. inner_product_computer_->mutable_result(),
  91. rhs_.data(),
  92. x,
  93. &summary.message);
  94. event_logger.AddEvent("SparseCholesky::FactorAndSolve");
  95. return summary;
  96. }
  97. } // namespace ceres::internal