1 // Ceres Solver - A fast non-linear least squares minimizer 2 // Copyright 2012 Google Inc. All rights reserved. 3 // http://code.google.com/p/ceres-solver/ 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 (at) google.com (Sameer Agarwal) 30 31 #include <algorithm> 32 #include "ceres/internal/port.h" 33 #include "ceres/suitesparse.h" 34 #include "ceres/triplet_sparse_matrix.h" 35 #include "glog/logging.h" 36 #include "gtest/gtest.h" 37 38 namespace ceres { 39 namespace internal { 40 41 TEST(SuiteSparse, BlockPermutationToScalarPermutation) { 42 vector<int> blocks; 43 // Block structure 44 // 0 --1- ---2--- ---3--- 4 45 // [0, 1, 2, 3, 4, 5, 6, 7, 8, 9] 46 blocks.push_back(1); 47 blocks.push_back(2); 48 blocks.push_back(3); 49 blocks.push_back(3); 50 blocks.push_back(1); 51 52 // Block ordering 53 // [1, 0, 2, 4, 5] 54 vector<int> block_ordering; 55 block_ordering.push_back(1); 56 block_ordering.push_back(0); 57 block_ordering.push_back(2); 58 block_ordering.push_back(4); 59 block_ordering.push_back(3); 60 61 // Expected ordering 62 // [1, 2, 0, 3, 4, 5, 9, 6, 7, 8] 63 vector<int> expected_scalar_ordering; 64 expected_scalar_ordering.push_back(1); 65 expected_scalar_ordering.push_back(2); 66 expected_scalar_ordering.push_back(0); 67 expected_scalar_ordering.push_back(3); 68 expected_scalar_ordering.push_back(4); 69 expected_scalar_ordering.push_back(5); 70 expected_scalar_ordering.push_back(9); 71 expected_scalar_ordering.push_back(6); 72 expected_scalar_ordering.push_back(7); 73 expected_scalar_ordering.push_back(8); 74 75 vector<int> scalar_ordering; 76 SuiteSparse::BlockOrderingToScalarOrdering(blocks, 77 block_ordering, 78 &scalar_ordering); 79 EXPECT_EQ(scalar_ordering.size(), expected_scalar_ordering.size()); 80 for (int i = 0; i < expected_scalar_ordering.size(); ++i) { 81 EXPECT_EQ(scalar_ordering[i], expected_scalar_ordering[i]); 82 } 83 } 84 85 // Helper function to fill the sparsity pattern of a TripletSparseMatrix. 86 int FillBlock(const vector<int>& row_blocks, 87 const vector<int>& col_blocks, 88 const int row_block_id, 89 const int col_block_id, 90 int* rows, 91 int* cols) { 92 int row_pos = 0; 93 for (int i = 0; i < row_block_id; ++i) { 94 row_pos += row_blocks[i]; 95 } 96 97 int col_pos = 0; 98 for (int i = 0; i < col_block_id; ++i) { 99 col_pos += col_blocks[i]; 100 } 101 102 int offset = 0; 103 for (int r = 0; r < row_blocks[row_block_id]; ++r) { 104 for (int c = 0; c < col_blocks[col_block_id]; ++c, ++offset) { 105 rows[offset] = row_pos + r; 106 cols[offset] = col_pos + c; 107 } 108 } 109 return offset; 110 } 111 112 TEST(SuiteSparse, ScalarMatrixToBlockMatrix) { 113 // Block sparsity. 114 // 115 // [1 2 3 2] 116 // [1] x x 117 // [2] x x 118 // [2] x x 119 // num_nonzeros = 1 + 3 + 4 + 4 + 1 + 2 = 15 120 121 vector<int> col_blocks; 122 col_blocks.push_back(1); 123 col_blocks.push_back(2); 124 col_blocks.push_back(3); 125 col_blocks.push_back(2); 126 127 vector<int> row_blocks; 128 row_blocks.push_back(1); 129 row_blocks.push_back(2); 130 row_blocks.push_back(2); 131 132 TripletSparseMatrix tsm(5, 8, 18); 133 int* rows = tsm.mutable_rows(); 134 int* cols = tsm.mutable_cols(); 135 fill(tsm.mutable_values(), tsm.mutable_values() + 18, 1.0); 136 int offset = 0; 137 138 #define CERES_TEST_FILL_BLOCK(row_block_id, col_block_id) \ 139 offset += FillBlock(row_blocks, col_blocks, \ 140 row_block_id, col_block_id, \ 141 rows + offset, cols + offset); 142 143 CERES_TEST_FILL_BLOCK(0, 0); 144 CERES_TEST_FILL_BLOCK(2, 0); 145 CERES_TEST_FILL_BLOCK(1, 1); 146 CERES_TEST_FILL_BLOCK(2, 1); 147 CERES_TEST_FILL_BLOCK(0, 2); 148 CERES_TEST_FILL_BLOCK(1, 3); 149 #undef CERES_TEST_FILL_BLOCK 150 151 tsm.set_num_nonzeros(offset); 152 153 SuiteSparse ss; 154 scoped_ptr<cholmod_sparse> ccsm(ss.CreateSparseMatrix(&tsm)); 155 156 vector<int> expected_block_rows; 157 expected_block_rows.push_back(0); 158 expected_block_rows.push_back(2); 159 expected_block_rows.push_back(1); 160 expected_block_rows.push_back(2); 161 expected_block_rows.push_back(0); 162 expected_block_rows.push_back(1); 163 164 vector<int> expected_block_cols; 165 expected_block_cols.push_back(0); 166 expected_block_cols.push_back(2); 167 expected_block_cols.push_back(4); 168 expected_block_cols.push_back(5); 169 expected_block_cols.push_back(6); 170 171 vector<int> block_rows; 172 vector<int> block_cols; 173 SuiteSparse::ScalarMatrixToBlockMatrix(ccsm.get(), 174 row_blocks, 175 col_blocks, 176 &block_rows, 177 &block_cols); 178 179 EXPECT_EQ(block_cols.size(), expected_block_cols.size()); 180 EXPECT_EQ(block_rows.size(), expected_block_rows.size()); 181 182 for (int i = 0; i < expected_block_cols.size(); ++i) { 183 EXPECT_EQ(block_cols[i], expected_block_cols[i]); 184 } 185 186 for (int i = 0; i < expected_block_rows.size(); ++i) { 187 EXPECT_EQ(block_rows[i], expected_block_rows[i]); 188 } 189 190 ss.Free(ccsm.release()); 191 } 192 193 } // namespace internal 194 } // namespace ceres 195