Home | History | Annotate | Download | only in ceres
      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/small_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