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