1 // Ceres Solver - A fast non-linear least squares minimizer
2 // Copyright 2015 Google Inc. All rights reserved.
3 // http://ceres-solver.org/
5 // Redistribution and use in source and binary forms, with or without
6 // modification, are permitted provided that the following conditions are met:
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.
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.
29 // Author: keir@google.com (Keir Mierle)
31 #include "ceres/small_blas.h"
34 #include "gtest/gtest.h"
35 #include "ceres/internal/eigen.h"
40 const double kTolerance = 3.0 * std::numeric_limits<double>::epsilon();
42 TEST(BLAS, MatrixMatrixMultiply) {
45 Matrix A(kRowA, kColA);
50 Matrix B(kRowB, kColB);
53 for (int row_stride_c = kRowA; row_stride_c < 3 * kRowA; ++row_stride_c) {
54 for (int col_stride_c = kColB; col_stride_c < 3 * kColB; ++col_stride_c) {
55 Matrix C(row_stride_c, col_stride_c);
62 Matrix C_plus_ref = C;
63 Matrix C_minus_ref = C;
64 Matrix C_assign_ref = C;
65 for (int start_row_c = 0; start_row_c + kRowA < row_stride_c; ++start_row_c) {
66 for (int start_col_c = 0; start_col_c + kColB < col_stride_c; ++start_col_c) {
67 C_plus_ref.block(start_row_c, start_col_c, kRowA, kColB) +=
70 MatrixMatrixMultiply<kRowA, kColA, kRowB, kColB, 1>(
71 A.data(), kRowA, kColA,
72 B.data(), kRowB, kColB,
73 C_plus.data(), start_row_c, start_col_c, row_stride_c, col_stride_c);
75 EXPECT_NEAR((C_plus_ref - C_plus).norm(), 0.0, kTolerance)
77 << "row_stride_c : " << row_stride_c << "\n"
78 << "col_stride_c : " << col_stride_c << "\n"
79 << "start_row_c : " << start_row_c << "\n"
80 << "start_col_c : " << start_col_c << "\n"
81 << "Cref : \n" << C_plus_ref << "\n"
85 C_minus_ref.block(start_row_c, start_col_c, kRowA, kColB) -=
88 MatrixMatrixMultiply<kRowA, kColA, kRowB, kColB, -1>(
89 A.data(), kRowA, kColA,
90 B.data(), kRowB, kColB,
91 C_minus.data(), start_row_c, start_col_c, row_stride_c, col_stride_c);
93 EXPECT_NEAR((C_minus_ref - C_minus).norm(), 0.0, kTolerance)
95 << "row_stride_c : " << row_stride_c << "\n"
96 << "col_stride_c : " << col_stride_c << "\n"
97 << "start_row_c : " << start_row_c << "\n"
98 << "start_col_c : " << start_col_c << "\n"
99 << "Cref : \n" << C_minus_ref << "\n"
100 << "C: \n" << C_minus;
102 C_assign_ref.block(start_row_c, start_col_c, kRowA, kColB) =
105 MatrixMatrixMultiply<kRowA, kColA, kRowB, kColB, 0>(
106 A.data(), kRowA, kColA,
107 B.data(), kRowB, kColB,
108 C_assign.data(), start_row_c, start_col_c, row_stride_c, col_stride_c);
110 EXPECT_NEAR((C_assign_ref - C_assign).norm(), 0.0, kTolerance)
112 << "row_stride_c : " << row_stride_c << "\n"
113 << "col_stride_c : " << col_stride_c << "\n"
114 << "start_row_c : " << start_row_c << "\n"
115 << "start_col_c : " << start_col_c << "\n"
116 << "Cref : \n" << C_assign_ref << "\n"
117 << "C: \n" << C_assign;
124 TEST(BLAS, MatrixTransposeMatrixMultiply) {
127 Matrix A(kRowA, kColA);
132 Matrix B(kRowB, kColB);
135 for (int row_stride_c = kColA; row_stride_c < 3 * kColA; ++row_stride_c) {
136 for (int col_stride_c = kColB; col_stride_c < 3 * kColB; ++col_stride_c) {
137 Matrix C(row_stride_c, col_stride_c);
144 Matrix C_plus_ref = C;
145 Matrix C_minus_ref = C;
146 Matrix C_assign_ref = C;
147 for (int start_row_c = 0; start_row_c + kColA < row_stride_c; ++start_row_c) {
148 for (int start_col_c = 0; start_col_c + kColB < col_stride_c; ++start_col_c) {
149 C_plus_ref.block(start_row_c, start_col_c, kColA, kColB) +=
152 MatrixTransposeMatrixMultiply<kRowA, kColA, kRowB, kColB, 1>(
153 A.data(), kRowA, kColA,
154 B.data(), kRowB, kColB,
155 C_plus.data(), start_row_c, start_col_c, row_stride_c, col_stride_c);
157 EXPECT_NEAR((C_plus_ref - C_plus).norm(), 0.0, kTolerance)
159 << "row_stride_c : " << row_stride_c << "\n"
160 << "col_stride_c : " << col_stride_c << "\n"
161 << "start_row_c : " << start_row_c << "\n"
162 << "start_col_c : " << start_col_c << "\n"
163 << "Cref : \n" << C_plus_ref << "\n"
164 << "C: \n" << C_plus;
166 C_minus_ref.block(start_row_c, start_col_c, kColA, kColB) -=
169 MatrixTransposeMatrixMultiply<kRowA, kColA, kRowB, kColB, -1>(
170 A.data(), kRowA, kColA,
171 B.data(), kRowB, kColB,
172 C_minus.data(), start_row_c, start_col_c, row_stride_c, col_stride_c);
174 EXPECT_NEAR((C_minus_ref - C_minus).norm(), 0.0, kTolerance)
176 << "row_stride_c : " << row_stride_c << "\n"
177 << "col_stride_c : " << col_stride_c << "\n"
178 << "start_row_c : " << start_row_c << "\n"
179 << "start_col_c : " << start_col_c << "\n"
180 << "Cref : \n" << C_minus_ref << "\n"
181 << "C: \n" << C_minus;
183 C_assign_ref.block(start_row_c, start_col_c, kColA, kColB) =
186 MatrixTransposeMatrixMultiply<kRowA, kColA, kRowB, kColB, 0>(
187 A.data(), kRowA, kColA,
188 B.data(), kRowB, kColB,
189 C_assign.data(), start_row_c, start_col_c, row_stride_c, col_stride_c);
191 EXPECT_NEAR((C_assign_ref - C_assign).norm(), 0.0, kTolerance)
193 << "row_stride_c : " << row_stride_c << "\n"
194 << "col_stride_c : " << col_stride_c << "\n"
195 << "start_row_c : " << start_row_c << "\n"
196 << "start_col_c : " << start_col_c << "\n"
197 << "Cref : \n" << C_assign_ref << "\n"
198 << "C: \n" << C_assign;
205 TEST(BLAS, MatrixVectorMultiply) {
208 Matrix A(kRowA, kColA);
221 Vector c_plus_ref = c;
222 Vector c_minus_ref = c;
223 Vector c_assign_ref = c;
226 MatrixVectorMultiply<kRowA, kColA, 1>(A.data(), kRowA, kColA,
229 EXPECT_NEAR((c_plus_ref - c_plus).norm(), 0.0, kTolerance)
231 << "c_ref : \n" << c_plus_ref << "\n"
232 << "c: \n" << c_plus;
234 c_minus_ref -= A * b;
235 MatrixVectorMultiply<kRowA, kColA, -1>(A.data(), kRowA, kColA,
238 EXPECT_NEAR((c_minus_ref - c_minus).norm(), 0.0, kTolerance)
240 << "c_ref : \n" << c_minus_ref << "\n"
241 << "c: \n" << c_minus;
243 c_assign_ref = A * b;
244 MatrixVectorMultiply<kRowA, kColA, 0>(A.data(), kRowA, kColA,
247 EXPECT_NEAR((c_assign_ref - c_assign).norm(), 0.0, kTolerance)
249 << "c_ref : \n" << c_assign_ref << "\n"
250 << "c: \n" << c_assign;
253 TEST(BLAS, MatrixTransposeVectorMultiply) {
256 Matrix A(kRowA, kColA);
269 Vector c_plus_ref = c;
270 Vector c_minus_ref = c;
271 Vector c_assign_ref = c;
273 c_plus_ref += A.transpose() * b;
274 MatrixTransposeVectorMultiply<kRowA, kColA, 1>(A.data(), kRowA, kColA,
277 EXPECT_NEAR((c_plus_ref - c_plus).norm(), 0.0, kTolerance)
279 << "c_ref : \n" << c_plus_ref << "\n"
280 << "c: \n" << c_plus;
282 c_minus_ref -= A.transpose() * b;
283 MatrixTransposeVectorMultiply<kRowA, kColA, -1>(A.data(), kRowA, kColA,
286 EXPECT_NEAR((c_minus_ref - c_minus).norm(), 0.0, kTolerance)
288 << "c_ref : \n" << c_minus_ref << "\n"
289 << "c: \n" << c_minus;
291 c_assign_ref = A.transpose() * b;
292 MatrixTransposeVectorMultiply<kRowA, kColA, 0>(A.data(), kRowA, kColA,
295 EXPECT_NEAR((c_assign_ref - c_assign).norm(), 0.0, kTolerance)
297 << "c_ref : \n" << c_assign_ref << "\n"
298 << "c: \n" << c_assign;
301 } // namespace internal