compressed_col_sparse_matrix_utils.cc 4.6 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107
  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/compressed_col_sparse_matrix_utils.h"
  31. #include <algorithm>
  32. #include <vector>
  33. #include "ceres/internal/export.h"
  34. #include "glog/logging.h"
  35. namespace ceres::internal {
  36. void CompressedColumnScalarMatrixToBlockMatrix(
  37. const int* scalar_rows,
  38. const int* scalar_cols,
  39. const std::vector<Block>& row_blocks,
  40. const std::vector<Block>& col_blocks,
  41. std::vector<int>* block_rows,
  42. std::vector<int>* block_cols) {
  43. CHECK(block_rows != nullptr);
  44. CHECK(block_cols != nullptr);
  45. block_rows->clear();
  46. block_cols->clear();
  47. const int num_col_blocks = col_blocks.size();
  48. // This loop extracts the block sparsity of the scalar sparse matrix
  49. // It does so by iterating over the columns, but only considering
  50. // the columns corresponding to the first element of each column
  51. // block. Within each column, the inner loop iterates over the rows,
  52. // and detects the presence of a row block by checking for the
  53. // presence of a non-zero entry corresponding to its first element.
  54. block_cols->push_back(0);
  55. int c = 0;
  56. for (int col_block = 0; col_block < num_col_blocks; ++col_block) {
  57. int column_size = 0;
  58. for (int idx = scalar_cols[c]; idx < scalar_cols[c + 1]; ++idx) {
  59. auto it = std::lower_bound(row_blocks.begin(),
  60. row_blocks.end(),
  61. scalar_rows[idx],
  62. [](const Block& block, double value) {
  63. return block.position < value;
  64. });
  65. // Since we are using lower_bound, it will return the row id where the row
  66. // block starts. For everything but the first row of the block, where
  67. // these values will be the same, we can skip, as we only need the first
  68. // row to detect the presence of the block.
  69. //
  70. // For rows all but the first row in the last row block, lower_bound will
  71. // return row_blocks_.end(), but those can be skipped like the rows in
  72. // other row blocks too.
  73. if (it == row_blocks.end() || it->position != scalar_rows[idx]) {
  74. continue;
  75. }
  76. block_rows->push_back(it - row_blocks.begin());
  77. ++column_size;
  78. }
  79. block_cols->push_back(block_cols->back() + column_size);
  80. c += col_blocks[col_block].size;
  81. }
  82. }
  83. void BlockOrderingToScalarOrdering(const std::vector<Block>& blocks,
  84. const std::vector<int>& block_ordering,
  85. std::vector<int>* scalar_ordering) {
  86. CHECK_EQ(blocks.size(), block_ordering.size());
  87. const int num_blocks = blocks.size();
  88. scalar_ordering->resize(NumScalarEntries(blocks));
  89. int cursor = 0;
  90. for (int i = 0; i < num_blocks; ++i) {
  91. const int block_id = block_ordering[i];
  92. const int block_size = blocks[block_id].size;
  93. int block_position = blocks[block_id].position;
  94. for (int j = 0; j < block_size; ++j) {
  95. (*scalar_ordering)[cursor++] = block_position++;
  96. }
  97. }
  98. }
  99. } // namespace ceres::internal