1 // Ceres Solver - A fast non-linear least squares minimizer
2 // Copyright 2010, 2011, 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@google.com (Sameer Agarwal)
30
31 #include "ceres/casts.h"
32 #include "ceres/compressed_row_sparse_matrix.h"
33 #include "ceres/internal/scoped_ptr.h"
34 #include "ceres/linear_least_squares_problems.h"
35 #include "ceres/linear_solver.h"
36 #include "ceres/triplet_sparse_matrix.h"
37 #include "ceres/types.h"
38 #include "glog/logging.h"
39 #include "gtest/gtest.h"
40
41
42 namespace ceres {
43 namespace internal {
44
45 class UnsymmetricLinearSolverTest : public ::testing::Test {
46 protected :
SetUp()47 virtual void SetUp() {
48 scoped_ptr<LinearLeastSquaresProblem> problem(
49 CreateLinearLeastSquaresProblemFromId(0));
50
51 CHECK_NOTNULL(problem.get());
52 A_.reset(down_cast<TripletSparseMatrix*>(problem->A.release()));
53 b_.reset(problem->b.release());
54 D_.reset(problem->D.release());
55 sol_unregularized_.reset(problem->x.release());
56 sol_regularized_.reset(problem->x_D.release());
57 }
58
TestSolver(const LinearSolver::Options & options)59 void TestSolver(const LinearSolver::Options& options) {
60
61
62 LinearSolver::PerSolveOptions per_solve_options;
63 LinearSolver::Summary unregularized_solve_summary;
64 LinearSolver::Summary regularized_solve_summary;
65 Vector x_unregularized(A_->num_cols());
66 Vector x_regularized(A_->num_cols());
67
68 scoped_ptr<SparseMatrix> transformed_A;
69
70 if (options.type == DENSE_QR ||
71 options.type == DENSE_NORMAL_CHOLESKY) {
72 transformed_A.reset(new DenseSparseMatrix(*A_));
73 } else if (options.type == SPARSE_NORMAL_CHOLESKY) {
74 CompressedRowSparseMatrix* crsm = new CompressedRowSparseMatrix(*A_);
75 // Add row/column blocks structure.
76 for (int i = 0; i < A_->num_rows(); ++i) {
77 crsm->mutable_row_blocks()->push_back(1);
78 }
79
80 for (int i = 0; i < A_->num_cols(); ++i) {
81 crsm->mutable_col_blocks()->push_back(1);
82 }
83 transformed_A.reset(crsm);
84 } else {
85 LOG(FATAL) << "Unknown linear solver : " << options.type;
86 }
87
88 // Unregularized
89 scoped_ptr<LinearSolver> solver(LinearSolver::Create(options));
90 unregularized_solve_summary =
91 solver->Solve(transformed_A.get(),
92 b_.get(),
93 per_solve_options,
94 x_unregularized.data());
95
96 // Sparsity structure is changing, reset the solver.
97 solver.reset(LinearSolver::Create(options));
98 // Regularized solution
99 per_solve_options.D = D_.get();
100 regularized_solve_summary =
101 solver->Solve(transformed_A.get(),
102 b_.get(),
103 per_solve_options,
104 x_regularized.data());
105
106 EXPECT_EQ(unregularized_solve_summary.termination_type,
107 LINEAR_SOLVER_SUCCESS);
108
109 for (int i = 0; i < A_->num_cols(); ++i) {
110 EXPECT_NEAR(sol_unregularized_[i], x_unregularized[i], 1e-8)
111 << "\nExpected: "
112 << ConstVectorRef(sol_unregularized_.get(), A_->num_cols()).transpose()
113 << "\nActual: " << x_unregularized.transpose();
114 }
115
116 EXPECT_EQ(regularized_solve_summary.termination_type,
117 LINEAR_SOLVER_SUCCESS);
118 for (int i = 0; i < A_->num_cols(); ++i) {
119 EXPECT_NEAR(sol_regularized_[i], x_regularized[i], 1e-8)
120 << "\nExpected: "
121 << ConstVectorRef(sol_regularized_.get(), A_->num_cols()).transpose()
122 << "\nActual: " << x_regularized.transpose();
123 }
124 }
125
126 scoped_ptr<TripletSparseMatrix> A_;
127 scoped_array<double> b_;
128 scoped_array<double> D_;
129 scoped_array<double> sol_unregularized_;
130 scoped_array<double> sol_regularized_;
131 };
132
TEST_F(UnsymmetricLinearSolverTest,EigenDenseQR)133 TEST_F(UnsymmetricLinearSolverTest, EigenDenseQR) {
134 LinearSolver::Options options;
135 options.type = DENSE_QR;
136 options.dense_linear_algebra_library_type = EIGEN;
137 TestSolver(options);
138 }
139
TEST_F(UnsymmetricLinearSolverTest,EigenDenseNormalCholesky)140 TEST_F(UnsymmetricLinearSolverTest, EigenDenseNormalCholesky) {
141 LinearSolver::Options options;
142 options.dense_linear_algebra_library_type = EIGEN;
143 options.type = DENSE_NORMAL_CHOLESKY;
144 TestSolver(options);
145 }
146
147 #ifndef CERES_NO_LAPACK
TEST_F(UnsymmetricLinearSolverTest,LAPACKDenseQR)148 TEST_F(UnsymmetricLinearSolverTest, LAPACKDenseQR) {
149 LinearSolver::Options options;
150 options.type = DENSE_QR;
151 options.dense_linear_algebra_library_type = LAPACK;
152 TestSolver(options);
153 }
154
TEST_F(UnsymmetricLinearSolverTest,LAPACKDenseNormalCholesky)155 TEST_F(UnsymmetricLinearSolverTest, LAPACKDenseNormalCholesky) {
156 LinearSolver::Options options;
157 options.dense_linear_algebra_library_type = LAPACK;
158 options.type = DENSE_NORMAL_CHOLESKY;
159 TestSolver(options);
160 }
161 #endif
162
163 #ifndef CERES_NO_SUITESPARSE
TEST_F(UnsymmetricLinearSolverTest,SparseNormalCholeskyUsingSuiteSparsePreOrdering)164 TEST_F(UnsymmetricLinearSolverTest,
165 SparseNormalCholeskyUsingSuiteSparsePreOrdering) {
166 LinearSolver::Options options;
167 options.sparse_linear_algebra_library_type = SUITE_SPARSE;
168 options.type = SPARSE_NORMAL_CHOLESKY;
169 options.use_postordering = false;
170 TestSolver(options);
171 }
172
TEST_F(UnsymmetricLinearSolverTest,SparseNormalCholeskyUsingSuiteSparsePostOrdering)173 TEST_F(UnsymmetricLinearSolverTest,
174 SparseNormalCholeskyUsingSuiteSparsePostOrdering) {
175 LinearSolver::Options options;
176 options.sparse_linear_algebra_library_type = SUITE_SPARSE;
177 options.type = SPARSE_NORMAL_CHOLESKY;
178 options.use_postordering = true;
179 TestSolver(options);
180 }
181
TEST_F(UnsymmetricLinearSolverTest,SparseNormalCholeskyUsingSuiteSparseDynamicSparsity)182 TEST_F(UnsymmetricLinearSolverTest,
183 SparseNormalCholeskyUsingSuiteSparseDynamicSparsity) {
184 LinearSolver::Options options;
185 options.sparse_linear_algebra_library_type = SUITE_SPARSE;
186 options.type = SPARSE_NORMAL_CHOLESKY;
187 options.dynamic_sparsity = true;
188 TestSolver(options);
189 }
190 #endif
191
192 #ifndef CERES_NO_CXSPARSE
TEST_F(UnsymmetricLinearSolverTest,SparseNormalCholeskyUsingCXSparsePreOrdering)193 TEST_F(UnsymmetricLinearSolverTest,
194 SparseNormalCholeskyUsingCXSparsePreOrdering) {
195 LinearSolver::Options options;
196 options.sparse_linear_algebra_library_type = CX_SPARSE;
197 options.type = SPARSE_NORMAL_CHOLESKY;
198 options.use_postordering = false;
199 TestSolver(options);
200 }
201
TEST_F(UnsymmetricLinearSolverTest,SparseNormalCholeskyUsingCXSparsePostOrdering)202 TEST_F(UnsymmetricLinearSolverTest,
203 SparseNormalCholeskyUsingCXSparsePostOrdering) {
204 LinearSolver::Options options;
205 options.sparse_linear_algebra_library_type = CX_SPARSE;
206 options.type = SPARSE_NORMAL_CHOLESKY;
207 options.use_postordering = true;
208 TestSolver(options);
209 }
210
TEST_F(UnsymmetricLinearSolverTest,SparseNormalCholeskyUsingCXSparseDynamicSparsity)211 TEST_F(UnsymmetricLinearSolverTest,
212 SparseNormalCholeskyUsingCXSparseDynamicSparsity) {
213 LinearSolver::Options options;
214 options.sparse_linear_algebra_library_type = CX_SPARSE;
215 options.type = SPARSE_NORMAL_CHOLESKY;
216 options.dynamic_sparsity = true;
217 TestSolver(options);
218 }
219 #endif
220
221 #ifdef CERES_USE_EIGEN_SPARSE
TEST_F(UnsymmetricLinearSolverTest,SparseNormalCholeskyUsingEigenPreOrdering)222 TEST_F(UnsymmetricLinearSolverTest,
223 SparseNormalCholeskyUsingEigenPreOrdering) {
224 LinearSolver::Options options;
225 options.sparse_linear_algebra_library_type = EIGEN_SPARSE;
226 options.type = SPARSE_NORMAL_CHOLESKY;
227 options.use_postordering = false;
228 TestSolver(options);
229 }
230
TEST_F(UnsymmetricLinearSolverTest,SparseNormalCholeskyUsingEigenPostOrdering)231 TEST_F(UnsymmetricLinearSolverTest,
232 SparseNormalCholeskyUsingEigenPostOrdering) {
233 LinearSolver::Options options;
234 options.sparse_linear_algebra_library_type = EIGEN_SPARSE;
235 options.type = SPARSE_NORMAL_CHOLESKY;
236 options.use_postordering = true;
237 TestSolver(options);
238 }
239
TEST_F(UnsymmetricLinearSolverTest,SparseNormalCholeskyUsingEigenDynamicSparsity)240 TEST_F(UnsymmetricLinearSolverTest,
241 SparseNormalCholeskyUsingEigenDynamicSparsity) {
242 LinearSolver::Options options;
243 options.sparse_linear_algebra_library_type = EIGEN_SPARSE;
244 options.type = SPARSE_NORMAL_CHOLESKY;
245 options.dynamic_sparsity = true;
246 TestSolver(options);
247 }
248
249 #endif // CERES_USE_EIGEN_SPARSE
250
251 } // namespace internal
252 } // namespace ceres
253