Add DenseQR Interface
1. Add EigenDenseQR & tests.
This implementation now uses an in place decomposition,
which means that we are not allocating, deallocating
memory every call.
2. Add LAPACKDenseQR and tests.
The LAPACK implementation instead of using dgels which is a
routine which does the factorization and solve in one
call, now uses dgeqrf for factorization and then
dormqr and dtrtrs for solving. This allows us to
have a factorize and solve interface like DenseCholesky.
And opens the door to iterative refinement and mixed
precision solves.
3. The refactor also allows us to simplify the interface to
DenseSparseMatrix considerably. The internals of this
class were complicated because we had the AppendDiagonal
and RemoveDiagonal methods and we did not want to allocate
deallocate memory every call. But since we pay the cost
of the copy anyways, we can just hold that buffer
in DenseQRSolver.
4. Delete lapack.cc/h
5. The net result is that everything seems to be a bit faster.
For LAPACK we are not doing some of the scaling work that
dgels was doing. For Eigen I think it maybe the inplace
decomposition.
Benchmark Time CPU Time Old Time New CPU Old CPU New
----------------------------------------------------------------------------------------------------------------------------------------------------------
BM_DenseSolver<ceres::EIGEN, ceres::DENSE_QR>/1/1 -0.1154 -0.1159 692 612 691 611
BM_DenseSolver<ceres::EIGEN, ceres::DENSE_QR>/2/1 -0.1601 -0.1553 717 603 712 601
BM_DenseSolver<ceres::EIGEN, ceres::DENSE_QR>/3/1 -0.1673 -0.1575 733 610 724 610
BM_DenseSolver<ceres::EIGEN, ceres::DENSE_QR>/6/2 -0.1008 -0.1003 886 797 884 796
BM_DenseSolver<ceres::EIGEN, ceres::DENSE_QR>/10/3 -0.1489 -0.1514 1283 1092 1281 1087
BM_DenseSolver<ceres::EIGEN, ceres::DENSE_QR>/12/4 -0.1040 -0.1104 1556 1394 1553 1381
BM_DenseSolver<ceres::EIGEN, ceres::DENSE_QR>/20/5 -0.0007 -0.0097 1911 1910 1908 1890
BM_DenseSolver<ceres::EIGEN, ceres::DENSE_QR>/40/5 -0.1033 -0.1022 2981 2673 2957 2655
BM_DenseSolver<ceres::EIGEN, ceres::DENSE_QR>/100/10 -0.0147 +0.0015 9275 9138 9026 9040
BM_DenseSolver<ceres::EIGEN, ceres::DENSE_QR>/200/10 -0.1408 -0.1284 15093 12968 14778 12880
BM_DenseSolver<ceres::EIGEN, ceres::DENSE_QR>/200/20 -0.0310 -0.0355 38973 37765 38837 37460
BM_DenseSolver<ceres::LAPACK, ceres::DENSE_QR>/1/1 -0.1228 -0.1256 736 646 731 640
BM_DenseSolver<ceres::LAPACK, ceres::DENSE_QR>/2/1 -0.1401 -0.1396 740 636 735 633
BM_DenseSolver<ceres::LAPACK, ceres::DENSE_QR>/3/1 -0.1731 -0.1695 744 615 738 613
BM_DenseSolver<ceres::LAPACK, ceres::DENSE_QR>/6/2 -0.1399 -0.1408 1121 965 1113 956
BM_DenseSolver<ceres::LAPACK, ceres::DENSE_QR>/10/3 -0.1110 -0.1145 1571 1397 1560 1382
BM_DenseSolver<ceres::LAPACK, ceres::DENSE_QR>/12/4 -0.1411 -0.1417 2006 1722 1993 1710
BM_DenseSolver<ceres::LAPACK, ceres::DENSE_QR>/20/5 -0.1740 -0.1729 2741 2264 2724 2253
BM_DenseSolver<ceres::LAPACK, ceres::DENSE_QR>/40/5 -0.0966 -0.1123 3462 3128 3425 3040
BM_DenseSolver<ceres::LAPACK, ceres::DENSE_QR>/100/10 -0.0387 -0.0998 10365 9964 10339 9307
BM_DenseSolver<ceres::LAPACK, ceres::DENSE_QR>/200/10 -0.2044 -0.2049 16031 12754 15998 12720
BM_DenseSolver<ceres::LAPACK, ceres::DENSE_QR>/200/20 -0.2391 -0.2386 35777 27223 35716 27193
Change-Id: I782f0d7664efe1435eebda92ddf47a0fe66c9c72
diff --git a/examples/nist.cc b/examples/nist.cc
index f12dbf3..3944586 100644
--- a/examples/nist.cc
+++ b/examples/nist.cc
@@ -99,6 +99,9 @@
"dense_qr",
"Options are: sparse_cholesky, dense_qr, dense_normal_cholesky "
"and cgnr");
+DEFINE_string(dense_linear_algebra_library,
+ "eigen",
+ "Options are: eigen and lapack.");
DEFINE_string(preconditioner, "jacobi", "Options are: identity, jacobi");
DEFINE_string(line_search,
"wolfe",
@@ -466,6 +469,9 @@
&options->minimizer_type));
CHECK(ceres::StringToLinearSolverType(CERES_GET_FLAG(FLAGS_linear_solver),
&options->linear_solver_type));
+ CHECK(StringToDenseLinearAlgebraLibraryType(
+ CERES_GET_FLAG(FLAGS_dense_linear_algebra_library),
+ &options->dense_linear_algebra_library_type));
CHECK(ceres::StringToPreconditionerType(CERES_GET_FLAG(FLAGS_preconditioner),
&options->preconditioner_type));
CHECK(ceres::StringToTrustRegionStrategyType(
diff --git a/internal/ceres/CMakeLists.txt b/internal/ceres/CMakeLists.txt
index 5d9e0ab..957d476 100644
--- a/internal/ceres/CMakeLists.txt
+++ b/internal/ceres/CMakeLists.txt
@@ -74,6 +74,7 @@
cxsparse.cc
dense_cholesky.cc
dense_normal_cholesky_solver.cc
+ dense_qr.cc
dense_qr_solver.cc
dense_sparse_matrix.cc
detect_structure.cc
@@ -97,7 +98,6 @@
iterative_refiner.cc
iterative_schur_complement_solver.cc
levenberg_marquardt_strategy.cc
- lapack.cc
line_search.cc
line_search_direction.cc
line_search_minimizer.cc
@@ -455,6 +455,7 @@
ceres_test(cubic_interpolation)
ceres_test(dense_linear_solver)
ceres_test(dense_cholesky)
+ ceres_test(dense_qr)
ceres_test(dense_sparse_matrix)
ceres_test(detect_structure)
ceres_test(dogleg_strategy)
diff --git a/internal/ceres/dense_jacobian_writer.h b/internal/ceres/dense_jacobian_writer.h
index 1d7187c..efecf35 100644
--- a/internal/ceres/dense_jacobian_writer.h
+++ b/internal/ceres/dense_jacobian_writer.h
@@ -59,8 +59,8 @@
}
SparseMatrix* CreateJacobian() const {
- return new DenseSparseMatrix(
- program_->NumResiduals(), program_->NumEffectiveParameters(), true);
+ return new DenseSparseMatrix(program_->NumResiduals(),
+ program_->NumEffectiveParameters());
}
void Write(int residual_id,
@@ -86,10 +86,10 @@
ConstMatrixRef parameter_jacobian(
jacobians[j], num_residuals, parameter_block_size);
- dense_jacobian->mutable_matrix().block(residual_offset,
- parameter_block->delta_offset(),
- num_residuals,
- parameter_block_size) =
+ dense_jacobian->mutable_matrix()->block(residual_offset,
+ parameter_block->delta_offset(),
+ num_residuals,
+ parameter_block_size) =
parameter_jacobian;
}
}
diff --git a/internal/ceres/dense_linear_solver_benchmark.cc b/internal/ceres/dense_linear_solver_benchmark.cc
index fa0d034..b8f1e00 100644
--- a/internal/ceres/dense_linear_solver_benchmark.cc
+++ b/internal/ceres/dense_linear_solver_benchmark.cc
@@ -43,10 +43,8 @@
static void BM_DenseSolver(benchmark::State& state) {
const int num_rows = state.range(0);
const int num_cols = state.range(1);
- constexpr bool kReserveDiagonal = true;
-
- DenseSparseMatrix jacobian(num_rows, num_cols, kReserveDiagonal);
- jacobian.mutable_matrix() = Eigen::MatrixXd::Random(num_rows, num_cols);
+ DenseSparseMatrix jacobian(num_rows, num_cols);
+ *jacobian.mutable_matrix() = Eigen::MatrixXd::Random(num_rows, num_cols);
Eigen::VectorXd rhs = Eigen::VectorXd::Random(num_rows, 1);
Eigen::VectorXd solution(num_cols);
diff --git a/internal/ceres/dense_linear_solver_test.cc b/internal/ceres/dense_linear_solver_test.cc
index 419ec3e..4dff5af 100644
--- a/internal/ceres/dense_linear_solver_test.cc
+++ b/internal/ceres/dense_linear_solver_test.cc
@@ -89,24 +89,22 @@
solver->Solve(&lhs, rhs.data(), per_solve_options, solution.data());
EXPECT_EQ(summary.termination_type, LINEAR_SOLVER_SUCCESS);
- // If solving for the regularized solution, add the diagonal to the
- // matrix. This makes subsequent computations simpler.
- if (testing::get<2>(param)) {
- lhs.AppendDiagonal(problem->D.get());
- };
+ Vector normal_rhs = lhs.matrix().transpose() * rhs.head(num_rows);
+ Matrix normal_lhs = lhs.matrix().transpose() * lhs.matrix();
- Vector tmp = Vector::Zero(num_rows + num_cols);
- lhs.RightMultiply(solution.data(), tmp.data());
- Vector actual_normal_rhs = Vector::Zero(num_cols);
- lhs.LeftMultiply(tmp.data(), actual_normal_rhs.data());
+ if (regularized) {
+ ConstVectorRef diagonal(problem->D.get(), num_cols);
+ normal_lhs += diagonal.array().square().matrix().asDiagonal();
+ }
- Vector expected_normal_rhs = Vector::Zero(num_cols);
- lhs.LeftMultiply(rhs.data(), expected_normal_rhs.data());
- const double residual = (expected_normal_rhs - actual_normal_rhs).norm() /
- expected_normal_rhs.norm();
+ Vector actual_normal_rhs = normal_lhs * solution;
- EXPECT_NEAR(residual, 0.0, 10 * std::numeric_limits<double>::epsilon())
- << "\nexpected: " << expected_normal_rhs.transpose()
+ const double normalized_residual =
+ (normal_rhs - actual_normal_rhs).norm() / normal_rhs.norm();
+
+ EXPECT_NEAR(
+ normalized_residual, 0.0, 10 * std::numeric_limits<double>::epsilon())
+ << "\nexpected: " << normal_rhs.transpose()
<< "\nactual: " << actual_normal_rhs.transpose();
}
diff --git a/internal/ceres/dense_normal_cholesky_solver.cc b/internal/ceres/dense_normal_cholesky_solver.cc
index 73e5bf8..b208d58 100644
--- a/internal/ceres/dense_normal_cholesky_solver.cc
+++ b/internal/ceres/dense_normal_cholesky_solver.cc
@@ -57,7 +57,6 @@
const int num_rows = A->num_rows();
const int num_cols = A->num_cols();
- ConstColMajorMatrixRef Aref = A->matrix();
Matrix lhs(num_cols, num_cols);
lhs.setZero();
@@ -68,10 +67,10 @@
// Using rankUpdate instead of GEMM, exposes the fact that its the
// same matrix being multiplied with itself and that the product is
// symmetric.
- lhs.selfadjointView<Eigen::Upper>().rankUpdate(Aref.transpose());
+ lhs.selfadjointView<Eigen::Upper>().rankUpdate(A->matrix().transpose());
// rhs = A'b
- Vector rhs = Aref.transpose() * ConstVectorRef(b, num_rows);
+ Vector rhs = A->matrix().transpose() * ConstVectorRef(b, num_rows);
if (per_solve_options.D != NULL) {
ConstVectorRef D(per_solve_options.D, num_cols);
diff --git a/internal/ceres/dense_qr.cc b/internal/ceres/dense_qr.cc
new file mode 100644
index 0000000..9c0b24e
--- /dev/null
+++ b/internal/ceres/dense_qr.cc
@@ -0,0 +1,297 @@
+// Ceres Solver - A fast non-linear least squares minimizer
+// Copyright 2022 Google Inc. All rights reserved.
+// http://ceres-solver.org/
+//
+// Redistribution and use in source and binary forms, with or without
+// modification, are permitted provided that the following conditions are met:
+//
+// * Redistributions of source code must retain the above copyright notice,
+// this list of conditions and the following disclaimer.
+// * Redistributions in binary form must reproduce the above copyright notice,
+// this list of conditions and the following disclaimer in the documentation
+// and/or other materials provided with the distribution.
+// * Neither the name of Google Inc. nor the names of its contributors may be
+// used to endorse or promote products derived from this software without
+// specific prior written permission.
+//
+// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
+// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+// POSSIBILITY OF SUCH DAMAGE.
+//
+// Author: sameeragarwal@google.com (Sameer Agarwal)
+
+#include "ceres/dense_qr.h"
+
+#include <algorithm>
+#include <memory>
+
+#ifndef CERES_NO_LAPACK
+
+// LAPACK routines for solving a linear least squares problem using QR
+// factorization. This is done in three stages:
+//
+// A * x = b
+// Q * R * x = b (dgeqrf)
+// R * x = Q' * b (dormqr)
+// x = R^{-1} * Q'* b (dtrtrs)
+
+// clang-format off
+
+// Compute the QR factorization of a.
+//
+// a is an m x n column major matrix (Denoted by "A" in the above description)
+// lda is the leading dimension of a. lda >= max(1, num_rows)
+// tau is an array of size min(m,n). It contains the scalar factors of the
+// elementary reflectors.
+// work is an array of size max(1,lwork). On exit, if info=0, work[0] contains
+// the optimal size of work.
+//
+// if lwork >= 1 it is the size of work. If lwork = -1, then a workspace query is assumed.
+// dgeqrf computes the optimal size of the work array and returns it as work[0].
+//
+// info = 0, successful exit.
+// info < 0, if info = -i, then the i^th argument had illegal value.
+extern "C" void dgeqrf_(const int* m, const int* n, double* a, const int* lda,
+ double* tau, double* work, const int* lwork, int* info);
+
+// Apply Q or Q' to b.
+//
+// b is a m times n column major matrix.
+// size = 'L' applies Q or Q' on the left, size = 'R' applies Q or Q' on the right.
+// trans = 'N', applies Q, trans = 'T', applies Q'.
+// k is the number of elementary reflectors whose product defines the matrix Q.
+// If size = 'L', m >= k >= 0 and if side = 'R', n >= k >= 0.
+// a is an lda x k column major matrix containing the reflectors as returned by dgeqrf.
+// ldb is the leading dimension of b.
+// work is an array of size max(1, lwork)
+// lwork if positive is the size of work. If lwork = -1, then a
+// workspace query is assumed.
+//
+// info = 0, successful exit.
+// info < 0, if info = -i, then the i^th argument had illegal value.
+extern "C" void dormqr_(const char* side, const char* trans, const int* m,
+ const int* n ,const int* k, double* a, const int* lda,
+ double* tau, double* b, const int* ldb, double* work,
+ const int* lwork, int* info);
+
+// Solve a triangular system of the form A * x = b
+//
+// uplo = 'U', A is upper triangular. uplo = 'L' is lower triangular.
+// trans = 'N', 'T', 'C' specifies the form - A, A^T, A^H.
+// DIAG = 'N', A is not unit triangular. 'U' is unit triangular.
+// n is the order of the matrix A.
+// nrhs number of columns of b.
+// a is a column major lda x n.
+// b is a column major matrix of ldb x nrhs
+//
+// info = 0 succesful.
+// = -i < 0 i^th argument is an illegal value.
+// = i > 0, i^th diagonal element of A is zero.
+extern "C" void dtrtrs_(const char* uplo, const char* trans, const char* diag,
+ const int* n, const int* nrhs, double* a, const int* lda,
+ double* b, const int* ldb, int* info);
+// clang-format on
+
+#endif
+
+namespace ceres {
+namespace internal {
+
+std::unique_ptr<DenseQR> DenseQR::Create(const LinearSolver::Options& options) {
+ std::unique_ptr<DenseQR> dense_qr;
+
+ switch (options.dense_linear_algebra_library_type) {
+ case EIGEN:
+ dense_qr = std::make_unique<EigenDenseQR>();
+ break;
+
+ case LAPACK:
+#ifndef CERES_NO_LAPACK
+ dense_qr = std::make_unique<LAPACKDenseQR>();
+ break;
+#else
+ LOG(FATAL) << "Ceres was compiled without support for LAPACK.";
+#endif
+
+ default:
+ LOG(FATAL) << "Unknown dense linear algebra library type : "
+ << DenseLinearAlgebraLibraryTypeToString(
+ options.dense_linear_algebra_library_type);
+ }
+ return dense_qr;
+}
+
+LinearSolverTerminationType DenseQR::FactorAndSolve(int num_rows,
+ int num_cols,
+ double* lhs,
+ const double* rhs,
+ double* solution,
+ std::string* message) {
+ LinearSolverTerminationType termination_type =
+ Factorize(num_rows, num_cols, lhs, message);
+ if (termination_type == LINEAR_SOLVER_SUCCESS) {
+ termination_type = Solve(rhs, solution, message);
+ }
+ return termination_type;
+}
+
+LinearSolverTerminationType EigenDenseQR::Factorize(int num_rows,
+ int num_cols,
+ double* lhs,
+ std::string* message) {
+ Eigen::Map<ColMajorMatrix> m(lhs, num_rows, num_cols);
+ qr_ = std::make_unique<QRType>(m);
+ *message = "Success.";
+ return LINEAR_SOLVER_SUCCESS;
+}
+
+LinearSolverTerminationType EigenDenseQR::Solve(const double* rhs,
+ double* solution,
+ std::string* message) {
+ VectorRef(solution, qr_->cols()) =
+ qr_->solve(ConstVectorRef(rhs, qr_->rows()));
+ *message = "Success.";
+ return LINEAR_SOLVER_SUCCESS;
+}
+
+#ifndef CERES_NO_LAPACK
+LinearSolverTerminationType LAPACKDenseQR::Factorize(int num_rows,
+ int num_cols,
+ double* lhs,
+ std::string* message) {
+ int lwork = -1;
+ double work_size;
+ int info = 0;
+
+ // Compute the size of the temporary workspace needed to compute the QR
+ // factorization in the dgeqrf call below.
+ dgeqrf_(&num_rows,
+ &num_cols,
+ lhs_,
+ &num_rows,
+ tau_.data(),
+ &work_size,
+ &lwork,
+ &info);
+ if (info < 0) {
+ LOG(FATAL) << "Congratulations, you found a bug in Ceres."
+ << "Please report it."
+ << "LAPACK::dgels fatal error."
+ << "Argument: " << -info << " is invalid.";
+ }
+
+ lhs_ = lhs;
+ num_rows_ = num_rows;
+ num_cols_ = num_cols;
+
+ lwork = static_cast<int>(work_size);
+
+ if (work_.size() < lwork) {
+ work_.resize(lwork);
+ }
+ if (tau_.size() < num_cols) {
+ tau_.resize(num_cols);
+ }
+
+ if (q_transpose_rhs_.size() < num_rows) {
+ q_transpose_rhs_.resize(num_rows);
+ }
+
+ // Factorize the lhs_ using the workspace that we just constructed above.
+ dgeqrf_(&num_rows,
+ &num_cols,
+ lhs_,
+ &num_rows,
+ tau_.data(),
+ work_.data(),
+ &lwork,
+ &info);
+
+ if (info < 0) {
+ LOG(FATAL) << "Congratulations, you found a bug in Ceres."
+ << "Please report it. dgeqrf fatal error."
+ << "Argument: " << -info << " is invalid.";
+ }
+
+ termination_type_ = LINEAR_SOLVER_SUCCESS;
+ *message = "Success.";
+ return termination_type_;
+}
+
+LinearSolverTerminationType LAPACKDenseQR::Solve(const double* rhs,
+ double* solution,
+ std::string* message) {
+ if (termination_type_ != LINEAR_SOLVER_SUCCESS) {
+ *message = "QR factorization failed and solve called.";
+ return termination_type_;
+ }
+
+ std::copy_n(rhs, num_rows_, q_transpose_rhs_.data());
+
+ const char side = 'L';
+ char trans = 'T';
+ const int num_c_cols = 1;
+ const int lwork = work_.size();
+ int info = 0;
+ dormqr_(&side,
+ &trans,
+ &num_rows_,
+ &num_c_cols,
+ &num_cols_,
+ lhs_,
+ &num_rows_,
+ tau_.data(),
+ q_transpose_rhs_.data(),
+ &num_rows_,
+ work_.data(),
+ &lwork,
+ &info);
+ if (info < 0) {
+ LOG(FATAL) << "Congratulations, you found a bug in Ceres."
+ << "Please report it. dormr fatal error."
+ << "Argument: " << -info << " is invalid.";
+ }
+
+ const char uplo = 'U';
+ trans = 'N';
+ const char diag = 'N';
+ dtrtrs_(&uplo,
+ &trans,
+ &diag,
+ &num_cols_,
+ &num_c_cols,
+ lhs_,
+ &num_rows_,
+ q_transpose_rhs_.data(),
+ &num_rows_,
+ &info);
+
+ if (info < 0) {
+ LOG(FATAL) << "Congratulations, you found a bug in Ceres."
+ << "Please report it. dormr fatal error."
+ << "Argument: " << -info << " is invalid.";
+ } else if (info > 0) {
+ *message =
+ "QR factorization failure. The factorization is not full rank. R has "
+ "zeros on the diagonal.";
+ termination_type_ = LINEAR_SOLVER_FAILURE;
+ } else {
+ std::copy_n(q_transpose_rhs_.data(), num_cols_, solution);
+ termination_type_ = LINEAR_SOLVER_SUCCESS;
+ }
+
+ return termination_type_;
+}
+
+#endif // CERES_NO_LAPACK
+
+} // namespace internal
+} // namespace ceres
diff --git a/internal/ceres/dense_qr.h b/internal/ceres/dense_qr.h
new file mode 100644
index 0000000..93d890b
--- /dev/null
+++ b/internal/ceres/dense_qr.h
@@ -0,0 +1,143 @@
+// Ceres Solver - A fast non-linear least squares minimizer
+// Copyright 2022 Google Inc. All rights reserved.
+// http://ceres-solver.org/
+//
+// Redistribution and use in source and binary forms, with or without
+// modification, are permitted provided that the following conditions are met:
+//
+// * Redistributions of source code must retain the above copyright notice,
+// this list of conditions and the following disclaimer.
+// * Redistributions in binary form must reproduce the above copyright notice,
+// this list of conditions and the following disclaimer in the documentation
+// and/or other materials provided with the distribution.
+// * Neither the name of Google Inc. nor the names of its contributors may be
+// used to endorse or promote products derived from this software without
+// specific prior written permission.
+//
+// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
+// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+// POSSIBILITY OF SUCH DAMAGE.
+//
+// Author: sameeragarwal@google.com (Sameer Agarwal)
+
+#ifndef CERES_INTERNAL_DENSE_QR_H_
+#define CERES_INTERNAL_DENSE_QR_H_
+
+// This include must come before any #ifndef check on Ceres compile options.
+// clang-format off
+#include "ceres/internal/port.h"
+// clang-format on
+
+#include <memory>
+#include <vector>
+
+#include "Eigen/Dense"
+#include "ceres/internal/eigen.h"
+#include "ceres/linear_solver.h"
+#include "glog/logging.h"
+
+namespace ceres {
+namespace internal {
+
+// An interface that abstracts away the internal details of various dense linear
+// algebra libraries and offers a simple API for solving dense linear systems
+// using a QR factorization.
+class CERES_EXPORT_INTERNAL DenseQR {
+ public:
+ static std::unique_ptr<DenseQR> Create(const LinearSolver::Options& options);
+
+ virtual ~DenseQR() = default;
+
+ // Computes the QR factorization of the given matrix.
+ //
+ // The input matrix lhs is assumed to be a column-major num_rows x num_cols
+ // matrix.
+ //
+ // The input matrix lhs may be modified by the implementation to store the
+ // factorization, irrespective of whether the factorization succeeds or not.
+ // As a result it is the user's responsibility to ensure that lhs is valid
+ // when Solve is called.
+ virtual LinearSolverTerminationType Factorize(int num_rows,
+ int num_cols,
+ double* lhs,
+ std::string* message) = 0;
+
+ // Computes the solution to the equation
+ //
+ // lhs * solution = rhs
+ //
+ // Calling Solve without calling Factorize is undefined behaviour. It is the
+ // user's responsibility to ensure that the input matrix lhs passed to
+ // Factorize has not been freed/modified when Solve is called.
+ virtual LinearSolverTerminationType Solve(const double* rhs,
+ double* solution,
+ std::string* message) = 0;
+
+ // Convenience method which combines a call to Factorize and Solve. Solve is
+ // only called if Factorize returns LINEAR_SOLVER_SUCCESS.
+ //
+ // The input matrix lhs may be modified by the implementation to store the
+ // factorization, irrespective of whether the method succeeds or not. It is
+ // the user's responsibility to ensure that lhs is valid if and when Solve is
+ // called again after this call.
+ LinearSolverTerminationType FactorAndSolve(int num_rows,
+ int num_cols,
+ double* lhs,
+ const double* rhs,
+ double* solution,
+ std::string* message);
+};
+
+class CERES_EXPORT_INTERNAL EigenDenseQR : public DenseQR {
+ public:
+ ~EigenDenseQR() override = default;
+
+ LinearSolverTerminationType Factorize(int num_rows,
+ int num_cols,
+ double* lhs,
+ std::string* message) override;
+ LinearSolverTerminationType Solve(const double* rhs,
+ double* solution,
+ std::string* message) override;
+
+ private:
+ using QRType = Eigen::HouseholderQR<Eigen::Ref<ColMajorMatrix>>;
+ std::unique_ptr<QRType> qr_;
+};
+
+#ifndef CERES_NO_LAPACK
+class CERES_EXPORT_INTERNAL LAPACKDenseQR : public DenseQR {
+ public:
+ ~LAPACKDenseQR() override = default;
+
+ LinearSolverTerminationType Factorize(int num_rows,
+ int num_cols,
+ double* lhs,
+ std::string* message) override;
+ LinearSolverTerminationType Solve(const double* rhs,
+ double* solution,
+ std::string* message) override;
+
+ private:
+ double* lhs_ = nullptr;
+ int num_rows_;
+ int num_cols_;
+ LinearSolverTerminationType termination_type_ = LINEAR_SOLVER_FATAL_ERROR;
+ Vector work_;
+ Vector tau_;
+ Vector q_transpose_rhs_;
+};
+#endif // CERES_NO_LAPACK
+
+} // namespace internal
+} // namespace ceres
+
+#endif // CERES_INTERNAL_DENSE_QR_H_
diff --git a/internal/ceres/dense_qr_solver.cc b/internal/ceres/dense_qr_solver.cc
index 44388f3..24cb25a 100644
--- a/internal/ceres/dense_qr_solver.cc
+++ b/internal/ceres/dense_qr_solver.cc
@@ -33,9 +33,9 @@
#include <cstddef>
#include "Eigen/Dense"
+#include "ceres/dense_qr.h"
#include "ceres/dense_sparse_matrix.h"
#include "ceres/internal/eigen.h"
-#include "ceres/lapack.h"
#include "ceres/linear_solver.h"
#include "ceres/types.h"
#include "ceres/wall_time.h"
@@ -44,124 +44,40 @@
namespace internal {
DenseQRSolver::DenseQRSolver(const LinearSolver::Options& options)
- : options_(options) {
- work_.resize(1);
-}
+ : options_(options), dense_qr_(DenseQR::Create(options)) {}
LinearSolver::Summary DenseQRSolver::SolveImpl(
DenseSparseMatrix* A,
const double* b,
const LinearSolver::PerSolveOptions& per_solve_options,
double* x) {
- if (options_.dense_linear_algebra_library_type == EIGEN) {
- return SolveUsingEigen(A, b, per_solve_options, x);
- } else {
- return SolveUsingLAPACK(A, b, per_solve_options, x);
- }
-}
-
-LinearSolver::Summary DenseQRSolver::SolveUsingLAPACK(
- DenseSparseMatrix* A,
- const double* b,
- const LinearSolver::PerSolveOptions& per_solve_options,
- double* x) {
EventLogger event_logger("DenseQRSolver::Solve");
const int num_rows = A->num_rows();
const int num_cols = A->num_cols();
+ const int num_augmented_rows =
+ num_rows + ((per_solve_options.D != nullptr) ? num_cols : 0);
- if (per_solve_options.D != NULL) {
- // Temporarily append a diagonal block to the A matrix, but undo
- // it before returning the matrix to the user.
- A->AppendDiagonal(per_solve_options.D);
+ if (lhs_.rows() != num_augmented_rows || lhs_.cols() != num_cols) {
+ lhs_.resize(num_augmented_rows, num_cols);
+ rhs_.resize(num_augmented_rows);
}
- // TODO(sameeragarwal): Since we are copying anyways, the diagonal
- // can be appended to the matrix instead of doing it on A.
- lhs_ = A->matrix();
-
- if (per_solve_options.D != NULL) {
- // Undo the modifications to the matrix A.
- A->RemoveDiagonal();
- }
-
- // rhs = [b;0] to account for the additional rows in the lhs.
- if (rhs_.rows() != lhs_.rows()) {
- rhs_.resize(lhs_.rows());
- }
- rhs_.setZero();
+ lhs_.topRows(num_rows) = A->matrix();
rhs_.head(num_rows) = ConstVectorRef(b, num_rows);
- if (work_.rows() == 1) {
- const int work_size =
- LAPACK::EstimateWorkSizeForQR(lhs_.rows(), lhs_.cols());
- VLOG(3) << "Working memory for Dense QR factorization: "
- << work_size * sizeof(double);
- work_.resize(work_size);
+ if (num_rows != num_augmented_rows) {
+ lhs_.bottomRows(num_cols) =
+ ConstVectorRef(per_solve_options.D, num_cols).asDiagonal();
+ rhs_.tail(num_cols).setZero();
}
LinearSolver::Summary summary;
+ summary.termination_type = dense_qr_->FactorAndSolve(
+ lhs_.rows(), lhs_.cols(), lhs_.data(), rhs_.data(), x, &summary.message);
summary.num_iterations = 1;
- summary.termination_type = LAPACK::SolveInPlaceUsingQR(lhs_.rows(),
- lhs_.cols(),
- lhs_.data(),
- work_.rows(),
- work_.data(),
- rhs_.data(),
- &summary.message);
- event_logger.AddEvent("Solve");
- if (summary.termination_type == LINEAR_SOLVER_SUCCESS) {
- VectorRef(x, num_cols) = rhs_.head(num_cols);
- }
-
- event_logger.AddEvent("TearDown");
- return summary;
-}
-
-LinearSolver::Summary DenseQRSolver::SolveUsingEigen(
- DenseSparseMatrix* A,
- const double* b,
- const LinearSolver::PerSolveOptions& per_solve_options,
- double* x) {
- EventLogger event_logger("DenseQRSolver::Solve");
-
- const int num_rows = A->num_rows();
- const int num_cols = A->num_cols();
-
- if (per_solve_options.D != NULL) {
- // Temporarily append a diagonal block to the A matrix, but undo
- // it before returning the matrix to the user.
- A->AppendDiagonal(per_solve_options.D);
- }
-
- // rhs = [b;0] to account for the additional rows in the lhs.
- const int augmented_num_rows =
- num_rows + ((per_solve_options.D != NULL) ? num_cols : 0);
- if (rhs_.rows() != augmented_num_rows) {
- rhs_.resize(augmented_num_rows);
- rhs_.setZero();
- }
- rhs_.head(num_rows) = ConstVectorRef(b, num_rows);
- event_logger.AddEvent("Setup");
-
- // Solve the system.
- VectorRef(x, num_cols) = A->matrix().householderQr().solve(rhs_);
event_logger.AddEvent("Solve");
- if (per_solve_options.D != NULL) {
- // Undo the modifications to the matrix A.
- A->RemoveDiagonal();
- }
-
- // We always succeed, since the QR solver returns the best solution
- // it can. It is the job of the caller to determine if the solution
- // is good enough or not.
- LinearSolver::Summary summary;
- summary.num_iterations = 1;
- summary.termination_type = LINEAR_SOLVER_SUCCESS;
- summary.message = "Success.";
-
- event_logger.AddEvent("TearDown");
return summary;
}
diff --git a/internal/ceres/dense_qr_solver.h b/internal/ceres/dense_qr_solver.h
index 980243b..927dc32 100644
--- a/internal/ceres/dense_qr_solver.h
+++ b/internal/ceres/dense_qr_solver.h
@@ -32,6 +32,9 @@
#ifndef CERES_INTERNAL_DENSE_QR_SOLVER_H_
#define CERES_INTERNAL_DENSE_QR_SOLVER_H_
+#include <memory>
+
+#include "ceres/dense_qr.h"
#include "ceres/internal/eigen.h"
#include "ceres/internal/port.h"
#include "ceres/linear_solver.h"
@@ -105,7 +108,7 @@
const LinearSolver::Options options_;
ColMajorMatrix lhs_;
Vector rhs_;
- Vector work_;
+ std::unique_ptr<DenseQR> dense_qr_;
};
} // namespace internal
diff --git a/internal/ceres/dense_qr_test.cc b/internal/ceres/dense_qr_test.cc
new file mode 100644
index 0000000..c71f40e
--- /dev/null
+++ b/internal/ceres/dense_qr_test.cc
@@ -0,0 +1,114 @@
+// Ceres Solver - A fast non-linear least squares minimizer
+// Copyright 2022 Google Inc. All rights reserved.
+// http://ceres-solver.org/
+//
+// Redistribution and use in source and binary forms, with or without
+// modification, are permitted provided that the following conditions are met:
+//
+// * Redistributions of source code must retain the above copyright notice,
+// this list of conditions and the following disclaimer.
+// * Redistributions in binary form must reproduce the above copyright notice,
+// this list of conditions and the following disclaimer in the documentation
+// and/or other materials provided with the distribution.
+// * Neither the name of Google Inc. nor the names of its contributors may be
+// used to endorse or promote products derived from this software without
+// specific prior written permission.
+//
+// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
+// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+// POSSIBILITY OF SUCH DAMAGE.
+//
+// Author: sameeragarwal@google.com (Sameer Agarwal)
+
+#include "ceres/dense_qr.h"
+
+#include <memory>
+#include <numeric>
+#include <string>
+#include <vector>
+
+#include "Eigen/Dense"
+#include "ceres/internal/eigen.h"
+#include "ceres/linear_solver.h"
+#include "glog/logging.h"
+#include "gmock/gmock.h"
+#include "gtest/gtest.h"
+
+namespace ceres {
+namespace internal {
+
+typedef DenseLinearAlgebraLibraryType Param;
+
+std::string ParamInfoToString(testing::TestParamInfo<Param> info) {
+ return DenseLinearAlgebraLibraryTypeToString(info.param);
+}
+
+class DenseQRTest : public ::testing::TestWithParam<Param> {};
+
+TEST_P(DenseQRTest, FactorAndSolve) {
+ // TODO(sameeragarwal): Convert these tests into type parameterized tests so
+ // that we can test the single and double precision solvers.
+
+ using Scalar = double;
+ using MatrixType = Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic>;
+ using VectorType = Eigen::Matrix<Scalar, Eigen::Dynamic, 1>;
+
+ LinearSolver::Options options;
+ options.dense_linear_algebra_library_type = GetParam();
+ std::unique_ptr<DenseQR> dense_qr = DenseQR::Create(options);
+
+ const int kNumTrials = 10;
+ const int kMinNumCols = 1;
+ const int kMaxNumCols = 10;
+ const int kMinRowsFactor = 1;
+ const int kMaxRowsFactor = 3;
+ for (int num_cols = kMinNumCols; num_cols < kMaxNumCols; ++num_cols) {
+ for (int num_rows = kMinRowsFactor * num_cols;
+ num_rows < kMaxRowsFactor * num_cols;
+ ++num_rows) {
+ for (int trial = 0; trial < kNumTrials; ++trial) {
+ MatrixType lhs = MatrixType::Random(num_rows, num_cols);
+ Vector x = VectorType::Random(num_cols);
+ Vector rhs = lhs * x;
+ Vector actual = Vector::Random(num_cols);
+
+ LinearSolver::Summary summary;
+ summary.termination_type = dense_qr->FactorAndSolve(num_rows,
+ num_cols,
+ lhs.data(),
+ rhs.data(),
+ actual.data(),
+ &summary.message);
+ EXPECT_EQ(summary.termination_type, LINEAR_SOLVER_SUCCESS);
+ EXPECT_NEAR((x - actual).norm() / x.norm(),
+ 0.0,
+ std::numeric_limits<double>::epsilon() * 100)
+ << "\nexpected: " << x.transpose()
+ << "\nactual : " << actual.transpose();
+ }
+ }
+ }
+}
+
+#ifndef CERES_NO_LAPACK
+INSTANTIATE_TEST_SUITE_P(_,
+ DenseQRTest,
+ ::testing::Values(EIGEN, LAPACK),
+ ParamInfoToString);
+#else
+INSTANTIATE_TEST_SUITE_P(_,
+ DenseQRTest,
+ ::testing::Values(EIGEN),
+ ParamInfoToString);
+#endif
+
+} // namespace internal
+} // namespace ceres
diff --git a/internal/ceres/dense_sparse_matrix.cc b/internal/ceres/dense_sparse_matrix.cc
index 53207fe..80ea658 100644
--- a/internal/ceres/dense_sparse_matrix.cc
+++ b/internal/ceres/dense_sparse_matrix.cc
@@ -41,28 +41,10 @@
namespace internal {
DenseSparseMatrix::DenseSparseMatrix(int num_rows, int num_cols)
- : has_diagonal_appended_(false), has_diagonal_reserved_(false) {
- m_.resize(num_rows, num_cols);
- m_.setZero();
-}
-
-DenseSparseMatrix::DenseSparseMatrix(int num_rows,
- int num_cols,
- bool reserve_diagonal)
- : has_diagonal_appended_(false), has_diagonal_reserved_(reserve_diagonal) {
- if (reserve_diagonal) {
- // Allocate enough space for the diagonal.
- m_.resize(num_rows + num_cols, num_cols);
- } else {
- m_.resize(num_rows, num_cols);
- }
- m_.setZero();
-}
+ : m_(Matrix(num_rows, num_cols)) {}
DenseSparseMatrix::DenseSparseMatrix(const TripletSparseMatrix& m)
- : m_(Eigen::MatrixXd::Zero(m.num_rows(), m.num_cols())),
- has_diagonal_appended_(false),
- has_diagonal_reserved_(false) {
+ : m_(Matrix::Zero(m.num_rows(), m.num_cols())) {
const double* values = m.values();
const int* rows = m.rows();
const int* cols = m.cols();
@@ -73,8 +55,7 @@
}
}
-DenseSparseMatrix::DenseSparseMatrix(const ColMajorMatrix& m)
- : m_(m), has_diagonal_appended_(false), has_diagonal_reserved_(false) {}
+DenseSparseMatrix::DenseSparseMatrix(const Matrix& m) : m_(m) {}
void DenseSparseMatrix::SetZero() { m_.setZero(); }
@@ -96,72 +77,22 @@
}
void DenseSparseMatrix::ToDenseMatrix(Matrix* dense_matrix) const {
- *dense_matrix = m_.block(0, 0, num_rows(), num_cols());
+ *dense_matrix = m_;
}
-void DenseSparseMatrix::AppendDiagonal(double* d) {
- CHECK(!has_diagonal_appended_);
- if (!has_diagonal_reserved_) {
- ColMajorMatrix tmp = m_;
- m_.resize(m_.rows() + m_.cols(), m_.cols());
- m_.setZero();
- m_.block(0, 0, tmp.rows(), tmp.cols()) = tmp;
- has_diagonal_reserved_ = true;
- }
-
- m_.bottomLeftCorner(m_.cols(), m_.cols()) =
- ConstVectorRef(d, m_.cols()).asDiagonal();
- has_diagonal_appended_ = true;
-}
-
-void DenseSparseMatrix::RemoveDiagonal() {
- CHECK(has_diagonal_appended_);
- has_diagonal_appended_ = false;
- // Leave the diagonal reserved.
-}
-
-int DenseSparseMatrix::num_rows() const {
- if (has_diagonal_reserved_ && !has_diagonal_appended_) {
- return m_.rows() - m_.cols();
- }
- return m_.rows();
-}
+int DenseSparseMatrix::num_rows() const { return m_.rows(); }
int DenseSparseMatrix::num_cols() const { return m_.cols(); }
-int DenseSparseMatrix::num_nonzeros() const {
- if (has_diagonal_reserved_ && !has_diagonal_appended_) {
- return (m_.rows() - m_.cols()) * m_.cols();
- }
- return m_.rows() * m_.cols();
-}
+int DenseSparseMatrix::num_nonzeros() const { return m_.rows() * m_.cols(); }
-ConstColMajorMatrixRef DenseSparseMatrix::matrix() const {
- return ConstColMajorMatrixRef(
- m_.data(),
- ((has_diagonal_reserved_ && !has_diagonal_appended_)
- ? m_.rows() - m_.cols()
- : m_.rows()),
- m_.cols(),
- Eigen::Stride<Eigen::Dynamic, 1>(m_.rows(), 1));
-}
+const Matrix& DenseSparseMatrix::matrix() const { return m_; }
-ColMajorMatrixRef DenseSparseMatrix::mutable_matrix() {
- return ColMajorMatrixRef(m_.data(),
- ((has_diagonal_reserved_ && !has_diagonal_appended_)
- ? m_.rows() - m_.cols()
- : m_.rows()),
- m_.cols(),
- Eigen::Stride<Eigen::Dynamic, 1>(m_.rows(), 1));
-}
+Matrix* DenseSparseMatrix::mutable_matrix() { return &m_; }
void DenseSparseMatrix::ToTextFile(FILE* file) const {
CHECK(file != nullptr);
- const int active_rows = (has_diagonal_reserved_ && !has_diagonal_appended_)
- ? (m_.rows() - m_.cols())
- : m_.rows();
-
- for (int r = 0; r < active_rows; ++r) {
+ for (int r = 0; r < m_.rows(); ++r) {
for (int c = 0; c < m_.cols(); ++c) {
fprintf(file, "% 10d % 10d %17f\n", r, c, m_(r, c));
}
diff --git a/internal/ceres/dense_sparse_matrix.h b/internal/ceres/dense_sparse_matrix.h
index 94064b3..e2f4a6e 100644
--- a/internal/ceres/dense_sparse_matrix.h
+++ b/internal/ceres/dense_sparse_matrix.h
@@ -48,12 +48,10 @@
// Build a matrix with the same content as the TripletSparseMatrix
// m. This assumes that m does not have any repeated entries.
explicit DenseSparseMatrix(const TripletSparseMatrix& m);
- explicit DenseSparseMatrix(const ColMajorMatrix& m);
-
+ explicit DenseSparseMatrix(const Matrix& m);
DenseSparseMatrix(int num_rows, int num_cols);
- DenseSparseMatrix(int num_rows, int num_cols, bool reserve_diagonal);
- virtual ~DenseSparseMatrix() {}
+ virtual ~DenseSparseMatrix() = default;
// SparseMatrix interface.
void SetZero() final;
@@ -69,37 +67,11 @@
const double* values() const final { return m_.data(); }
double* mutable_values() final { return m_.data(); }
- ConstColMajorMatrixRef matrix() const;
- ColMajorMatrixRef mutable_matrix();
-
- // Only one diagonal can be appended at a time. The diagonal is appended to
- // as a new set of rows, e.g.
- //
- // Original matrix:
- //
- // x x x
- // x x x
- // x x x
- //
- // After append diagonal (1, 2, 3):
- //
- // x x x
- // x x x
- // x x x
- // 1 0 0
- // 0 2 0
- // 0 0 3
- //
- // Calling RemoveDiagonal removes the block. It is a fatal error to append a
- // diagonal to a matrix that already has an appended diagonal, and it is also
- // a fatal error to remove a diagonal from a matrix that has none.
- void AppendDiagonal(double* d);
- void RemoveDiagonal();
+ const Matrix& matrix() const;
+ Matrix* mutable_matrix();
private:
- ColMajorMatrix m_;
- bool has_diagonal_appended_;
- bool has_diagonal_reserved_;
+ Matrix m_;
};
} // namespace internal
diff --git a/internal/ceres/lapack.cc b/internal/ceres/lapack.cc
deleted file mode 100644
index 7ce00a5..0000000
--- a/internal/ceres/lapack.cc
+++ /dev/null
@@ -1,132 +0,0 @@
-// Ceres Solver - A fast non-linear least squares minimizer
-// Copyright 2015 Google Inc. All rights reserved.
-// http://ceres-solver.org/
-//
-// Redistribution and use in source and binary forms, with or without
-// modification, are permitted provided that the following conditions are met:
-//
-// * Redistributions of source code must retain the above copyright notice,
-// this list of conditions and the following disclaimer.
-// * Redistributions in binary form must reproduce the above copyright notice,
-// this list of conditions and the following disclaimer in the documentation
-// and/or other materials provided with the distribution.
-// * Neither the name of Google Inc. nor the names of its contributors may be
-// used to endorse or promote products derived from this software without
-// specific prior written permission.
-//
-// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
-// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
-// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
-// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
-// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
-// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
-// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
-// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
-// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
-// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
-// POSSIBILITY OF SUCH DAMAGE.
-//
-// Author: sameeragarwal@google.com (Sameer Agarwal)
-
-#include "ceres/lapack.h"
-
-#include "ceres/internal/port.h"
-#include "ceres/linear_solver.h"
-#include "glog/logging.h"
-
-#ifndef CERES_NO_LAPACK
-extern "C" void dgels_(char* uplo,
- int* m,
- int* n,
- int* nrhs,
- double* a,
- int* lda,
- double* b,
- int* ldb,
- double* work,
- int* lwork,
- int* info);
-#endif
-
-namespace ceres {
-namespace internal {
-
-int LAPACK::EstimateWorkSizeForQR(int num_rows, int num_cols) {
-#ifdef CERES_NO_LAPACK
- LOG(FATAL) << "Ceres was built without a LAPACK library.";
- return -1;
-#else
- char trans = 'N';
- int nrhs = 1;
- int lwork = -1;
- double work;
- int info = 0;
- dgels_(&trans,
- &num_rows,
- &num_cols,
- &nrhs,
- NULL,
- &num_rows,
- NULL,
- &num_rows,
- &work,
- &lwork,
- &info);
-
- if (info < 0) {
- LOG(FATAL) << "Congratulations, you found a bug in Ceres."
- << "Please report it."
- << "LAPACK::dgels fatal error."
- << "Argument: " << -info << " is invalid.";
- }
- return static_cast<int>(work);
-#endif
-}
-
-LinearSolverTerminationType LAPACK::SolveInPlaceUsingQR(
- int num_rows,
- int num_cols,
- const double* in_lhs,
- int work_size,
- double* work,
- double* rhs_and_solution,
- std::string* message) {
-#ifdef CERES_NO_LAPACK
- LOG(FATAL) << "Ceres was built without a LAPACK library.";
- return LINEAR_SOLVER_FATAL_ERROR;
-#else
- char trans = 'N';
- int m = num_rows;
- int n = num_cols;
- int nrhs = 1;
- int lda = num_rows;
- int ldb = num_rows;
- int info = 0;
- double* lhs = const_cast<double*>(in_lhs);
-
- dgels_(&trans,
- &m,
- &n,
- &nrhs,
- lhs,
- &lda,
- rhs_and_solution,
- &ldb,
- work,
- &work_size,
- &info);
-
- if (info < 0) {
- LOG(FATAL) << "Congratulations, you found a bug in Ceres."
- << "Please report it."
- << "LAPACK::dgels fatal error."
- << "Argument: " << -info << " is invalid.";
- }
-
- *message = "Success.";
- return LINEAR_SOLVER_SUCCESS;
-#endif
-}
-
-} // namespace internal
-} // namespace ceres
diff --git a/internal/ceres/lapack.h b/internal/ceres/lapack.h
deleted file mode 100644
index 5c5bf8b..0000000
--- a/internal/ceres/lapack.h
+++ /dev/null
@@ -1,101 +0,0 @@
-// Ceres Solver - A fast non-linear least squares minimizer
-// Copyright 2015 Google Inc. All rights reserved.
-// http://ceres-solver.org/
-//
-// Redistribution and use in source and binary forms, with or without
-// modification, are permitted provided that the following conditions are met:
-//
-// * Redistributions of source code must retain the above copyright notice,
-// this list of conditions and the following disclaimer.
-// * Redistributions in binary form must reproduce the above copyright notice,
-// this list of conditions and the following disclaimer in the documentation
-// and/or other materials provided with the distribution.
-// * Neither the name of Google Inc. nor the names of its contributors may be
-// used to endorse or promote products derived from this software without
-// specific prior written permission.
-//
-// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
-// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
-// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
-// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
-// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
-// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
-// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
-// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
-// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
-// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
-// POSSIBILITY OF SUCH DAMAGE.
-//
-// Author: sameeragarwal@google.com (Sameer Agarwal)
-
-#ifndef CERES_INTERNAL_LAPACK_H_
-#define CERES_INTERNAL_LAPACK_H_
-
-#include <string>
-
-#include "ceres/internal/port.h"
-#include "ceres/linear_solver.h"
-
-namespace ceres {
-namespace internal {
-
-class LAPACK {
- public:
- // Solve
- //
- // lhs * solution = rhs
- //
- // using a Cholesky factorization. Here
- // lhs is a symmetric positive definite matrix. It is assumed to be
- // column major and only the lower triangular part of the matrix is
- // referenced.
- //
- // This function uses the LAPACK dpotrf and dpotrs routines.
- //
- // The return value and the message string together describe whether
- // the solver terminated successfully or not and if so, what was the
- // reason for failure.
- static LinearSolverTerminationType SolveInPlaceUsingCholesky(
- int num_rows,
- const double* lhs,
- double* rhs_and_solution,
- std::string* message);
-
- // The SolveUsingQR function requires a buffer for its temporary
- // computation. This function given the size of the lhs matrix will
- // return the size of the buffer needed.
- static int EstimateWorkSizeForQR(int num_rows, int num_cols);
-
- // Solve
- //
- // lhs * solution = rhs
- //
- // using a dense QR factorization. lhs is an arbitrary (possibly
- // rectangular) matrix with full column rank.
- //
- // work is an array of size work_size that this routine uses for its
- // temporary storage. The optimal size of this array can be obtained
- // by calling EstimateWorkSizeForQR.
- //
- // When calling, rhs_and_solution contains the rhs, and upon return
- // the first num_col entries are the solution.
- //
- // This function uses the LAPACK dgels routine.
- //
- // The return value and the message string together describe whether
- // the solver terminated successfully or not and if so, what was the
- // reason for failure.
- static LinearSolverTerminationType SolveInPlaceUsingQR(
- int num_rows,
- int num_cols,
- const double* lhs,
- int work_size,
- double* work,
- double* rhs_and_solution,
- std::string* message);
-};
-
-} // namespace internal
-} // namespace ceres
-
-#endif // CERES_INTERNAL_LAPACK_H_
diff --git a/internal/ceres/trust_region_minimizer_test.cc b/internal/ceres/trust_region_minimizer_test.cc
index cad30eb..8f544d6 100644
--- a/internal/ceres/trust_region_minimizer_test.cc
+++ b/internal/ceres/trust_region_minimizer_test.cc
@@ -131,7 +131,7 @@
dense_jacobian = down_cast<DenseSparseMatrix*>(jacobian);
dense_jacobian->SetZero();
- ColMajorMatrixRef jacobian_matrix = dense_jacobian->mutable_matrix();
+ Matrix& jacobian_matrix = *(dense_jacobian->mutable_matrix());
CHECK_EQ(jacobian_matrix.cols(), num_active_cols_);
int column_index = 0;
