internal/ceres/lapack.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: sameeragarwal@google.com (Sameer Agarwal)
  30
  31 #include "ceres/lapack.h"
  32
  33 #include "ceres/internal/port.h"
  34 #include "ceres/linear_solver.h"
  35 #include "glog/logging.h"
  36
  37 // C interface to the LAPACK Cholesky factorization and triangular solve.
  38 extern "C" void dpotrf_(char* uplo,
  39                        int* n,
  40                        double* a,
  41                        int* lda,
  42                        int* info);
  43
  44 extern "C" void dpotrs_(char* uplo,
  45                         int* n,
  46                         int* nrhs,
  47                         double* a,
  48                         int* lda,
  49                         double* b,
  50                         int* ldb,
  51                         int* info);
  52
  53 extern "C" void dgels_(char* uplo,
  54                        int* m,
  55                        int* n,
  56                        int* nrhs,
  57                        double* a,
  58                        int* lda,
  59                        double* b,
  60                        int* ldb,
  61                        double* work,
  62                        int* lwork,
  63                        int* info);
  64
  65
  66 namespace ceres {
  67 namespace internal {
  68
  69 LinearSolverTerminationType LAPACK::SolveInPlaceUsingCholesky(
  70     int num_rows,
  71     const double* in_lhs,
  72     double* rhs_and_solution,
  73     std::string* message) {
  74 #ifdef CERES_NO_LAPACK
  75   LOG(FATAL) << "Ceres was built without a BLAS library.";
  76   return LINEAR_SOLVER_FATAL_ERROR;
  77 #else
  78   char uplo = 'L';
  79   int n = num_rows;
  80   int info = 0;
  81   int nrhs = 1;
  82   double* lhs = const_cast<double*>(in_lhs);
  83
  84   dpotrf_(&uplo, &n, lhs, &n, &info);
  85   if (info < 0) {
  86     LOG(FATAL) << "Congratulations, you found a bug in Ceres."
  87                << "Please report it."
  88                << "LAPACK::dpotrf fatal error."
  89                << "Argument: " << -info << " is invalid.";
  90     return LINEAR_SOLVER_FATAL_ERROR;
  91   }
  92
  93   if (info > 0) {
  94     *message =
  95         StringPrintf(
  96             "LAPACK::dpotrf numerical failure. "
  97              "The leading minor of order %d is not positive definite.", info);
  98     return LINEAR_SOLVER_FAILURE;
  99   }
 100
 101   dpotrs_(&uplo, &n, &nrhs, lhs, &n, rhs_and_solution, &n, &info);
 102   if (info < 0) {
 103     LOG(FATAL) << "Congratulations, you found a bug in Ceres."
 104                << "Please report it."
 105                << "LAPACK::dpotrs fatal error."
 106                << "Argument: " << -info << " is invalid.";
 107     return LINEAR_SOLVER_FATAL_ERROR;
 108   }
 109
 110   *message = "Success";
 111   return LINEAR_SOLVER_SUCCESS;
 112 #endif
 113 }
 114
 115 int LAPACK::EstimateWorkSizeForQR(int num_rows, int num_cols) {
 116 #ifdef CERES_NO_LAPACK
 117   LOG(FATAL) << "Ceres was built without a LAPACK library.";
 118   return -1;
 119 #else
 120   char trans = 'N';
 121   int nrhs = 1;
 122   int lwork = -1;
 123   double work;
 124   int info = 0;
 125   dgels_(&trans,
 126          &num_rows,
 127          &num_cols,
 128          &nrhs,
 129          NULL,
 130          &num_rows,
 131          NULL,
 132          &num_rows,
 133          &work,
 134          &lwork,
 135          &info);
 136
 137   if (info < 0) {
 138     LOG(FATAL) << "Congratulations, you found a bug in Ceres."
 139                << "Please report it."
 140                << "LAPACK::dgels fatal error."
 141                << "Argument: " << -info << " is invalid.";
 142   }
 143   return static_cast<int>(work);
 144 #endif
 145 }
 146
 147 LinearSolverTerminationType LAPACK::SolveInPlaceUsingQR(
 148     int num_rows,
 149     int num_cols,
 150     const double* in_lhs,
 151     int work_size,
 152     double* work,
 153     double* rhs_and_solution,
 154     std::string* message) {
 155 #ifdef CERES_NO_LAPACK
 156   LOG(FATAL) << "Ceres was built without a LAPACK library.";
 157   return LINEAR_SOLVER_FATAL_ERROR;
 158 #else
 159   char trans = 'N';
 160   int m = num_rows;
 161   int n = num_cols;
 162   int nrhs = 1;
 163   int lda = num_rows;
 164   int ldb = num_rows;
 165   int info = 0;
 166   double* lhs = const_cast<double*>(in_lhs);
 167
 168   dgels_(&trans,
 169          &m,
 170          &n,
 171          &nrhs,
 172          lhs,
 173          &lda,
 174          rhs_and_solution,
 175          &ldb,
 176          work,
 177          &work_size,
 178          &info);
 179
 180   if (info < 0) {
 181     LOG(FATAL) << "Congratulations, you found a bug in Ceres."
 182                << "Please report it."
 183                << "LAPACK::dgels fatal error."
 184                << "Argument: " << -info << " is invalid.";
 185   }
 186
 187   *message = "Success.";
 188   return LINEAR_SOLVER_SUCCESS;
 189 #endif
 190 }
 191
 192 }  // namespace internal
 193 }  // namespace ceres