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#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.
38extern "C" void dpotrf_(char* uplo,
39                       int* n,
40                       double* a,
41                       int* lda,
42                       int* info);
43
44extern "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
53extern "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
66namespace ceres {
67namespace internal {
68
69LinearSolverTerminationType LAPACK::SolveInPlaceUsingCholesky(
70    int num_rows,
71    const double* in_lhs,
72    double* rhs_and_solution,
73    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
115int 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
147LinearSolverTerminationType 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    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
194