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@google.com (Keir Mierle)
30
31#include "ceres/small_blas.h"
32
33#include "gtest/gtest.h"
34#include "ceres/internal/eigen.h"
35
36namespace ceres {
37namespace internal {
38
39TEST(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
122TEST(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
204TEST(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
253TEST(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