Home | History | Annotate | Download | only in ceres 
      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 (at) 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 
     40 namespace ceres {
     41 namespace internal {
     42 
     43 class BlockSparseMatrix;
     44 class SparseMatrix;
     45 
     46 class 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.
    130 template <typename MatrixType>
    131 class 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.
    144 typedef TypedPreconditioner<SparseMatrix>              SparseMatrixPreconditioner;               // NOLINT
    145 typedef TypedPreconditioner<BlockSparseMatrix>         BlockSparseMatrixPreconditioner;          // NOLINT
    146 typedef TypedPreconditioner<CompressedRowSparseMatrix> CompressedRowSparseMatrixPreconditioner;  // NOLINT
    147 
    148 // Wrap a SparseMatrix object as a preconditioner.
    149 class 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