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