schur_complement_solver_test.cc revision 1d2624a10e2c559f8ba9ef89eaa30832c0a83a96 
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::SparseLinearAlgebraLibraryType sparse_linear_algebra_library,
91      bool use_postordering) {
92    SetUpFromProblemId(problem_id);
93    LinearSolver::Options options;
94    options.elimination_groups.push_back(num_eliminate_blocks);
95    options.elimination_groups.push_back(
96        A->block_structure()->cols.size() - num_eliminate_blocks);
97    options.type = linear_solver_type;
98    options.sparse_linear_algebra_library = sparse_linear_algebra_library;
99    options.use_postordering = use_postordering;
100
101    scoped_ptr<LinearSolver> solver(LinearSolver::Create(options));
102
103    LinearSolver::PerSolveOptions per_solve_options;
104    LinearSolver::Summary summary;
105    if (regularization) {
106      per_solve_options.D = D.get();
107    }
108
109    summary = solver->Solve(A.get(), b.get(), per_solve_options, x.get());
110
111    if (regularization) {
112      for (int i = 0; i < num_cols; ++i) {
113        ASSERT_NEAR(sol_d.get()[i], x[i], 1e-10);
114      }
115    } else {
116      for (int i = 0; i < num_cols; ++i) {
117        ASSERT_NEAR(sol.get()[i], x[i], 1e-10);
118      }
119    }
120  }
121
122  int num_rows;
123  int num_cols;
124  int num_eliminate_blocks;
125
126  scoped_ptr<BlockSparseMatrix> A;
127  scoped_array<double> b;
128  scoped_array<double> x;
129  scoped_array<double> D;
130  scoped_array<double> sol;
131  scoped_array<double> sol_d;
132};
133
134TEST_F(SchurComplementSolverTest, DenseSchurWithSmallProblem) {
135  ComputeAndCompareSolutions(2, false, DENSE_SCHUR, SUITE_SPARSE, true);
136  ComputeAndCompareSolutions(2, true, DENSE_SCHUR, SUITE_SPARSE, true);
137}
138
139TEST_F(SchurComplementSolverTest, DenseSchurWithLargeProblem) {
140  ComputeAndCompareSolutions(3, false, DENSE_SCHUR, SUITE_SPARSE, true);
141  ComputeAndCompareSolutions(3, true, DENSE_SCHUR, SUITE_SPARSE, true);
142}
143
144#ifndef CERES_NO_SUITESPARSE
145TEST_F(SchurComplementSolverTest,
146       SparseSchurWithSuiteSparseSmallProblemNoPostOrdering) {
147  ComputeAndCompareSolutions(2, false, SPARSE_SCHUR, SUITE_SPARSE, false);
148  ComputeAndCompareSolutions(2, true, SPARSE_SCHUR, SUITE_SPARSE, false);
149}
150
151TEST_F(SchurComplementSolverTest,
152       SparseSchurWithSuiteSparseSmallProblemPostOrdering) {
153  ComputeAndCompareSolutions(2, false, SPARSE_SCHUR, SUITE_SPARSE, true);
154  ComputeAndCompareSolutions(2, true, SPARSE_SCHUR, SUITE_SPARSE, true);
155}
156
157TEST_F(SchurComplementSolverTest,
158       SparseSchurWithSuiteSparseLargeProblemNoPostOrdering) {
159  ComputeAndCompareSolutions(3, false, SPARSE_SCHUR, SUITE_SPARSE, false);
160  ComputeAndCompareSolutions(3, true, SPARSE_SCHUR, SUITE_SPARSE, false);
161}
162
163TEST_F(SchurComplementSolverTest,
164       SparseSchurWithSuiteSparseLargeProblemPostOrdering) {
165  ComputeAndCompareSolutions(3, false, SPARSE_SCHUR, SUITE_SPARSE, true);
166  ComputeAndCompareSolutions(3, true, SPARSE_SCHUR, SUITE_SPARSE, true);
167}
168#endif  // CERES_NO_SUITESPARSE
169
170#ifndef CERES_NO_CXSPARSE
171TEST_F(SchurComplementSolverTest,
172       SparseSchurWithSuiteSparseSmallProblem) {
173  ComputeAndCompareSolutions(2, false, SPARSE_SCHUR, SUITE_SPARSE, true);
174  ComputeAndCompareSolutions(2, true, SPARSE_SCHUR, SUITE_SPARSE, true);
175}
176
177TEST_F(SchurComplementSolverTest,
178       SparseSchurWithSuiteSparseLargeProblem) {
179  ComputeAndCompareSolutions(3, false, SPARSE_SCHUR, SUITE_SPARSE, true);
180  ComputeAndCompareSolutions(3, true, SPARSE_SCHUR, SUITE_SPARSE, true);
181}
182#endif  // CERES_NO_CXSPARSE
183
184}  // namespace internal
185}  // namespace ceres
186