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 
41 namespace ceres {
42 namespace internal {
43 
44 class BlockSparseMatrix;
45 class SparseMatrix;
46 
47 class Preconditioner : public LinearOperator {
48  public:
49   struct Options {
OptionsOptions50     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;
LeftMultiply(const double * x,double * y)127   virtual void LeftMultiply(const double* x, double* y) const {
128     return RightMultiply(x, y);
129   }
130 
131   virtual int num_rows() const = 0;
num_cols()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.
140 template <typename MatrixType>
141 class TypedPreconditioner : public Preconditioner {
142  public:
~TypedPreconditioner()143   virtual ~TypedPreconditioner() {}
Update(const LinearOperator & A,const double * D)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.
154 typedef TypedPreconditioner<SparseMatrix>              SparseMatrixPreconditioner;               // NOLINT
155 typedef TypedPreconditioner<BlockSparseMatrix>         BlockSparseMatrixPreconditioner;          // NOLINT
156 typedef TypedPreconditioner<CompressedRowSparseMatrix> CompressedRowSparseMatrixPreconditioner;  // NOLINT
157 
158 // Wrap a SparseMatrix object as a preconditioner.
159 class 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