block_sparse_matrix.h 8.3 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210
  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. //
  31. // Implementation of the SparseMatrix interface for block sparse
  32. // matrices.
  33. #ifndef CERES_INTERNAL_BLOCK_SPARSE_MATRIX_H_
  34. #define CERES_INTERNAL_BLOCK_SPARSE_MATRIX_H_
  35. #include <memory>
  36. #include <random>
  37. #include "ceres/block_structure.h"
  38. #include "ceres/compressed_row_sparse_matrix.h"
  39. #include "ceres/context_impl.h"
  40. #include "ceres/internal/disable_warnings.h"
  41. #include "ceres/internal/eigen.h"
  42. #include "ceres/internal/export.h"
  43. #include "ceres/sparse_matrix.h"
  44. namespace ceres::internal {
  45. class TripletSparseMatrix;
  46. // This class implements the SparseMatrix interface for storing and
  47. // manipulating block sparse matrices. The block structure is stored
  48. // in the CompressedRowBlockStructure object and one is needed to
  49. // initialize the matrix. For details on how the blocks structure of
  50. // the matrix is stored please see the documentation
  51. //
  52. // internal/ceres/block_structure.h
  53. //
  54. class CERES_NO_EXPORT BlockSparseMatrix final : public SparseMatrix {
  55. public:
  56. // Construct a block sparse matrix with a fully initialized
  57. // CompressedRowBlockStructure objected. The matrix takes over
  58. // ownership of this object and destroys it upon destruction.
  59. //
  60. // TODO(sameeragarwal): Add a function which will validate legal
  61. // CompressedRowBlockStructure objects.
  62. explicit BlockSparseMatrix(CompressedRowBlockStructure* block_structure,
  63. bool use_page_locked_memory = false);
  64. ~BlockSparseMatrix();
  65. BlockSparseMatrix(const BlockSparseMatrix&) = delete;
  66. void operator=(const BlockSparseMatrix&) = delete;
  67. // Implementation of SparseMatrix interface.
  68. void SetZero() override final;
  69. void SetZero(ContextImpl* context, int num_threads) override final;
  70. void RightMultiplyAndAccumulate(const double* x, double* y) const final;
  71. void RightMultiplyAndAccumulate(const double* x,
  72. double* y,
  73. ContextImpl* context,
  74. int num_threads) const final;
  75. void LeftMultiplyAndAccumulate(const double* x, double* y) const final;
  76. void LeftMultiplyAndAccumulate(const double* x,
  77. double* y,
  78. ContextImpl* context,
  79. int num_threads) const final;
  80. void SquaredColumnNorm(double* x) const final;
  81. void SquaredColumnNorm(double* x,
  82. ContextImpl* context,
  83. int num_threads) const final;
  84. void ScaleColumns(const double* scale) final;
  85. void ScaleColumns(const double* scale,
  86. ContextImpl* context,
  87. int num_threads) final;
  88. // Convert to CompressedRowSparseMatrix
  89. std::unique_ptr<CompressedRowSparseMatrix> ToCompressedRowSparseMatrix()
  90. const;
  91. // Create CompressedRowSparseMatrix corresponding to transposed matrix
  92. std::unique_ptr<CompressedRowSparseMatrix>
  93. ToCompressedRowSparseMatrixTranspose() const;
  94. // Copy values to CompressedRowSparseMatrix that has compatible structure
  95. void UpdateCompressedRowSparseMatrix(
  96. CompressedRowSparseMatrix* crs_matrix) const;
  97. // Copy values to CompressedRowSparseMatrix that has structure of transposed
  98. // matrix
  99. void UpdateCompressedRowSparseMatrixTranspose(
  100. CompressedRowSparseMatrix* crs_matrix) const;
  101. void ToDenseMatrix(Matrix* dense_matrix) const final;
  102. void ToTextFile(FILE* file) const final;
  103. void AddTransposeBlockStructure();
  104. // clang-format off
  105. int num_rows() const final { return num_rows_; }
  106. int num_cols() const final { return num_cols_; }
  107. int num_nonzeros() const final { return num_nonzeros_; }
  108. const double* values() const final { return values_; }
  109. double* mutable_values() final { return values_; }
  110. // clang-format on
  111. void ToTripletSparseMatrix(TripletSparseMatrix* matrix) const;
  112. const CompressedRowBlockStructure* block_structure() const;
  113. const CompressedRowBlockStructure* transpose_block_structure() const;
  114. // Append the contents of m to the bottom of this matrix. m must
  115. // have the same column blocks structure as this matrix.
  116. void AppendRows(const BlockSparseMatrix& m);
  117. // Delete the bottom delta_rows_blocks.
  118. void DeleteRowBlocks(int delta_row_blocks);
  119. static std::unique_ptr<BlockSparseMatrix> CreateDiagonalMatrix(
  120. const double* diagonal, const std::vector<Block>& column_blocks);
  121. struct RandomMatrixOptions {
  122. int num_row_blocks = 0;
  123. int min_row_block_size = 0;
  124. int max_row_block_size = 0;
  125. int num_col_blocks = 0;
  126. int min_col_block_size = 0;
  127. int max_col_block_size = 0;
  128. // 0 < block_density <= 1 is the probability of a block being
  129. // present in the matrix. A given random matrix will not have
  130. // precisely this density.
  131. double block_density = 0.0;
  132. // If col_blocks is non-empty, then the generated random matrix
  133. // has this block structure and the column related options in this
  134. // struct are ignored.
  135. std::vector<Block> col_blocks;
  136. };
  137. // Create a random BlockSparseMatrix whose entries are normally
  138. // distributed and whose structure is determined by
  139. // RandomMatrixOptions.
  140. static std::unique_ptr<BlockSparseMatrix> CreateRandomMatrix(
  141. const RandomMatrixOptions& options,
  142. std::mt19937& prng,
  143. bool use_page_locked_memory = false);
  144. private:
  145. double* AllocateValues(int size);
  146. void FreeValues(double*& values);
  147. const bool use_page_locked_memory_;
  148. int num_rows_;
  149. int num_cols_;
  150. int num_nonzeros_;
  151. int max_num_nonzeros_;
  152. double* values_;
  153. std::unique_ptr<CompressedRowBlockStructure> block_structure_;
  154. std::unique_ptr<CompressedRowBlockStructure> transpose_block_structure_;
  155. };
  156. // A number of algorithms like the SchurEliminator do not need
  157. // access to the full BlockSparseMatrix interface. They only
  158. // need read only access to the values array and the block structure.
  159. //
  160. // BlockSparseDataMatrix a struct that carries these two bits of
  161. // information
  162. class CERES_NO_EXPORT BlockSparseMatrixData {
  163. public:
  164. explicit BlockSparseMatrixData(const BlockSparseMatrix& m)
  165. : block_structure_(m.block_structure()), values_(m.values()){};
  166. BlockSparseMatrixData(const CompressedRowBlockStructure* block_structure,
  167. const double* values)
  168. : block_structure_(block_structure), values_(values) {}
  169. const CompressedRowBlockStructure* block_structure() const {
  170. return block_structure_;
  171. }
  172. const double* values() const { return values_; }
  173. private:
  174. const CompressedRowBlockStructure* block_structure_;
  175. const double* values_;
  176. };
  177. std::unique_ptr<CompressedRowBlockStructure> CreateTranspose(
  178. const CompressedRowBlockStructure& bs);
  179. } // namespace ceres::internal
  180. #include "ceres/internal/reenable_warnings.h"
  181. #endif // CERES_INTERNAL_BLOCK_SPARSE_MATRIX_H_