block_jacobian_writer.h 6.1 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146
  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: keir@google.com (Keir Mierle)
  30. //
  31. // A jacobian writer that writes to block sparse matrices. The "writer" name is
  32. // misleading, since the Write() operation on the block jacobian writer does not
  33. // write anything. Instead, the Prepare() method on the BlockEvaluatePreparers
  34. // makes a jacobians array which has direct pointers into the block sparse
  35. // jacobian. When the cost function is evaluated, the jacobian blocks get placed
  36. // directly in their final location.
  37. #ifndef CERES_INTERNAL_BLOCK_JACOBIAN_WRITER_H_
  38. #define CERES_INTERNAL_BLOCK_JACOBIAN_WRITER_H_
  39. #include <memory>
  40. #include <vector>
  41. #include "ceres/evaluator.h"
  42. #include "ceres/internal/export.h"
  43. namespace ceres::internal {
  44. class BlockEvaluatePreparer;
  45. class Program;
  46. class SparseMatrix;
  47. // TODO(sameeragarwal): This class needs documentation.
  48. class CERES_NO_EXPORT BlockJacobianWriter {
  49. public:
  50. // Pre-computes positions of cells in block-sparse jacobian.
  51. // Two possible memory layouts are implemented:
  52. // - Non-partitioned case
  53. // - Partitioned case (for Schur type linear solver)
  54. //
  55. // In non-partitioned case, cells are stored sequentially in the
  56. // lexicographic order of (row block id, column block id).
  57. //
  58. // In the case of partitoned matrix, cells of each sub-matrix (E and F) are
  59. // stored sequentially in the lexicographic order of (row block id, column
  60. // block id) and cells from E sub-matrix precede cells from F sub-matrix.
  61. BlockJacobianWriter(const Evaluator::Options& options, Program* program);
  62. // JacobianWriter interface.
  63. // Create evaluate preparers that point directly into the final jacobian.
  64. // This makes the final Write() a nop.
  65. std::unique_ptr<BlockEvaluatePreparer[]> CreateEvaluatePreparers(
  66. unsigned num_threads);
  67. std::unique_ptr<SparseMatrix> CreateJacobian() const;
  68. void Write(int /* residual_id */,
  69. int /* residual_offset */,
  70. double** /* jacobians */,
  71. SparseMatrix* /* jacobian */) {
  72. // This is a noop since the blocks were written directly into their final
  73. // position by the outside evaluate call, thanks to the jacobians array
  74. // prepared by the BlockEvaluatePreparers.
  75. }
  76. private:
  77. Evaluator::Options options_;
  78. Program* program_;
  79. // Stores the position of each residual / parameter jacobian.
  80. //
  81. // The block sparse matrix that this writer writes to is stored as a set of
  82. // contiguous dense blocks, one after each other; see BlockSparseMatrix. The
  83. // "double* values_" member of the block sparse matrix contains all of these
  84. // blocks. Given a pointer to the first element of a block and the size of
  85. // that block, it's possible to write to it.
  86. //
  87. // In the case of a block sparse jacobian, the jacobian writer needs a way to
  88. // find the offset in the values_ array of each residual/parameter jacobian
  89. // block.
  90. //
  91. // That is the purpose of jacobian_layout_.
  92. //
  93. // In particular, jacobian_layout_[i][j] is the offset in the values_ array of
  94. // the derivative of residual block i with respect to the parameter block at
  95. // active argument position j.
  96. //
  97. // The active qualifier means that non-active parameters do not count. Care
  98. // must be taken when indexing into jacobian_layout_ to account for this.
  99. // Consider a single residual example:
  100. //
  101. // r(x, y, z)
  102. //
  103. // with r in R^3, x in R^4, y in R^2, and z in R^5.
  104. // Take y as a constant (non-active) parameter.
  105. // Take r as residual number 0.
  106. //
  107. // In this case, the active arguments are only (x, z), so the active argument
  108. // position for x is 0, and the active argument position for z is 1. This is
  109. // similar to thinking of r as taking only 2 parameters:
  110. //
  111. // r(x, z)
  112. //
  113. // There are only 2 jacobian blocks: dr/dx and dr/dz. jacobian_layout_ would
  114. // have the following contents:
  115. //
  116. // jacobian_layout_[0] = { 0, 12 }
  117. //
  118. // which indicates that dr/dx is located at values_[0], and dr/dz is at
  119. // values_[12]. See BlockEvaluatePreparer::Prepare()'s comments about 'j'.
  120. std::vector<int*> jacobian_layout_;
  121. // The pointers in jacobian_layout_ point directly into this vector.
  122. std::vector<int> jacobian_layout_storage_;
  123. // The constructor computes the layout of the Jacobian, and this bool keeps
  124. // track of whether the computation of the layout completed successfully or
  125. // not, if it is false, then jacobian_layout and jacobian_layout_storage are
  126. // both in an invalid state.
  127. bool jacobian_layout_is_valid_ = false;
  128. };
  129. } // namespace ceres::internal
  130. #endif // CERES_INTERNAL_BLOCK_JACOBIAN_WRITER_H_