1// Ceres Solver - A fast non-linear least squares minimizer
2// Copyright 2013 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#ifndef CERES_INTERNAL_PRECONDITIONER_H_
32#define CERES_INTERNAL_PRECONDITIONER_H_
33
34#include <vector>
35#include "ceres/casts.h"
36#include "ceres/compressed_row_sparse_matrix.h"
37#include "ceres/linear_operator.h"
38#include "ceres/sparse_matrix.h"
39
40namespace ceres {
41namespace internal {
42
43class BlockSparseMatrix;
44class SparseMatrix;
45
46class Preconditioner : public LinearOperator {
47 public:
48  struct Options {
49    Options()
50        : type(JACOBI),
51          sparse_linear_algebra_library_type(SUITE_SPARSE),
52          num_threads(1),
53          row_block_size(Eigen::Dynamic),
54          e_block_size(Eigen::Dynamic),
55          f_block_size(Eigen::Dynamic) {
56    }
57
58    PreconditionerType type;
59
60    SparseLinearAlgebraLibraryType sparse_linear_algebra_library_type;
61
62    // If possible, how many threads the preconditioner can use.
63    int num_threads;
64
65    // Hints about the order in which the parameter blocks should be
66    // eliminated by the linear solver.
67    //
68    // For example if elimination_groups is a vector of size k, then
69    // the linear solver is informed that it should eliminate the
70    // parameter blocks 0 ... elimination_groups[0] - 1 first, and
71    // then elimination_groups[0] ... elimination_groups[1] - 1 and so
72    // on. Within each elimination group, the linear solver is free to
73    // choose how the parameter blocks are ordered. Different linear
74    // solvers have differing requirements on elimination_groups.
75    //
76    // The most common use is for Schur type solvers, where there
77    // should be at least two elimination groups and the first
78    // elimination group must form an independent set in the normal
79    // equations. The first elimination group corresponds to the
80    // num_eliminate_blocks in the Schur type solvers.
81    vector<int> elimination_groups;
82
83    // If the block sizes in a BlockSparseMatrix are fixed, then in
84    // some cases the Schur complement based solvers can detect and
85    // specialize on them.
86    //
87    // It is expected that these parameters are set programmatically
88    // rather than manually.
89    //
90    // Please see schur_complement_solver.h and schur_eliminator.h for
91    // more details.
92    int row_block_size;
93    int e_block_size;
94    int f_block_size;
95  };
96
97  virtual ~Preconditioner();
98
99  // Update the numerical value of the preconditioner for the linear
100  // system:
101  //
102  //  |   A   | x = |b|
103  //  |diag(D)|     |0|
104  //
105  // for some vector b. It is important that the matrix A have the
106  // same block structure as the one used to construct this object.
107  //
108  // D can be NULL, in which case its interpreted as a diagonal matrix
109  // of size zero.
110  virtual bool Update(const LinearOperator& A, const double* D) = 0;
111
112  // LinearOperator interface. Since the operator is symmetric,
113  // LeftMultiply and num_cols are just calls to RightMultiply and
114  // num_rows respectively. Update() must be called before
115  // RightMultiply can be called.
116  virtual void RightMultiply(const double* x, double* y) const = 0;
117  virtual void LeftMultiply(const double* x, double* y) const {
118    return RightMultiply(x, y);
119  }
120
121  virtual int num_rows() const = 0;
122  virtual int num_cols() const {
123    return num_rows();
124  }
125};
126
127// This templated subclass of Preconditioner serves as a base class for
128// other preconditioners that depend on the particular matrix layout of
129// the underlying linear operator.
130template <typename MatrixType>
131class TypedPreconditioner : public Preconditioner {
132 public:
133  virtual ~TypedPreconditioner() {}
134  virtual bool Update(const LinearOperator& A, const double* D) {
135    return UpdateImpl(*down_cast<const MatrixType*>(&A), D);
136  }
137
138 private:
139  virtual bool UpdateImpl(const MatrixType& A, const double* D) = 0;
140};
141
142// Preconditioners that depend on acccess to the low level structure
143// of a SparseMatrix.
144typedef TypedPreconditioner<SparseMatrix>              SparseMatrixPreconditioner;               // NOLINT
145typedef TypedPreconditioner<BlockSparseMatrix>         BlockSparseMatrixPreconditioner;          // NOLINT
146typedef TypedPreconditioner<CompressedRowSparseMatrix> CompressedRowSparseMatrixPreconditioner;  // NOLINT
147
148// Wrap a SparseMatrix object as a preconditioner.
149class SparseMatrixPreconditionerWrapper : public SparseMatrixPreconditioner {
150 public:
151  // Wrapper does NOT take ownership of the matrix pointer.
152  explicit SparseMatrixPreconditionerWrapper(const SparseMatrix* matrix);
153  virtual ~SparseMatrixPreconditionerWrapper();
154
155  // Preconditioner interface
156  virtual void RightMultiply(const double* x, double* y) const;
157  virtual int num_rows() const;
158
159 private:
160  virtual bool UpdateImpl(const SparseMatrix& A, const double* D);
161  const SparseMatrix* matrix_;
162};
163
164}  // namespace internal
165}  // namespace ceres
166
167#endif  // CERES_INTERNAL_PRECONDITIONER_H_
168