internal/ceres/cgnr_solver.cc

   1 // Ceres Solver - A fast non-linear least squares minimizer
   2 // Copyright 2015 Google Inc. All rights reserved.
   3 // http://ceres-solver.org/
   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: keir@google.com (Keir Mierle)
  30
  31 #include "ceres/cgnr_solver.h"
  32
  33 #include "ceres/block_jacobi_preconditioner.h"
  34 #include "ceres/cgnr_linear_operator.h"
  35 #include "ceres/conjugate_gradients_solver.h"
  36 #include "ceres/internal/eigen.h"
  37 #include "ceres/linear_solver.h"
  38 #include "ceres/wall_time.h"
  39 #include "glog/logging.h"
  40
  41 namespace ceres {
  42 namespace internal {
  43
  44 CgnrSolver::CgnrSolver(const LinearSolver::Options& options)
  45   : options_(options),
  46     preconditioner_(NULL) {
  47   if (options_.preconditioner_type != JACOBI &&
  48       options_.preconditioner_type != IDENTITY) {
  49     LOG(FATAL) << "CGNR only supports IDENTITY and JACOBI preconditioners.";
  50   }
  51 }
  52
  53 LinearSolver::Summary CgnrSolver::SolveImpl(
  54     BlockSparseMatrix* A,
  55     const double* b,
  56     const LinearSolver::PerSolveOptions& per_solve_options,
  57     double* x) {
  58   EventLogger event_logger("CgnrSolver::Solve");
  59
  60   // Form z = Atb.
  61   Vector z(A->num_cols());
  62   z.setZero();
  63   A->LeftMultiply(b, z.data());
  64
  65   // Precondition if necessary.
  66   LinearSolver::PerSolveOptions cg_per_solve_options = per_solve_options;
  67   if (options_.preconditioner_type == JACOBI) {
  68     if (preconditioner_.get() == NULL) {
  69       preconditioner_.reset(new BlockJacobiPreconditioner(*A));
  70     }
  71     preconditioner_->Update(*A, per_solve_options.D);
  72     cg_per_solve_options.preconditioner = preconditioner_.get();
  73   }
  74
  75   // Solve (AtA + DtD)x = z (= Atb).
  76   VectorRef(x, A->num_cols()).setZero();
  77   CgnrLinearOperator lhs(*A, per_solve_options.D);
  78   event_logger.AddEvent("Setup");
  79
  80   ConjugateGradientsSolver conjugate_gradient_solver(options_);
  81   LinearSolver::Summary summary =
  82       conjugate_gradient_solver.Solve(&lhs, z.data(), cg_per_solve_options, x);
  83   event_logger.AddEvent("Solve");
  84   return summary;
  85 }
  86
  87 }  // namespace internal
  88 }  // namespace ceres