block_structure.h 7.6 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199
  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. // Block structure objects are used to carry information about the
  32. // dense block structure of sparse matrices. The BlockSparseMatrix
  33. // object uses the BlockStructure objects to keep track of the matrix
  34. // structure and operate upon it. This allows us to use more cache
  35. // friendly block oriented linear algebra operations on the matrix
  36. // instead of accessing it one scalar entry at a time.
  37. #ifndef CERES_INTERNAL_BLOCK_STRUCTURE_H_
  38. #define CERES_INTERNAL_BLOCK_STRUCTURE_H_
  39. #include <cstdint>
  40. #include <vector>
  41. #include "ceres/internal/export.h"
  42. // This file is being included into source files that are compiled with nvcc.
  43. // nvcc shipped with ubuntu 20.04 does not support some features of c++17,
  44. // including nested namespace definitions
  45. namespace ceres {
  46. namespace internal {
  47. using BlockSize = int32_t;
  48. struct CERES_NO_EXPORT Block {
  49. Block() = default;
  50. Block(int size_, int position_) noexcept : size(size_), position(position_) {}
  51. BlockSize size{-1};
  52. int position{-1}; // Position along the row/column.
  53. };
  54. inline bool operator==(const Block& left, const Block& right) noexcept {
  55. return (left.size == right.size) && (left.position == right.position);
  56. }
  57. struct CERES_NO_EXPORT Cell {
  58. Cell() = default;
  59. Cell(int block_id_, int position_) noexcept
  60. : block_id(block_id_), position(position_) {}
  61. // Column or row block id as the case maybe.
  62. int block_id{-1};
  63. // Where in the values array of the jacobian is this cell located.
  64. int position{-1};
  65. };
  66. // Order cell by their block_id;
  67. CERES_NO_EXPORT bool CellLessThan(const Cell& lhs, const Cell& rhs);
  68. struct CERES_NO_EXPORT CompressedList {
  69. CompressedList() = default;
  70. // Construct a CompressedList with the cells containing num_cells
  71. // entries.
  72. explicit CompressedList(int num_cells) noexcept : cells(num_cells) {}
  73. Block block;
  74. std::vector<Cell> cells;
  75. // Number of non-zeros in cells of this row block
  76. int nnz{-1};
  77. // Number of non-zeros in cells of this and every preceeding row block in
  78. // block-sparse matrix
  79. int cumulative_nnz{-1};
  80. };
  81. using CompressedRow = CompressedList;
  82. using CompressedColumn = CompressedList;
  83. // CompressedRowBlockStructure specifies the storage structure of a row block
  84. // sparse matrix.
  85. //
  86. // Consider the following matrix A:
  87. // A = [A_11 A_12 ...
  88. // A_21 A_22 ...
  89. // ...
  90. // A_m1 A_m2 ... ]
  91. //
  92. // A row block sparse matrix is a matrix where the following properties hold:
  93. // 1. The number of rows in every block A_ij and A_ik are the same.
  94. // 2. The number of columns in every block A_ij and A_kj are the same.
  95. // 3. The number of rows in A_ij and A_kj may be different (i != k).
  96. // 4. The number of columns in A_ij and A_ik may be different (j != k).
  97. // 5. Any block A_ij may be all 0s, in which case the block is not stored.
  98. //
  99. // The structure of the matrix is stored as follows:
  100. //
  101. // The `rows' array contains the following information for each row block:
  102. // - rows[i].block.size: The number of rows in each block A_ij in the row block.
  103. // - rows[i].block.position: The starting row in the full matrix A of the
  104. // row block i.
  105. // - rows[i].cells[j].block_id: The index into the `cols' array corresponding to
  106. // the non-zero blocks A_ij.
  107. // - rows[i].cells[j].position: The index in the `values' array for the contents
  108. // of block A_ij.
  109. //
  110. // The `cols' array contains the following information for block:
  111. // - cols[.].size: The number of columns spanned by the block.
  112. // - cols[.].position: The starting column in the full matrix A of the block.
  113. //
  114. //
  115. // Example of a row block sparse matrix:
  116. // block_id: | 0 |1|2 |3 |
  117. // rows[0]: [ 1 2 0 3 4 0 ]
  118. // [ 5 6 0 7 8 0 ]
  119. // rows[1]: [ 0 0 9 0 0 0 ]
  120. //
  121. // This matrix is stored as follows:
  122. //
  123. // There are four column blocks:
  124. // cols[0].size = 2
  125. // cols[0].position = 0
  126. // cols[1].size = 1
  127. // cols[1].position = 2
  128. // cols[2].size = 2
  129. // cols[2].position = 3
  130. // cols[3].size = 1
  131. // cols[3].position = 5
  132. // The first row block spans two rows, starting at row 0:
  133. // rows[0].block.size = 2 // This row block spans two rows.
  134. // rows[0].block.position = 0 // It starts at row 0.
  135. // rows[0] has two cells, at column blocks 0 and 2:
  136. // rows[0].cells[0].block_id = 0 // This cell is in column block 0.
  137. // rows[0].cells[0].position = 0 // See below for an explanation of this.
  138. // rows[0].cells[1].block_id = 2 // This cell is in column block 2.
  139. // rows[0].cells[1].position = 4 // See below for an explanation of this.
  140. //
  141. // The second row block spans two rows, starting at row 2:
  142. // rows[1].block.size = 1 // This row block spans one row.
  143. // rows[1].block.position = 2 // It starts at row 2.
  144. // rows[1] has one cell at column block 1:
  145. // rows[1].cells[0].block_id = 1 // This cell is in column block 1.
  146. // rows[1].cells[0].position = 8 // See below for an explanation of this.
  147. //
  148. // The values in each blocks are stored contiguously in row major order.
  149. // However, there is no unique way to order the blocks -- it is usually
  150. // optimized to promote cache coherent access, e.g. ordering it so that
  151. // Jacobian blocks of parameters of the same type are stored nearby.
  152. // This is one possible way to store the values of the blocks in a values array:
  153. // values = { 1, 2, 5, 6, 3, 4, 7, 8, 9 }
  154. // | | | | // The three blocks.
  155. // ^ rows[0].cells[0].position = 0
  156. // ^ rows[0].cells[1].position = 4
  157. // ^ rows[1].cells[0].position = 8
  158. struct CERES_NO_EXPORT CompressedRowBlockStructure {
  159. std::vector<Block> cols;
  160. std::vector<CompressedRow> rows;
  161. };
  162. struct CERES_NO_EXPORT CompressedColumnBlockStructure {
  163. std::vector<Block> rows;
  164. std::vector<CompressedColumn> cols;
  165. };
  166. inline int NumScalarEntries(const std::vector<Block>& blocks) {
  167. if (blocks.empty()) {
  168. return 0;
  169. }
  170. auto& block = blocks.back();
  171. return block.position + block.size;
  172. }
  173. std::vector<Block> Tail(const std::vector<Block>& blocks, int n);
  174. int SumSquaredSizes(const std::vector<Block>& blocks);
  175. } // namespace internal
  176. } // namespace ceres
  177. #endif // CERES_INTERNAL_BLOCK_STRUCTURE_H_