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 #include "ceres/lapack.h"
     32 #include "glog/logging.h"
     33 
     34 // C interface to the LAPACK Cholesky factorization and triangular solve.
     35 extern "C" void dpotrf_(char* uplo,
     36                        int* n,
     37                        double* a,
     38                        int* lda,
     39                        int* info);
     40 
     41 extern "C" void dpotrs_(char* uplo,
     42                         int* n,
     43                         int* nrhs,
     44                         double* a,
     45                         int* lda,
     46                         double* b,
     47                         int* ldb,
     48                         int* info);
     49 
     50 extern "C" void dgels_(char* uplo,
     51                        int* m,
     52                        int* n,
     53                        int* nrhs,
     54                        double* a,
     55                        int* lda,
     56                        double* b,
     57                        int* ldb,
     58                        double* work,
     59                        int* lwork,
     60                        int* info);
     61 
     62 
     63 namespace ceres {
     64 namespace internal {
     65 
     66 int LAPACK::SolveInPlaceUsingCholesky(int num_rows,
     67                                       const double* in_lhs,
     68                                       double* rhs_and_solution) {
     69 #ifdef CERES_NO_LAPACK
     70   LOG(FATAL) << "Ceres was built without a BLAS library.";
     71   return -1;
     72 #else
     73   char uplo = 'L';
     74   int n = num_rows;
     75   int info = 0;
     76   int nrhs = 1;
     77   double* lhs = const_cast<double*>(in_lhs);
     78 
     79   dpotrf_(&uplo, &n, lhs, &n, &info);
     80   if (info != 0) {
     81     LOG(INFO) << "Cholesky factorization (dpotrf) failed: " << info;
     82     return info;
     83   }
     84 
     85   dpotrs_(&uplo, &n, &nrhs, lhs, &n, rhs_and_solution, &n, &info);
     86   if (info != 0) {
     87     LOG(INFO) << "Triangular solve (dpotrs) failed: " << info;
     88   }
     89 
     90   return info;
     91 #endif
     92 };
     93 
     94 int LAPACK::EstimateWorkSizeForQR(int num_rows, int num_cols) {
     95 #ifdef CERES_NO_LAPACK
     96   LOG(FATAL) << "Ceres was built without a LAPACK library.";
     97   return -1;
     98 #else
     99   char trans = 'N';
    100   int nrhs = 1;
    101   int lwork = -1;
    102   double work;
    103   int info = 0;
    104   dgels_(&trans,
    105          &num_rows,
    106          &num_cols,
    107          &nrhs,
    108          NULL,
    109          &num_rows,
    110          NULL,
    111          &num_rows,
    112          &work,
    113          &lwork,
    114          &info);
    115 
    116   CHECK_EQ(info, 0);
    117   return work;
    118 #endif
    119 }
    120 
    121 int LAPACK::SolveUsingQR(int num_rows,
    122                          int num_cols,
    123                          const double* in_lhs,
    124                          int work_size,
    125                          double* work,
    126                          double* rhs_and_solution) {
    127 #ifdef CERES_NO_LAPACK
    128   LOG(FATAL) << "Ceres was built without a LAPACK library.";
    129   return -1;
    130 #else
    131   char trans = 'N';
    132   int m = num_rows;
    133   int n = num_cols;
    134   int nrhs = 1;
    135   int lda = num_rows;
    136   int ldb = num_rows;
    137   int info = 0;
    138   double* lhs = const_cast<double*>(in_lhs);
    139 
    140   dgels_(&trans,
    141          &m,
    142          &n,
    143          &nrhs,
    144          lhs,
    145          &lda,
    146          rhs_and_solution,
    147          &ldb,
    148          work,
    149          &work_size,
    150          &info);
    151 
    152   return info;
    153 #endif
    154 }
    155 
    156 }  // namespace internal
    157 }  // namespace ceres
    158