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 <cstddef>
32#include "ceres/block_sparse_matrix.h"
33#include "ceres/block_structure.h"
34#include "ceres/casts.h"
35#include "ceres/internal/scoped_ptr.h"
36#include "ceres/linear_least_squares_problems.h"
37#include "ceres/linear_solver.h"
38#include "ceres/schur_complement_solver.h"
39#include "ceres/triplet_sparse_matrix.h"
40#include "ceres/types.h"
41#include "glog/logging.h"
42#include "gtest/gtest.h"
43
44namespace ceres {
45namespace internal {
46
47class SchurComplementSolverTest : public ::testing::Test {
48 protected:
49  void SetUpFromProblemId(int problem_id) {
50    scoped_ptr<LinearLeastSquaresProblem> problem(
51        CreateLinearLeastSquaresProblemFromId(problem_id));
52
53    CHECK_NOTNULL(problem.get());
54    A.reset(down_cast<BlockSparseMatrix*>(problem->A.release()));
55    b.reset(problem->b.release());
56    D.reset(problem->D.release());
57
58    num_cols = A->num_cols();
59    num_rows = A->num_rows();
60    num_eliminate_blocks = problem->num_eliminate_blocks;
61
62    x.reset(new double[num_cols]);
63    sol.reset(new double[num_cols]);
64    sol_d.reset(new double[num_cols]);
65
66    LinearSolver::Options options;
67    options.type = DENSE_QR;
68
69    scoped_ptr<LinearSolver> qr(LinearSolver::Create(options));
70
71    TripletSparseMatrix triplet_A(A->num_rows(),
72                                  A->num_cols(),
73                                  A->num_nonzeros());
74    A->ToTripletSparseMatrix(&triplet_A);
75
76    // Gold standard solutions using dense QR factorization.
77    DenseSparseMatrix dense_A(triplet_A);
78    qr->Solve(&dense_A, b.get(), LinearSolver::PerSolveOptions(), sol.get());
79
80    // Gold standard solution with appended diagonal.
81    LinearSolver::PerSolveOptions per_solve_options;
82    per_solve_options.D = D.get();
83    qr->Solve(&dense_A, b.get(), per_solve_options, sol_d.get());
84  }
85
86  void ComputeAndCompareSolutions(
87      int problem_id,
88      bool regularization,
89      ceres::LinearSolverType linear_solver_type,
90      ceres::DenseLinearAlgebraLibraryType dense_linear_algebra_library_type,
91      ceres::SparseLinearAlgebraLibraryType sparse_linear_algebra_library_type,
92      bool use_postordering) {
93    SetUpFromProblemId(problem_id);
94    LinearSolver::Options options;
95    options.elimination_groups.push_back(num_eliminate_blocks);
96    options.elimination_groups.push_back(
97        A->block_structure()->cols.size() - num_eliminate_blocks);
98    options.type = linear_solver_type;
99    options.dense_linear_algebra_library_type =
100        dense_linear_algebra_library_type;
101    options.sparse_linear_algebra_library_type =
102        sparse_linear_algebra_library_type;
103    options.use_postordering = use_postordering;
104
105    scoped_ptr<LinearSolver> solver(LinearSolver::Create(options));
106
107    LinearSolver::PerSolveOptions per_solve_options;
108    LinearSolver::Summary summary;
109    if (regularization) {
110      per_solve_options.D = D.get();
111    }
112
113    summary = solver->Solve(A.get(), b.get(), per_solve_options, x.get());
114
115    if (regularization) {
116      for (int i = 0; i < num_cols; ++i) {
117        ASSERT_NEAR(sol_d.get()[i], x[i], 1e-10);
118      }
119    } else {
120      for (int i = 0; i < num_cols; ++i) {
121        ASSERT_NEAR(sol.get()[i], x[i], 1e-10);
122      }
123    }
124  }
125
126  int num_rows;
127  int num_cols;
128  int num_eliminate_blocks;
129
130  scoped_ptr<BlockSparseMatrix> A;
131  scoped_array<double> b;
132  scoped_array<double> x;
133  scoped_array<double> D;
134  scoped_array<double> sol;
135  scoped_array<double> sol_d;
136};
137
138TEST_F(SchurComplementSolverTest, EigenBasedDenseSchurWithSmallProblem) {
139  ComputeAndCompareSolutions(2, false, DENSE_SCHUR, EIGEN, SUITE_SPARSE, true);
140  ComputeAndCompareSolutions(2, true, DENSE_SCHUR, EIGEN, SUITE_SPARSE, true);
141}
142
143TEST_F(SchurComplementSolverTest, EigenBasedDenseSchurWithLargeProblem) {
144  ComputeAndCompareSolutions(3, false, DENSE_SCHUR, EIGEN, SUITE_SPARSE, true);
145  ComputeAndCompareSolutions(3, true, DENSE_SCHUR, EIGEN, SUITE_SPARSE, true);
146}
147
148#ifndef CERES_NO_LAPACK
149TEST_F(SchurComplementSolverTest, LAPACKBasedDenseSchurWithSmallProblem) {
150  ComputeAndCompareSolutions(2, false, DENSE_SCHUR, LAPACK, SUITE_SPARSE, true);
151  ComputeAndCompareSolutions(2, true, DENSE_SCHUR, LAPACK, SUITE_SPARSE, true);
152}
153
154TEST_F(SchurComplementSolverTest, LAPACKBasedDenseSchurWithLargeProblem) {
155  ComputeAndCompareSolutions(3, false, DENSE_SCHUR, LAPACK, SUITE_SPARSE, true);
156  ComputeAndCompareSolutions(3, true, DENSE_SCHUR, LAPACK, SUITE_SPARSE, true);
157}
158#endif
159
160#ifndef CERES_NO_SUITESPARSE
161TEST_F(SchurComplementSolverTest,
162       SparseSchurWithSuiteSparseSmallProblemNoPostOrdering) {
163  ComputeAndCompareSolutions(
164      2, false, SPARSE_SCHUR, EIGEN, SUITE_SPARSE, false);
165  ComputeAndCompareSolutions(2, true, SPARSE_SCHUR, EIGEN, SUITE_SPARSE, false);
166}
167
168TEST_F(SchurComplementSolverTest,
169       SparseSchurWithSuiteSparseSmallProblemPostOrdering) {
170  ComputeAndCompareSolutions(2, false, SPARSE_SCHUR, EIGEN, SUITE_SPARSE, true);
171  ComputeAndCompareSolutions(2, true, SPARSE_SCHUR, EIGEN, SUITE_SPARSE, true);
172}
173
174TEST_F(SchurComplementSolverTest,
175       SparseSchurWithSuiteSparseLargeProblemNoPostOrdering) {
176  ComputeAndCompareSolutions(
177      3, false, SPARSE_SCHUR, EIGEN, SUITE_SPARSE, false);
178  ComputeAndCompareSolutions(3, true, SPARSE_SCHUR, EIGEN, SUITE_SPARSE, false);
179}
180
181TEST_F(SchurComplementSolverTest,
182       SparseSchurWithSuiteSparseLargeProblemPostOrdering) {
183  ComputeAndCompareSolutions(3, false, SPARSE_SCHUR, EIGEN, SUITE_SPARSE, true);
184  ComputeAndCompareSolutions(3, true, SPARSE_SCHUR, EIGEN, SUITE_SPARSE, true);
185}
186#endif  // CERES_NO_SUITESPARSE
187
188#ifndef CERES_NO_CXSPARSE
189TEST_F(SchurComplementSolverTest,
190       SparseSchurWithCXSparseSmallProblem) {
191  ComputeAndCompareSolutions(2, false, SPARSE_SCHUR, EIGEN, CX_SPARSE, true);
192  ComputeAndCompareSolutions(2, true, SPARSE_SCHUR, EIGEN, CX_SPARSE, true);
193}
194
195TEST_F(SchurComplementSolverTest,
196       SparseSchurWithCXSparseLargeProblem) {
197  ComputeAndCompareSolutions(3, false, SPARSE_SCHUR, EIGEN, CX_SPARSE, true);
198  ComputeAndCompareSolutions(3, true, SPARSE_SCHUR, EIGEN, CX_SPARSE, true);
199}
200#endif  // CERES_NO_CXSPARSE
201
202#ifdef CERES_USE_EIGEN_SPARSE
203TEST_F(SchurComplementSolverTest,
204       SparseSchurWithEigenSparseSmallProblem) {
205  ComputeAndCompareSolutions(2, false, SPARSE_SCHUR, EIGEN, EIGEN_SPARSE, true);
206  ComputeAndCompareSolutions(2, true, SPARSE_SCHUR, EIGEN, EIGEN_SPARSE, true);
207}
208
209TEST_F(SchurComplementSolverTest,
210       SparseSchurWithEigenSparseLargeProblem) {
211  ComputeAndCompareSolutions(3, false, SPARSE_SCHUR, EIGEN, EIGEN_SPARSE, true);
212  ComputeAndCompareSolutions(3, true, SPARSE_SCHUR, EIGEN, EIGEN_SPARSE, true);
213}
214#endif  // CERES_USE_EIGEN_SPARSE
215
216}  // namespace internal
217}  // namespace ceres
218