1 // Ceres Solver - A fast non-linear least squares minimizer 2 // Copyright 2013 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: keir (at) google.com (Keir Mierle) 30 31 #include "ceres/blas.h" 32 33 #include "gtest/gtest.h" 34 #include "ceres/internal/eigen.h" 35 36 namespace ceres { 37 namespace internal { 38 39 TEST(BLAS, MatrixMatrixMultiply) { 40 const double kTolerance = 1e-16; 41 const int kRowA = 3; 42 const int kColA = 5; 43 Matrix A(kRowA, kColA); 44 A.setOnes(); 45 46 const int kRowB = 5; 47 const int kColB = 7; 48 Matrix B(kRowB, kColB); 49 B.setOnes(); 50 51 for (int row_stride_c = kRowA; row_stride_c < 3 * kRowA; ++row_stride_c) { 52 for (int col_stride_c = kColB; col_stride_c < 3 * kColB; ++col_stride_c) { 53 Matrix C(row_stride_c, col_stride_c); 54 C.setOnes(); 55 56 Matrix C_plus = C; 57 Matrix C_minus = C; 58 Matrix C_assign = C; 59 60 Matrix C_plus_ref = C; 61 Matrix C_minus_ref = C; 62 Matrix C_assign_ref = C; 63 for (int start_row_c = 0; start_row_c + kRowA < row_stride_c; ++start_row_c) { 64 for (int start_col_c = 0; start_col_c + kColB < col_stride_c; ++start_col_c) { 65 C_plus_ref.block(start_row_c, start_col_c, kRowA, kColB) += 66 A * B; 67 68 MatrixMatrixMultiply<kRowA, kColA, kRowB, kColB, 1>( 69 A.data(), kRowA, kColA, 70 B.data(), kRowB, kColB, 71 C_plus.data(), start_row_c, start_col_c, row_stride_c, col_stride_c); 72 73 EXPECT_NEAR((C_plus_ref - C_plus).norm(), 0.0, kTolerance) 74 << "C += A * B \n" 75 << "row_stride_c : " << row_stride_c << "\n" 76 << "col_stride_c : " << col_stride_c << "\n" 77 << "start_row_c : " << start_row_c << "\n" 78 << "start_col_c : " << start_col_c << "\n" 79 << "Cref : \n" << C_plus_ref << "\n" 80 << "C: \n" << C_plus; 81 82 83 C_minus_ref.block(start_row_c, start_col_c, kRowA, kColB) -= 84 A * B; 85 86 MatrixMatrixMultiply<kRowA, kColA, kRowB, kColB, -1>( 87 A.data(), kRowA, kColA, 88 B.data(), kRowB, kColB, 89 C_minus.data(), start_row_c, start_col_c, row_stride_c, col_stride_c); 90 91 EXPECT_NEAR((C_minus_ref - C_minus).norm(), 0.0, kTolerance) 92 << "C -= A * B \n" 93 << "row_stride_c : " << row_stride_c << "\n" 94 << "col_stride_c : " << col_stride_c << "\n" 95 << "start_row_c : " << start_row_c << "\n" 96 << "start_col_c : " << start_col_c << "\n" 97 << "Cref : \n" << C_minus_ref << "\n" 98 << "C: \n" << C_minus; 99 100 C_assign_ref.block(start_row_c, start_col_c, kRowA, kColB) = 101 A * B; 102 103 MatrixMatrixMultiply<kRowA, kColA, kRowB, kColB, 0>( 104 A.data(), kRowA, kColA, 105 B.data(), kRowB, kColB, 106 C_assign.data(), start_row_c, start_col_c, row_stride_c, col_stride_c); 107 108 EXPECT_NEAR((C_assign_ref - C_assign).norm(), 0.0, kTolerance) 109 << "C = A * B \n" 110 << "row_stride_c : " << row_stride_c << "\n" 111 << "col_stride_c : " << col_stride_c << "\n" 112 << "start_row_c : " << start_row_c << "\n" 113 << "start_col_c : " << start_col_c << "\n" 114 << "Cref : \n" << C_assign_ref << "\n" 115 << "C: \n" << C_assign; 116 } 117 } 118 } 119 } 120 } 121 122 TEST(BLAS, MatrixTransposeMatrixMultiply) { 123 const double kTolerance = 1e-16; 124 const int kRowA = 5; 125 const int kColA = 3; 126 Matrix A(kRowA, kColA); 127 A.setOnes(); 128 129 const int kRowB = 5; 130 const int kColB = 7; 131 Matrix B(kRowB, kColB); 132 B.setOnes(); 133 134 for (int row_stride_c = kColA; row_stride_c < 3 * kColA; ++row_stride_c) { 135 for (int col_stride_c = kColB; col_stride_c < 3 * kColB; ++col_stride_c) { 136 Matrix C(row_stride_c, col_stride_c); 137 C.setOnes(); 138 139 Matrix C_plus = C; 140 Matrix C_minus = C; 141 Matrix C_assign = C; 142 143 Matrix C_plus_ref = C; 144 Matrix C_minus_ref = C; 145 Matrix C_assign_ref = C; 146 for (int start_row_c = 0; start_row_c + kColA < row_stride_c; ++start_row_c) { 147 for (int start_col_c = 0; start_col_c + kColB < col_stride_c; ++start_col_c) { 148 C_plus_ref.block(start_row_c, start_col_c, kColA, kColB) += 149 A.transpose() * B; 150 151 MatrixTransposeMatrixMultiply<kRowA, kColA, kRowB, kColB, 1>( 152 A.data(), kRowA, kColA, 153 B.data(), kRowB, kColB, 154 C_plus.data(), start_row_c, start_col_c, row_stride_c, col_stride_c); 155 156 EXPECT_NEAR((C_plus_ref - C_plus).norm(), 0.0, kTolerance) 157 << "C += A' * B \n" 158 << "row_stride_c : " << row_stride_c << "\n" 159 << "col_stride_c : " << col_stride_c << "\n" 160 << "start_row_c : " << start_row_c << "\n" 161 << "start_col_c : " << start_col_c << "\n" 162 << "Cref : \n" << C_plus_ref << "\n" 163 << "C: \n" << C_plus; 164 165 C_minus_ref.block(start_row_c, start_col_c, kColA, kColB) -= 166 A.transpose() * B; 167 168 MatrixTransposeMatrixMultiply<kRowA, kColA, kRowB, kColB, -1>( 169 A.data(), kRowA, kColA, 170 B.data(), kRowB, kColB, 171 C_minus.data(), start_row_c, start_col_c, row_stride_c, col_stride_c); 172 173 EXPECT_NEAR((C_minus_ref - C_minus).norm(), 0.0, kTolerance) 174 << "C -= A' * B \n" 175 << "row_stride_c : " << row_stride_c << "\n" 176 << "col_stride_c : " << col_stride_c << "\n" 177 << "start_row_c : " << start_row_c << "\n" 178 << "start_col_c : " << start_col_c << "\n" 179 << "Cref : \n" << C_minus_ref << "\n" 180 << "C: \n" << C_minus; 181 182 C_assign_ref.block(start_row_c, start_col_c, kColA, kColB) = 183 A.transpose() * B; 184 185 MatrixTransposeMatrixMultiply<kRowA, kColA, kRowB, kColB, 0>( 186 A.data(), kRowA, kColA, 187 B.data(), kRowB, kColB, 188 C_assign.data(), start_row_c, start_col_c, row_stride_c, col_stride_c); 189 190 EXPECT_NEAR((C_assign_ref - C_assign).norm(), 0.0, kTolerance) 191 << "C = A' * B \n" 192 << "row_stride_c : " << row_stride_c << "\n" 193 << "col_stride_c : " << col_stride_c << "\n" 194 << "start_row_c : " << start_row_c << "\n" 195 << "start_col_c : " << start_col_c << "\n" 196 << "Cref : \n" << C_assign_ref << "\n" 197 << "C: \n" << C_assign; 198 } 199 } 200 } 201 } 202 } 203 204 TEST(BLAS, MatrixVectorMultiply) { 205 const double kTolerance = 1e-16; 206 const int kRowA = 5; 207 const int kColA = 3; 208 Matrix A(kRowA, kColA); 209 A.setOnes(); 210 211 Vector b(kColA); 212 b.setOnes(); 213 214 Vector c(kRowA); 215 c.setOnes(); 216 217 Vector c_plus = c; 218 Vector c_minus = c; 219 Vector c_assign = c; 220 221 Vector c_plus_ref = c; 222 Vector c_minus_ref = c; 223 Vector c_assign_ref = c; 224 225 c_plus_ref += A * b; 226 MatrixVectorMultiply<kRowA, kColA, 1>(A.data(), kRowA, kColA, 227 b.data(), 228 c_plus.data()); 229 EXPECT_NEAR((c_plus_ref - c_plus).norm(), 0.0, kTolerance) 230 << "c += A * b \n" 231 << "c_ref : \n" << c_plus_ref << "\n" 232 << "c: \n" << c_plus; 233 234 c_minus_ref -= A * b; 235 MatrixVectorMultiply<kRowA, kColA, -1>(A.data(), kRowA, kColA, 236 b.data(), 237 c_minus.data()); 238 EXPECT_NEAR((c_minus_ref - c_minus).norm(), 0.0, kTolerance) 239 << "c += A * b \n" 240 << "c_ref : \n" << c_minus_ref << "\n" 241 << "c: \n" << c_minus; 242 243 c_assign_ref = A * b; 244 MatrixVectorMultiply<kRowA, kColA, 0>(A.data(), kRowA, kColA, 245 b.data(), 246 c_assign.data()); 247 EXPECT_NEAR((c_assign_ref - c_assign).norm(), 0.0, kTolerance) 248 << "c += A * b \n" 249 << "c_ref : \n" << c_assign_ref << "\n" 250 << "c: \n" << c_assign; 251 } 252 253 TEST(BLAS, MatrixTransposeVectorMultiply) { 254 const double kTolerance = 1e-16; 255 const int kRowA = 5; 256 const int kColA = 3; 257 Matrix A(kRowA, kColA); 258 A.setRandom(); 259 260 Vector b(kRowA); 261 b.setRandom(); 262 263 Vector c(kColA); 264 c.setOnes(); 265 266 Vector c_plus = c; 267 Vector c_minus = c; 268 Vector c_assign = c; 269 270 Vector c_plus_ref = c; 271 Vector c_minus_ref = c; 272 Vector c_assign_ref = c; 273 274 c_plus_ref += A.transpose() * b; 275 MatrixTransposeVectorMultiply<kRowA, kColA, 1>(A.data(), kRowA, kColA, 276 b.data(), 277 c_plus.data()); 278 EXPECT_NEAR((c_plus_ref - c_plus).norm(), 0.0, kTolerance) 279 << "c += A' * b \n" 280 << "c_ref : \n" << c_plus_ref << "\n" 281 << "c: \n" << c_plus; 282 283 c_minus_ref -= A.transpose() * b; 284 MatrixTransposeVectorMultiply<kRowA, kColA, -1>(A.data(), kRowA, kColA, 285 b.data(), 286 c_minus.data()); 287 EXPECT_NEAR((c_minus_ref - c_minus).norm(), 0.0, kTolerance) 288 << "c += A' * b \n" 289 << "c_ref : \n" << c_minus_ref << "\n" 290 << "c: \n" << c_minus; 291 292 c_assign_ref = A.transpose() * b; 293 MatrixTransposeVectorMultiply<kRowA, kColA, 0>(A.data(), kRowA, kColA, 294 b.data(), 295 c_assign.data()); 296 EXPECT_NEAR((c_assign_ref - c_assign).norm(), 0.0, kTolerance) 297 << "c += A' * b \n" 298 << "c_ref : \n" << c_assign_ref << "\n" 299 << "c: \n" << c_assign; 300 } 301 302 } // namespace internal 303 } // namespace ceres 304
