compressed_row_sparse_matrix.h 9.7 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249
  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_COMPRESSED_ROW_SPARSE_MATRIX_H_
  31. #define CERES_INTERNAL_COMPRESSED_ROW_SPARSE_MATRIX_H_
  32. #include <memory>
  33. #include <random>
  34. #include <vector>
  35. #include "ceres/block_structure.h"
  36. #include "ceres/internal/disable_warnings.h"
  37. #include "ceres/internal/export.h"
  38. #include "ceres/sparse_matrix.h"
  39. #include "ceres/types.h"
  40. #include "glog/logging.h"
  41. namespace ceres {
  42. struct CRSMatrix;
  43. namespace internal {
  44. class ContextImpl;
  45. class TripletSparseMatrix;
  46. class CERES_NO_EXPORT CompressedRowSparseMatrix : public SparseMatrix {
  47. public:
  48. enum class StorageType {
  49. UNSYMMETRIC,
  50. // Matrix is assumed to be symmetric but only the lower triangular
  51. // part of the matrix is stored.
  52. LOWER_TRIANGULAR,
  53. // Matrix is assumed to be symmetric but only the upper triangular
  54. // part of the matrix is stored.
  55. UPPER_TRIANGULAR
  56. };
  57. // Create a matrix with the same content as the TripletSparseMatrix
  58. // input. We assume that input does not have any repeated
  59. // entries.
  60. //
  61. // The storage type of the matrix is set to UNSYMMETRIC.
  62. static std::unique_ptr<CompressedRowSparseMatrix> FromTripletSparseMatrix(
  63. const TripletSparseMatrix& input);
  64. // Create a matrix with the same content as the TripletSparseMatrix
  65. // input transposed. We assume that input does not have any repeated
  66. // entries.
  67. //
  68. // The storage type of the matrix is set to UNSYMMETRIC.
  69. static std::unique_ptr<CompressedRowSparseMatrix>
  70. FromTripletSparseMatrixTransposed(const TripletSparseMatrix& input);
  71. // Use this constructor only if you know what you are doing. This
  72. // creates a "blank" matrix with the appropriate amount of memory
  73. // allocated. However, the object itself is in an inconsistent state
  74. // as the rows and cols matrices do not match the values of
  75. // num_rows, num_cols and max_num_nonzeros.
  76. //
  77. // The use case for this constructor is that when the user knows the
  78. // size of the matrix to begin with and wants to update the layout
  79. // manually, instead of going via the indirect route of first
  80. // constructing a TripletSparseMatrix, which leads to more than
  81. // double the peak memory usage.
  82. //
  83. // The storage type is set to UNSYMMETRIC.
  84. CompressedRowSparseMatrix(int num_rows, int num_cols, int max_num_nonzeros);
  85. // Build a square sparse diagonal matrix with num_rows rows and
  86. // columns. The diagonal m(i,i) = diagonal(i);
  87. //
  88. // The storage type is set to UNSYMMETRIC
  89. CompressedRowSparseMatrix(const double* diagonal, int num_rows);
  90. // SparseMatrix interface.
  91. ~CompressedRowSparseMatrix() override;
  92. void SetZero() final;
  93. void RightMultiplyAndAccumulate(const double* x, double* y) const final;
  94. void RightMultiplyAndAccumulate(const double* x,
  95. double* y,
  96. ContextImpl* context,
  97. int num_threads) const final;
  98. void LeftMultiplyAndAccumulate(const double* x, double* y) const final;
  99. void SquaredColumnNorm(double* x) const final;
  100. void ScaleColumns(const double* scale) final;
  101. void ToDenseMatrix(Matrix* dense_matrix) const final;
  102. void ToTextFile(FILE* file) const final;
  103. int num_rows() const final { return num_rows_; }
  104. int num_cols() const final { return num_cols_; }
  105. int num_nonzeros() const final { return rows_[num_rows_]; }
  106. const double* values() const final { return values_.data(); }
  107. double* mutable_values() final { return values_.data(); }
  108. // Delete the bottom delta_rows.
  109. // num_rows -= delta_rows
  110. void DeleteRows(int delta_rows);
  111. // Append the contents of m to the bottom of this matrix. m must
  112. // have the same number of columns as this matrix.
  113. void AppendRows(const CompressedRowSparseMatrix& m);
  114. void ToCRSMatrix(CRSMatrix* matrix) const;
  115. std::unique_ptr<CompressedRowSparseMatrix> Transpose() const;
  116. // Destructive array resizing method.
  117. void SetMaxNumNonZeros(int num_nonzeros);
  118. // Non-destructive array resizing method.
  119. void set_num_rows(const int num_rows) { num_rows_ = num_rows; }
  120. void set_num_cols(const int num_cols) { num_cols_ = num_cols; }
  121. // Low level access methods that expose the structure of the matrix.
  122. const int* cols() const { return cols_.data(); }
  123. int* mutable_cols() { return cols_.data(); }
  124. const int* rows() const { return rows_.data(); }
  125. int* mutable_rows() { return rows_.data(); }
  126. StorageType storage_type() const { return storage_type_; }
  127. void set_storage_type(const StorageType storage_type) {
  128. storage_type_ = storage_type;
  129. }
  130. const std::vector<Block>& row_blocks() const { return row_blocks_; }
  131. std::vector<Block>* mutable_row_blocks() { return &row_blocks_; }
  132. const std::vector<Block>& col_blocks() const { return col_blocks_; }
  133. std::vector<Block>* mutable_col_blocks() { return &col_blocks_; }
  134. // Create a block diagonal CompressedRowSparseMatrix with the given
  135. // block structure. The individual blocks are assumed to be laid out
  136. // contiguously in the diagonal array, one block at a time.
  137. static std::unique_ptr<CompressedRowSparseMatrix> CreateBlockDiagonalMatrix(
  138. const double* diagonal, const std::vector<Block>& blocks);
  139. // Options struct to control the generation of random block sparse
  140. // matrices in compressed row sparse format.
  141. //
  142. // The random matrix generation proceeds as follows.
  143. //
  144. // First the row and column block structure is determined by
  145. // generating random row and column block sizes that lie within the
  146. // given bounds.
  147. //
  148. // Then we walk the block structure of the resulting matrix, and with
  149. // probability block_density determine whether they are structurally
  150. // zero or not. If the answer is no, then we generate entries for the
  151. // block which are distributed normally.
  152. struct RandomMatrixOptions {
  153. // Type of matrix to create.
  154. //
  155. // If storage_type is UPPER_TRIANGULAR (LOWER_TRIANGULAR), then
  156. // create a square symmetric matrix with just the upper triangular
  157. // (lower triangular) part. In this case, num_col_blocks,
  158. // min_col_block_size and max_col_block_size will be ignored and
  159. // assumed to be equal to the corresponding row settings.
  160. StorageType storage_type = StorageType::UNSYMMETRIC;
  161. int num_row_blocks = 0;
  162. int min_row_block_size = 0;
  163. int max_row_block_size = 0;
  164. int num_col_blocks = 0;
  165. int min_col_block_size = 0;
  166. int max_col_block_size = 0;
  167. // 0 < block_density <= 1 is the probability of a block being
  168. // present in the matrix. A given random matrix will not have
  169. // precisely this density.
  170. double block_density = 0.0;
  171. };
  172. // Create a random CompressedRowSparseMatrix whose entries are
  173. // normally distributed and whose structure is determined by
  174. // RandomMatrixOptions.
  175. static std::unique_ptr<CompressedRowSparseMatrix> CreateRandomMatrix(
  176. RandomMatrixOptions options, std::mt19937& prng);
  177. private:
  178. static std::unique_ptr<CompressedRowSparseMatrix> FromTripletSparseMatrix(
  179. const TripletSparseMatrix& input, bool transpose);
  180. int num_rows_;
  181. int num_cols_;
  182. std::vector<int> rows_;
  183. std::vector<int> cols_;
  184. std::vector<double> values_;
  185. StorageType storage_type_;
  186. // If the matrix has an underlying block structure, then it can also
  187. // carry with it row and column block sizes. This is auxiliary and
  188. // optional information for use by algorithms operating on the
  189. // matrix. The class itself does not make use of this information in
  190. // any way.
  191. std::vector<Block> row_blocks_;
  192. std::vector<Block> col_blocks_;
  193. };
  194. inline std::ostream& operator<<(std::ostream& s,
  195. CompressedRowSparseMatrix::StorageType type) {
  196. switch (type) {
  197. case CompressedRowSparseMatrix::StorageType::UNSYMMETRIC:
  198. s << "UNSYMMETRIC";
  199. break;
  200. case CompressedRowSparseMatrix::StorageType::UPPER_TRIANGULAR:
  201. s << "UPPER_TRIANGULAR";
  202. break;
  203. case CompressedRowSparseMatrix::StorageType::LOWER_TRIANGULAR:
  204. s << "LOWER_TRIANGULAR";
  205. break;
  206. default:
  207. s << "UNKNOWN CompressedRowSparseMatrix::StorageType";
  208. }
  209. return s;
  210. }
  211. } // namespace internal
  212. } // namespace ceres
  213. #include "ceres/internal/reenable_warnings.h"
  214. #endif // CERES_INTERNAL_COMPRESSED_ROW_SPARSE_MATRIX_H_