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
42namespace ceres {
43namespace internal {
44
45class UnsymmetricLinearSolverTest : public ::testing::Test {
46 protected :
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
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
133TEST_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
140TEST_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
148TEST_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
155TEST_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
164TEST_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
173TEST_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
182TEST_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
193TEST_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
202TEST_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
211TEST_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
222TEST_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
231TEST_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
240TEST_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