Better error checking and reporting for linear solvers.
A lot of error checking cruft has accumulated over the years
in the various linear solvers. This change makes the error reporting
more robust and consistent across the various solvers.
Preconditioners are not covered by this change and will be the
subject of a future change.
Change-Id: Ibeb2572a1e67758953dde8d12e3abc6d1df9052d
diff --git a/internal/ceres/cgnr_solver.cc b/internal/ceres/cgnr_solver.cc
index 9b8f980..88e61d9 100644
--- a/internal/ceres/cgnr_solver.cc
+++ b/internal/ceres/cgnr_solver.cc
@@ -33,6 +33,7 @@
#include "ceres/block_jacobi_preconditioner.h"
#include "ceres/cgnr_linear_operator.h"
#include "ceres/conjugate_gradients_solver.h"
+#include "ceres/internal/eigen.h"
#include "ceres/linear_solver.h"
#include "ceres/wall_time.h"
#include "glog/logging.h"
@@ -43,6 +44,10 @@
CgnrSolver::CgnrSolver(const LinearSolver::Options& options)
: options_(options),
preconditioner_(NULL) {
+ if (options_.preconditioner_type != JACOBI &&
+ options_.preconditioner_type != IDENTITY) {
+ LOG(FATAL) << "CGNR only supports IDENTITY and JACOBI preconditioners.";
+ }
}
LinearSolver::Summary CgnrSolver::SolveImpl(
@@ -53,9 +58,9 @@
EventLogger event_logger("CgnrSolver::Solve");
// Form z = Atb.
- scoped_array<double> z(new double[A->num_cols()]);
- std::fill(z.get(), z.get() + A->num_cols(), 0.0);
- A->LeftMultiply(b, z.get());
+ Vector z(A->num_cols());
+ z.setZero();
+ A->LeftMultiply(b, z.data());
// Precondition if necessary.
LinearSolver::PerSolveOptions cg_per_solve_options = per_solve_options;
@@ -65,20 +70,17 @@
}
preconditioner_->Update(*A, per_solve_options.D);
cg_per_solve_options.preconditioner = preconditioner_.get();
- } else if (options_.preconditioner_type != IDENTITY) {
- LOG(FATAL) << "CGNR only supports IDENTITY and JACOBI preconditioners.";
}
// Solve (AtA + DtD)x = z (= Atb).
- std::fill(x, x + A->num_cols(), 0.0);
+ VectorRef(x, A->num_cols()).setZero();
CgnrLinearOperator lhs(*A, per_solve_options.D);
event_logger.AddEvent("Setup");
ConjugateGradientsSolver conjugate_gradient_solver(options_);
LinearSolver::Summary summary =
- conjugate_gradient_solver.Solve(&lhs, z.get(), cg_per_solve_options, x);
+ conjugate_gradient_solver.Solve(&lhs, z.data(), cg_per_solve_options, x);
event_logger.AddEvent("Solve");
-
return summary;
}
diff --git a/internal/ceres/conjugate_gradients_solver.cc b/internal/ceres/conjugate_gradients_solver.cc
index ae8e877..ffea501 100644
--- a/internal/ceres/conjugate_gradients_solver.cc
+++ b/internal/ceres/conjugate_gradients_solver.cc
@@ -44,6 +44,7 @@
#include "ceres/fpclassify.h"
#include "ceres/internal/eigen.h"
#include "ceres/linear_operator.h"
+#include "ceres/stringprintf.h"
#include "ceres/types.h"
#include "glog/logging.h"
@@ -55,9 +56,6 @@
return ((x == 0.0) || (IsInfinite(x)));
}
-// Constant used in the MATLAB implementation ~ 2 * eps.
-const double kEpsilon = 2.2204e-16;
-
} // namespace
ConjugateGradientsSolver::ConjugateGradientsSolver(
@@ -77,16 +75,18 @@
LinearSolver::Summary summary;
summary.termination_type = MAX_ITERATIONS;
+ summary.status = "Maximum number of iterations reached.";
summary.num_iterations = 0;
- int num_cols = A->num_cols();
+ const int num_cols = A->num_cols();
VectorRef xref(x, num_cols);
ConstVectorRef bref(b, num_cols);
- double norm_b = bref.norm();
+ const double norm_b = bref.norm();
if (norm_b == 0.0) {
xref.setZero();
summary.termination_type = TOLERANCE;
+ summary.status = "Convergence. |b| = 0.";
return summary;
}
@@ -95,15 +95,16 @@
Vector z(num_cols);
Vector tmp(num_cols);
- double tol_r = per_solve_options.r_tolerance * norm_b;
+ const double tol_r = per_solve_options.r_tolerance * norm_b;
tmp.setZero();
A->RightMultiply(x, tmp.data());
r = bref - tmp;
double norm_r = r.norm();
-
if (norm_r <= tol_r) {
summary.termination_type = TOLERANCE;
+ summary.status =
+ StringPrintf("Convergence. |r| = %e <= %e.", norm_r, tol_r);
return summary;
}
@@ -115,8 +116,6 @@
for (summary.num_iterations = 1;
summary.num_iterations < options_.max_num_iterations;
++summary.num_iterations) {
- VLOG(3) << "cg iteration " << summary.num_iterations;
-
// Apply preconditioner
if (per_solve_options.preconditioner != NULL) {
z.setZero();
@@ -127,10 +126,9 @@
double last_rho = rho;
rho = r.dot(z);
-
if (IsZeroOrInfinity(rho)) {
- LOG(ERROR) << "Numerical failure. rho = " << rho;
summary.termination_type = FAILURE;
+ summary.status = StringPrintf("Numerical failure. rho = r'z = %e.", rho);
break;
};
@@ -139,8 +137,9 @@
} else {
double beta = rho / last_rho;
if (IsZeroOrInfinity(beta)) {
- LOG(ERROR) << "Numerical failure. beta = " << beta;
summary.termination_type = FAILURE;
+ summary.status = StringPrintf(
+ "Numerical failure. beta = rho_n / rho_{n-1} = %e.", beta);
break;
}
p = z + beta * p;
@@ -149,18 +148,18 @@
Vector& q = z;
q.setZero();
A->RightMultiply(p.data(), q.data());
- double pq = p.dot(q);
-
+ const double pq = p.dot(q);
if ((pq <= 0) || IsInfinite(pq)) {
- LOG(ERROR) << "Numerical failure. pq = " << pq;
summary.termination_type = FAILURE;
+ summary.status = StringPrintf("Numerical failure. p'q = %e.", pq);
break;
}
- double alpha = rho / pq;
+ const double alpha = rho / pq;
if (IsInfinite(alpha)) {
- LOG(ERROR) << "Numerical failure. alpha " << alpha;
summary.termination_type = FAILURE;
+ summary.status =
+ StringPrintf("Numerical failure. alpha = rho / pq = %e", alpha);
break;
}
@@ -183,7 +182,7 @@
// Quadratic model based termination.
// Q1 = x'Ax - 2 * b' x.
- double Q1 = -1.0 * xref.dot(bref + r);
+ const double Q1 = -1.0 * xref.dot(bref + r);
// For PSD matrices A, let
//
@@ -207,21 +206,23 @@
// Journal of Computational and Applied Mathematics,
// 124(1-2), 45-59, 2000.
//
- double zeta = summary.num_iterations * (Q1 - Q0) / Q1;
- VLOG(3) << "Q termination: zeta " << zeta
- << " " << per_solve_options.q_tolerance;
+ const double zeta = summary.num_iterations * (Q1 - Q0) / Q1;
if (zeta < per_solve_options.q_tolerance) {
summary.termination_type = TOLERANCE;
+ summary.status =
+ StringPrintf("Convergence: zeta = %e < %e",
+ zeta,
+ per_solve_options.q_tolerance);
break;
}
Q0 = Q1;
// Residual based termination.
norm_r = r. norm();
- VLOG(3) << "R termination: norm_r " << norm_r
- << " " << tol_r;
if (norm_r <= tol_r) {
summary.termination_type = TOLERANCE;
+ summary.status =
+ StringPrintf("Convergence. |r| = %e <= %e.", norm_r, tol_r);
break;
}
}
diff --git a/internal/ceres/covariance_impl.cc b/internal/ceres/covariance_impl.cc
index 6fe7a80..6c79fa5 100644
--- a/internal/ceres/covariance_impl.cc
+++ b/internal/ceres/covariance_impl.cc
@@ -165,9 +165,9 @@
}
if (offset == row_size) {
- LOG(WARNING) << "Unable to find covariance block for "
- << original_parameter_block1 << " "
- << original_parameter_block2;
+ LOG(ERROR) << "Unable to find covariance block for "
+ << original_parameter_block1 << " "
+ << original_parameter_block2;
return false;
}
@@ -441,18 +441,24 @@
cholmod_jacobian_view.sorted = 1;
cholmod_jacobian_view.packed = 1;
- cholmod_factor* factor = ss.AnalyzeCholesky(&cholmod_jacobian_view);
+ string status;
+ cholmod_factor* factor = ss.AnalyzeCholesky(&cholmod_jacobian_view, &status);
event_logger.AddEvent("Symbolic Factorization");
if (factor == NULL) {
+ LOG(ERROR) << "Covariance estimation failed. "
+ << "CHOLMOD symbolic cholesky factorization returned with: "
+ << status;
return false;
}
LinearSolverTerminationType termination_type =
- ss.Cholesky(&cholmod_jacobian_view, factor);
+ ss.Cholesky(&cholmod_jacobian_view, factor, &status);
event_logger.AddEvent("Numeric Factorization");
-
if (termination_type != TOLERANCE) {
- LOG(WARNING) << "Cholesky factorization failed.";
+ LOG(ERROR) << "Covariance estimation failed. "
+ << "CHOLMOD numeric cholesky factorization returned with: "
+ << status;
+ ss.Free(factor);
return false;
}
@@ -461,11 +467,11 @@
if (reciprocal_condition_number <
options_.min_reciprocal_condition_number) {
- LOG(WARNING) << "Cholesky factorization of J'J is not reliable. "
- << "Reciprocal condition number: "
- << reciprocal_condition_number << " "
- << "min_reciprocal_condition_number : "
- << options_.min_reciprocal_condition_number;
+ LOG(ERROR) << "Cholesky factorization of J'J is not reliable. "
+ << "Reciprocal condition number: "
+ << reciprocal_condition_number << " "
+ << "min_reciprocal_condition_number : "
+ << options_.min_reciprocal_condition_number;
ss.Free(factor);
return false;
}
@@ -687,9 +693,9 @@
CHECK_NOTNULL(R);
if (rank < cholmod_jacobian.ncol) {
- LOG(WARNING) << "Jacobian matrix is rank deficient."
- << "Number of columns: " << cholmod_jacobian.ncol
- << " rank: " << rank;
+ LOG(ERROR) << "Jacobian matrix is rank deficient. "
+ << "Number of columns: " << cholmod_jacobian.ncol
+ << " rank: " << rank;
free(permutation);
cholmod_l_free_sparse(&R, &cc);
cholmod_l_finish(&cc);
@@ -813,11 +819,11 @@
if (automatic_truncation) {
break;
} else {
- LOG(WARNING) << "Cholesky factorization of J'J is not reliable. "
- << "Reciprocal condition number: "
- << singular_value_ratio * singular_value_ratio << " "
- << "min_reciprocal_condition_number : "
- << options_.min_reciprocal_condition_number;
+ LOG(ERROR) << "Cholesky factorization of J'J is not reliable. "
+ << "Reciprocal condition number: "
+ << singular_value_ratio * singular_value_ratio << " "
+ << "min_reciprocal_condition_number : "
+ << options_.min_reciprocal_condition_number;
return false;
}
}
diff --git a/internal/ceres/dense_normal_cholesky_solver.cc b/internal/ceres/dense_normal_cholesky_solver.cc
index fbf3cbe..e0fe86e 100644
--- a/internal/ceres/dense_normal_cholesky_solver.cc
+++ b/internal/ceres/dense_normal_cholesky_solver.cc
@@ -96,8 +96,17 @@
LinearSolver::Summary summary;
summary.num_iterations = 1;
summary.termination_type = TOLERANCE;
- VectorRef(x, num_cols) =
- lhs.selfadjointView<Eigen::Upper>().llt().solve(rhs);
+ Eigen::LLT<Matrix, Eigen::Upper> llt = lhs.selfadjointView<Eigen::Upper>().llt();
+
+ if (llt.info() != Eigen::Success) {
+ summary.termination_type = FAILURE;
+ summary.status = "Eigen LLT decomposition failed.";
+ } else {
+ summary.termination_type = TOLERANCE;
+ summary.status = "Success.";
+ }
+
+ VectorRef(x, num_cols) = llt.solve(rhs);
event_logger.AddEvent("Solve");
return summary;
}
@@ -142,14 +151,13 @@
A->matrix().transpose() * ConstVectorRef(b, A->num_rows());
event_logger.AddEvent("Product");
- const int info = LAPACK::SolveInPlaceUsingCholesky(num_cols, lhs.data(), x);
- event_logger.AddEvent("Solve");
-
LinearSolver::Summary summary;
summary.num_iterations = 1;
- summary.termination_type = info == 0 ? TOLERANCE : FAILURE;
-
- event_logger.AddEvent("TearDown");
+ summary.termination_type = LAPACK::SolveInPlaceUsingCholesky(num_cols,
+ lhs.data(),
+ x,
+ &summary.status);
+ event_logger.AddEvent("Solve");
return summary;
}
} // namespace internal
diff --git a/internal/ceres/dense_qr_solver.cc b/internal/ceres/dense_qr_solver.cc
index d76d58b..fcc87d2 100644
--- a/internal/ceres/dense_qr_solver.cc
+++ b/internal/ceres/dense_qr_solver.cc
@@ -60,6 +60,7 @@
return SolveUsingLAPACK(A, b, per_solve_options, x);
}
}
+
LinearSolver::Summary DenseQRSolver::SolveUsingLAPACK(
DenseSparseMatrix* A,
const double* b,
@@ -100,21 +101,18 @@
work_.resize(work_size);
}
- const int info = LAPACK::SolveUsingQR(lhs_.rows(),
- lhs_.cols(),
- lhs_.data(),
- work_.rows(),
- work_.data(),
- rhs_.data());
- event_logger.AddEvent("Solve");
-
LinearSolver::Summary summary;
summary.num_iterations = 1;
- if (info == 0) {
+ summary.termination_type = LAPACK::SolveInPlaceUsingQR(lhs_.rows(),
+ lhs_.cols(),
+ lhs_.data(),
+ work_.rows(),
+ work_.data(),
+ rhs_.data(),
+ &summary.status);
+ event_logger.AddEvent("Solve");
+ if (summary.termination_type == TOLERANCE) {
VectorRef(x, num_cols) = rhs_.head(num_cols);
- summary.termination_type = TOLERANCE;
- } else {
- summary.termination_type = FAILURE;
}
event_logger.AddEvent("TearDown");
@@ -162,6 +160,7 @@
LinearSolver::Summary summary;
summary.num_iterations = 1;
summary.termination_type = TOLERANCE;
+ summary.status = "Success.";
event_logger.AddEvent("TearDown");
return summary;
diff --git a/internal/ceres/iterative_schur_complement_solver.cc b/internal/ceres/iterative_schur_complement_solver.cc
index e9c54c8..4777d02 100644
--- a/internal/ceres/iterative_schur_complement_solver.cc
+++ b/internal/ceres/iterative_schur_complement_solver.cc
@@ -153,26 +153,26 @@
preconditioner_->Update(*A, per_solve_options.D);
cg_per_solve_options.preconditioner = preconditioner_.get();
}
-
event_logger.AddEvent("Setup");
LinearSolver::Summary cg_summary;
cg_summary.num_iterations = 0;
cg_summary.termination_type = FAILURE;
+ // TODO(sameeragarwal): Refactor preconditioners to return a more
+ // sane status.
+ cg_summary.status = "Preconditioner update failed.";
if (preconditioner_update_was_successful) {
cg_summary = cg_solver.Solve(schur_complement_.get(),
schur_complement_->rhs().data(),
cg_per_solve_options,
reduced_linear_system_solution_.data());
- if (cg_summary.termination_type != FAILURE) {
+ if (cg_summary.termination_type != FAILURE &&
+ cg_summary.termination_type != FATAL_ERROR) {
schur_complement_->BackSubstitute(
reduced_linear_system_solution_.data(), x);
}
}
-
- VLOG(2) << "CG Iterations : " << cg_summary.num_iterations;
-
event_logger.AddEvent("Solve");
return cg_summary;
}
diff --git a/internal/ceres/lapack.cc b/internal/ceres/lapack.cc
index e93e05f..0061433 100644
--- a/internal/ceres/lapack.cc
+++ b/internal/ceres/lapack.cc
@@ -29,6 +29,9 @@
// Author: sameeragarwal@google.com (Sameer Agarwal)
#include "ceres/lapack.h"
+
+#include "ceres/internal/port.h"
+#include "ceres/linear_solver.h"
#include "glog/logging.h"
// C interface to the LAPACK Cholesky factorization and triangular solve.
@@ -63,12 +66,14 @@
namespace ceres {
namespace internal {
-int LAPACK::SolveInPlaceUsingCholesky(int num_rows,
- const double* in_lhs,
- double* rhs_and_solution) {
+LinearSolverTerminationType LAPACK::SolveInPlaceUsingCholesky(
+ int num_rows,
+ const double* in_lhs,
+ double* rhs_and_solution,
+ string* status) {
#ifdef CERES_NO_LAPACK
LOG(FATAL) << "Ceres was built without a BLAS library.";
- return -1;
+ return FATAL_ERROR;
#else
char uplo = 'L';
int n = num_rows;
@@ -77,17 +82,33 @@
double* lhs = const_cast<double*>(in_lhs);
dpotrf_(&uplo, &n, lhs, &n, &info);
- if (info != 0) {
- LOG(INFO) << "Cholesky factorization (dpotrf) failed: " << info;
- return info;
+ if (info < 0) {
+ LOG(FATAL) << "Congratulations, you found a bug in Ceres."
+ << "Please report it."
+ << "LAPACK::dpotrf fatal error."
+ << "Argument: " << -info << " is invalid.";
+ return FATAL_ERROR;
+ }
+
+ if (info > 0) {
+ *status =
+ StringPrintf(
+ "LAPACK::dpotrf numerical failure. "
+ "The leading minor of order %d is not positive definite.", info);
+ return FAILURE;
}
dpotrs_(&uplo, &n, &nrhs, lhs, &n, rhs_and_solution, &n, &info);
- if (info != 0) {
- LOG(INFO) << "Triangular solve (dpotrs) failed: " << info;
+ if (info < 0) {
+ LOG(FATAL) << "Congratulations, you found a bug in Ceres."
+ << "Please report it."
+ << "LAPACK::dpotrs fatal error."
+ << "Argument: " << -info << " is invalid.";
+ return FATAL_ERROR;
}
- return info;
+ *status = "Success";
+ return TOLERANCE;
#endif
};
@@ -113,20 +134,27 @@
&lwork,
&info);
- CHECK_EQ(info, 0);
+ 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
}
-int LAPACK::SolveUsingQR(int num_rows,
- int num_cols,
- const double* in_lhs,
- int work_size,
- double* work,
- double* rhs_and_solution) {
+LinearSolverTerminationType LAPACK::SolveInPlaceUsingQR(
+ int num_rows,
+ int num_cols,
+ const double* in_lhs,
+ int work_size,
+ double* work,
+ double* rhs_and_solution,
+ string* status) {
#ifdef CERES_NO_LAPACK
LOG(FATAL) << "Ceres was built without a LAPACK library.";
- return -1;
+ return FATAL_ERROR;
#else
char trans = 'N';
int m = num_rows;
@@ -149,7 +177,15 @@
&work_size,
&info);
- return info;
+ if (info < 0) {
+ LOG(FATAL) << "Congratulations, you found a bug in Ceres."
+ << "Please report it."
+ << "LAPACK::dgels fatal error."
+ << "Argument: " << -info << " is invalid.";
+ }
+
+ *status = "Success.";
+ return TOLERANCE;
#endif
}
diff --git a/internal/ceres/lapack.h b/internal/ceres/lapack.h
index 4f3a88c..53a33e1 100644
--- a/internal/ceres/lapack.h
+++ b/internal/ceres/lapack.h
@@ -31,6 +31,10 @@
#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 {
@@ -47,10 +51,14 @@
//
// This function uses the LAPACK dpotrf and dpotrs routines.
//
- // The return value is zero if the solve is successful.
- static int SolveInPlaceUsingCholesky(int num_rows,
- const double* lhs,
- double* rhs_and_solution);
+ // The return value and the status 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,
+ string* status);
// The SolveUsingQR function requires a buffer for its temporary
// computation. This function given the size of the lhs matrix will
@@ -73,13 +81,17 @@
//
// This function uses the LAPACK dgels routine.
//
- // The return value is zero if the solve is successful.
- static int SolveUsingQR(int num_rows,
- int num_cols,
- const double* lhs,
- int work_size,
- double* work,
- double* rhs_and_solution);
+ // The return value and the status 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,
+ string* status);
};
} // namespace internal
diff --git a/internal/ceres/linear_solver.h b/internal/ceres/linear_solver.h
index 4af586a..fb50d7e 100644
--- a/internal/ceres/linear_solver.h
+++ b/internal/ceres/linear_solver.h
@@ -112,10 +112,8 @@
}
LinearSolverType type;
-
PreconditionerType preconditioner_type;
VisibilityClusteringType visibility_clustering_type;
-
DenseLinearAlgebraLibraryType dense_linear_algebra_library_type;
SparseLinearAlgebraLibraryType sparse_linear_algebra_library_type;
diff --git a/internal/ceres/schur_complement_solver.cc b/internal/ceres/schur_complement_solver.cc
index 91dd455..ff51ab2 100644
--- a/internal/ceres/schur_complement_solver.cc
+++ b/internal/ceres/schur_complement_solver.cc
@@ -75,24 +75,19 @@
fill(x, x + A->num_cols(), 0.0);
event_logger.AddEvent("Setup");
- LinearSolver::Summary summary;
- summary.num_iterations = 1;
- summary.termination_type = FAILURE;
eliminator_->Eliminate(A, b, per_solve_options.D, lhs_.get(), rhs_.get());
event_logger.AddEvent("Eliminate");
double* reduced_solution = x + A->num_cols() - lhs_->num_cols();
- summary.termination_type = SolveReducedLinearSystem(reduced_solution);
+ const LinearSolver::Summary summary =
+ SolveReducedLinearSystem(reduced_solution);
event_logger.AddEvent("ReducedSolve");
- if (summary.termination_type != TOLERANCE) {
- return summary;
+ if (summary.termination_type == TOLERANCE) {
+ eliminator_->BackSubstitute(A, b, per_solve_options.D, reduced_solution, x);
+ event_logger.AddEvent("BackSubstitute");
}
- eliminator_->BackSubstitute(A, b, per_solve_options.D, reduced_solution, x);
- summary.termination_type = TOLERANCE;
-
- event_logger.AddEvent("BackSubstitute");
return summary;
}
@@ -117,8 +112,13 @@
// Solve the system Sx = r, assuming that the matrix S is stored in a
// BlockRandomAccessDenseMatrix. The linear system is solved using
// Eigen's Cholesky factorization.
-LinearSolverTerminationType
+LinearSolver::Summary
DenseSchurComplementSolver::SolveReducedLinearSystem(double* solution) {
+ LinearSolver::Summary summary;
+ summary.num_iterations = 0;
+ summary.termination_type = TOLERANCE;
+ summary.status = "Success.";
+
const BlockRandomAccessDenseMatrix* m =
down_cast<const BlockRandomAccessDenseMatrix*>(lhs());
const int num_rows = m->num_rows();
@@ -126,29 +126,33 @@
// The case where there are no f blocks, and the system is block
// diagonal.
if (num_rows == 0) {
- return TOLERANCE;
+ return summary;
}
+ summary.num_iterations = 1;
+
if (options().dense_linear_algebra_library_type == EIGEN) {
- // TODO(sameeragarwal): Add proper error handling; this completely ignores
- // the quality of the solution to the solve.
- VectorRef(solution, num_rows) =
+ Eigen::LLT<Matrix, Eigen::Upper> llt =
ConstMatrixRef(m->values(), num_rows, num_rows)
.selfadjointView<Eigen::Upper>()
- .llt()
- .solve(ConstVectorRef(rhs(), num_rows));
- return TOLERANCE;
+ .llt();
+ if (llt.info() != Eigen::Success) {
+ summary.termination_type = FAILURE;
+ summary.status = "Eigen LLT decomposition failed.";
+ return summary;
+ }
+
+ VectorRef(solution, num_rows) = llt.solve(ConstVectorRef(rhs(), num_rows));
+ } else {
+ VectorRef(solution, num_rows) = ConstVectorRef(rhs(), num_rows);
+ summary.termination_type =
+ LAPACK::SolveInPlaceUsingCholesky(num_rows,
+ m->values(),
+ solution,
+ &summary.status);
}
- VectorRef(solution, num_rows) = ConstVectorRef(rhs(), num_rows);
- const int info = LAPACK::SolveInPlaceUsingCholesky(num_rows,
- m->values(),
- solution);
- if (info == 0) {
- return TOLERANCE;
- } else {
- return FAILURE;
- }
+ return summary;
}
#if !defined(CERES_NO_SUITESPARSE) || !defined(CERES_NO_CXSPARE)
@@ -247,7 +251,7 @@
set_rhs(new double[lhs()->num_rows()]);
}
-LinearSolverTerminationType
+LinearSolver::Summary
SparseSchurComplementSolver::SolveReducedLinearSystem(double* solution) {
switch (options().sparse_linear_algebra_library_type) {
case SUITE_SPARSE:
@@ -259,30 +263,33 @@
<< options().sparse_linear_algebra_library_type;
}
- LOG(FATAL) << "Unknown sparse linear algebra library : "
- << options().sparse_linear_algebra_library_type;
- return FATAL_ERROR;
+ return LinearSolver::Summary();
}
#ifndef CERES_NO_SUITESPARSE
// Solve the system Sx = r, assuming that the matrix S is stored in a
// BlockRandomAccessSparseMatrix. The linear system is solved using
// CHOLMOD's sparse cholesky factorization routines.
-LinearSolverTerminationType
+LinearSolver::Summary
SparseSchurComplementSolver::SolveReducedLinearSystemUsingSuiteSparse(
double* solution) {
+ LinearSolver::Summary summary;
+ summary.num_iterations = 0;
+ summary.termination_type = TOLERANCE;
+ summary.status = "Success.";
+
TripletSparseMatrix* tsm =
const_cast<TripletSparseMatrix*>(
down_cast<const BlockRandomAccessSparseMatrix*>(lhs())->matrix());
-
const int num_rows = tsm->num_rows();
// The case where there are no f blocks, and the system is block
// diagonal.
if (num_rows == 0) {
- return TOLERANCE;
+ return summary;
}
+ summary.num_iterations = 1;
cholmod_sparse* cholmod_lhs = NULL;
if (options().use_postordering) {
// If we are going to do a full symbolic analysis of the schur
@@ -295,7 +302,10 @@
cholmod_lhs->stype = 1;
if (factor_ == NULL) {
- factor_ = ss_.BlockAnalyzeCholesky(cholmod_lhs, blocks_, blocks_);
+ factor_ = ss_.BlockAnalyzeCholesky(cholmod_lhs,
+ blocks_,
+ blocks_,
+ &summary.status);
}
} else {
// If we are going to use the natural ordering (i.e. rely on the
@@ -308,43 +318,47 @@
cholmod_lhs->stype = -1;
if (factor_ == NULL) {
- factor_ = ss_.AnalyzeCholeskyWithNaturalOrdering(cholmod_lhs);
+ factor_ = ss_.AnalyzeCholeskyWithNaturalOrdering(cholmod_lhs,
+ &summary.status);
}
}
if (factor_ == NULL) {
ss_.Free(cholmod_lhs);
- return FATAL_ERROR;
+ summary.termination_type = FATAL_ERROR;
+ return summary;
+ }
+
+ summary.termination_type =
+ ss_.Cholesky(cholmod_lhs, factor_, &summary.status);
+ if (summary.termination_type != TOLERANCE) {
+ return summary;
}
cholmod_dense* cholmod_rhs =
ss_.CreateDenseVector(const_cast<double*>(rhs()), num_rows, num_rows);
-
- LinearSolverTerminationType status = ss_.Cholesky(cholmod_lhs, factor_);
- if (status != TOLERANCE) {
- return status;
- }
-
- cholmod_dense* cholmod_solution = ss_.Solve(factor_, cholmod_rhs);
+ cholmod_dense* cholmod_solution = ss_.Solve(factor_,
+ cholmod_rhs,
+ &summary.status);
ss_.Free(cholmod_lhs);
ss_.Free(cholmod_rhs);
if (cholmod_solution == NULL) {
- LOG(WARNING) << "CHOLMOD solve failed.";
- return FAILURE;
+ summary.termination_type = FAILURE;
+ return summary;
}
VectorRef(solution, num_rows)
= VectorRef(static_cast<double*>(cholmod_solution->x), num_rows);
ss_.Free(cholmod_solution);
- return TOLERANCE;
+ return summary;
}
#else
-LinearSolverTerminationType
+LinearSolver::Summary
SparseSchurComplementSolver::SolveReducedLinearSystemUsingSuiteSparse(
double* solution) {
LOG(FATAL) << "No SuiteSparse support in Ceres.";
- return FATAL_ERROR;
+ return LinearSolver::Summary();
}
#endif // CERES_NO_SUITESPARSE
@@ -352,20 +366,24 @@
// Solve the system Sx = r, assuming that the matrix S is stored in a
// BlockRandomAccessSparseMatrix. The linear system is solved using
// CXSparse's sparse cholesky factorization routines.
-LinearSolverTerminationType
+LinearSolver::Summary
SparseSchurComplementSolver::SolveReducedLinearSystemUsingCXSparse(
double* solution) {
+ LinearSolver::Summary summary;
+ summary.num_iterations = 0;
+ summary.termination_type = TOLERANCE;
+ summary.status = "Success.";
+
// Extract the TripletSparseMatrix that is used for actually storing S.
TripletSparseMatrix* tsm =
const_cast<TripletSparseMatrix*>(
down_cast<const BlockRandomAccessSparseMatrix*>(lhs())->matrix());
-
const int num_rows = tsm->num_rows();
// The case where there are no f blocks, and the system is block
// diagonal.
if (num_rows == 0) {
- return TOLERANCE;
+ return summary;
}
cs_di* lhs = CHECK_NOTNULL(cxsparse_.CreateSparseMatrix(tsm));
@@ -373,26 +391,26 @@
// Compute symbolic factorization if not available.
if (cxsparse_factor_ == NULL) {
- cxsparse_factor_ =
- CHECK_NOTNULL(cxsparse_.BlockAnalyzeCholesky(lhs, blocks_, blocks_));
+ cxsparse_factor_ = cxsparse_.BlockAnalyzeCholesky(lhs, blocks_, blocks_);
}
- // Solve the linear system.
- bool ok = cxsparse_.SolveCholesky(lhs, cxsparse_factor_, solution);
+ if (cxsparse_factor_ == NULL) {
+ summary.termination_type = FATAL_ERROR;
+ summary.status = "CXSparse failure. Unable to find symbolic factorization.";
+ } else if (!cxsparse_.SolveCholesky(lhs, cxsparse_factor_, solution)) {
+ summary.termination_type = FAILURE;
+ summary.status = "CXSparse::SolveCholesky failed.";
+ }
cxsparse_.Free(lhs);
- if (ok) {
- return TOLERANCE;
- } else {
- return FAILURE;
- }
+ return summary;
}
#else
-LinearSolverTerminationType
+LinearSolver::Summary
SparseSchurComplementSolver::SolveReducedLinearSystemUsingCXSparse(
double* solution) {
LOG(FATAL) << "No CXSparse support in Ceres.";
- return FATAL_ERROR;
+ return LinearSolver::Summary();
}
#endif // CERES_NO_CXPARSE
diff --git a/internal/ceres/schur_complement_solver.h b/internal/ceres/schur_complement_solver.h
index 3303205..a997851 100644
--- a/internal/ceres/schur_complement_solver.h
+++ b/internal/ceres/schur_complement_solver.h
@@ -126,7 +126,7 @@
private:
virtual void InitStorage(const CompressedRowBlockStructure* bs) = 0;
- virtual LinearSolverTerminationType SolveReducedLinearSystem(
+ virtual LinearSolver::Summary SolveReducedLinearSystem(
double* solution) = 0;
LinearSolver::Options options_;
@@ -147,7 +147,7 @@
private:
virtual void InitStorage(const CompressedRowBlockStructure* bs);
- virtual LinearSolverTerminationType SolveReducedLinearSystem(
+ virtual LinearSolver::Summary SolveReducedLinearSystem(
double* solution);
CERES_DISALLOW_COPY_AND_ASSIGN(DenseSchurComplementSolver);
@@ -162,11 +162,11 @@
private:
virtual void InitStorage(const CompressedRowBlockStructure* bs);
- virtual LinearSolverTerminationType SolveReducedLinearSystem(
+ virtual LinearSolver::Summary SolveReducedLinearSystem(
double* solution);
- LinearSolverTerminationType SolveReducedLinearSystemUsingSuiteSparse(
+ LinearSolver::Summary SolveReducedLinearSystemUsingSuiteSparse(
double* solution);
- LinearSolverTerminationType SolveReducedLinearSystemUsingCXSparse(
+ LinearSolver::Summary SolveReducedLinearSystemUsingCXSparse(
double* solution);
// Size of the blocks in the Schur complement.
diff --git a/internal/ceres/sparse_normal_cholesky_solver.cc b/internal/ceres/sparse_normal_cholesky_solver.cc
index 697adc1..2d03c44 100644
--- a/internal/ceres/sparse_normal_cholesky_solver.cc
+++ b/internal/ceres/sparse_normal_cholesky_solver.cc
@@ -102,6 +102,9 @@
LinearSolver::Summary summary;
summary.num_iterations = 1;
+ summary.termination_type = TOLERANCE;
+ summary.status = "Success.";
+
const int num_cols = A->num_cols();
Vector Atb = Vector::Zero(num_cols);
A->LeftMultiply(b, Atb.data());
@@ -137,21 +140,22 @@
// Compute symbolic factorization if not available.
if (cxsparse_factor_ == NULL) {
if (options_.use_postordering) {
- cxsparse_factor_ =
- CHECK_NOTNULL(cxsparse_.BlockAnalyzeCholesky(AtA,
- A->col_blocks(),
- A->col_blocks()));
+ cxsparse_factor_ = cxsparse_.BlockAnalyzeCholesky(AtA,
+ A->col_blocks(),
+ A->col_blocks());
} else {
- cxsparse_factor_ =
- CHECK_NOTNULL(cxsparse_.AnalyzeCholeskyWithNaturalOrdering(AtA));
+ cxsparse_factor_ = cxsparse_.AnalyzeCholeskyWithNaturalOrdering(AtA);
}
}
event_logger.AddEvent("Analysis");
- // Solve the linear system.
- if (cxsparse_.SolveCholesky(AtA, cxsparse_factor_, Atb.data())) {
+ if (cxsparse_factor_ == NULL) {
+ summary.termination_type = FATAL_ERROR;
+ summary.status = "CXSparse failure. Unable to find symbolic factorization.";
+ } else if (cxsparse_.SolveCholesky(AtA, cxsparse_factor_, Atb.data())) {
VectorRef(x, Atb.rows()) = Atb;
- summary.termination_type = TOLERANCE;
+ } else {
+ summary.termination_type = FAILURE;
}
event_logger.AddEvent("Solve");
@@ -179,32 +183,34 @@
const LinearSolver::PerSolveOptions& per_solve_options,
double * x) {
EventLogger event_logger("SparseNormalCholeskySolver::SuiteSparse::Solve");
+ LinearSolver::Summary summary;
+ summary.termination_type = TOLERANCE;
+ summary.num_iterations = 1;
+ summary.status = "Success.";
const int num_cols = A->num_cols();
- LinearSolver::Summary summary;
Vector Atb = Vector::Zero(num_cols);
A->LeftMultiply(b, Atb.data());
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.
+ // Temporarily append a diagonal block to the A matrix, but undo
+ // it before returning the matrix to the user.
CompressedRowSparseMatrix D(per_solve_options.D, num_cols);
A->AppendRows(D);
}
VectorRef(x, num_cols).setZero();
-
cholmod_sparse lhs = ss_.CreateSparseMatrixTransposeView(A);
-
event_logger.AddEvent("Setup");
if (factor_ == NULL) {
if (options_.use_postordering) {
factor_ = ss_.BlockAnalyzeCholesky(&lhs,
A->col_blocks(),
- A->row_blocks());
+ A->row_blocks(),
+ &summary.status);
} else {
- factor_ = ss_.AnalyzeCholeskyWithNaturalOrdering(&lhs);
+ factor_ = ss_.AnalyzeCholeskyWithNaturalOrdering(&lhs, &summary.status);
}
}
event_logger.AddEvent("Analysis");
@@ -213,23 +219,20 @@
if (per_solve_options.D != NULL) {
A->DeleteRows(num_cols);
}
-
summary.termination_type = FATAL_ERROR;
return summary;
}
- const LinearSolverTerminationType status = ss_.Cholesky(&lhs, factor_);
- if (status != TOLERANCE) {
+ summary.termination_type = ss_.Cholesky(&lhs, factor_, &summary.status);
+ if (summary.termination_type != TOLERANCE) {
if (per_solve_options.D != NULL) {
A->DeleteRows(num_cols);
}
-
- summary.termination_type = FATAL_ERROR;
return summary;
}
cholmod_dense* rhs = ss_.CreateDenseVector(Atb.data(), num_cols, num_cols);
- cholmod_dense* sol = ss_.Solve(factor_, rhs);
+ cholmod_dense* sol = ss_.Solve(factor_, rhs, &summary.status);
event_logger.AddEvent("Solve");
ss_.Free(rhs);
@@ -237,14 +240,11 @@
A->DeleteRows(num_cols);
}
- summary.num_iterations = 1;
-
if (sol != NULL) {
memcpy(x, sol->x, num_cols * sizeof(*x));
-
ss_.Free(sol);
- sol = NULL;
- summary.termination_type = TOLERANCE;
+ } else {
+ summary.termination_type = FAILURE;
}
event_logger.AddEvent("Teardown");
diff --git a/internal/ceres/suitesparse.cc b/internal/ceres/suitesparse.cc
index ec93b1a..399b4ae 100644
--- a/internal/ceres/suitesparse.cc
+++ b/internal/ceres/suitesparse.cc
@@ -121,7 +121,8 @@
return v;
}
-cholmod_factor* SuiteSparse::AnalyzeCholesky(cholmod_sparse* A) {
+cholmod_factor* SuiteSparse::AnalyzeCholesky(cholmod_sparse* A,
+ string* status) {
// Cholmod can try multiple re-ordering strategies to find a fill
// reducing ordering. Here we just tell it use AMD with automatic
// matrix dependence choice of supernodal versus simplicial
@@ -131,34 +132,34 @@
cc_.supernodal = CHOLMOD_AUTO;
cholmod_factor* factor = cholmod_analyze(A, &cc_);
- if (cc_.status != CHOLMOD_OK) {
- LOG(ERROR) << "cholmod_analyze failed. error code: " << cc_.status;
- return NULL;
- }
-
- CHECK_NOTNULL(factor);
-
if (VLOG_IS_ON(2)) {
cholmod_print_common(const_cast<char*>("Symbolic Analysis"), &cc_);
}
- return factor;
+ if (cc_.status != CHOLMOD_OK) {
+ *status = StringPrintf("cholmod_analyze failed. error code: %d", cc_.status);
+ return NULL;
+ }
+
+ return CHECK_NOTNULL(factor);
}
cholmod_factor* SuiteSparse::BlockAnalyzeCholesky(
cholmod_sparse* A,
const vector<int>& row_blocks,
- const vector<int>& col_blocks) {
+ const vector<int>& col_blocks,
+ string* status) {
vector<int> ordering;
if (!BlockAMDOrdering(A, row_blocks, col_blocks, &ordering)) {
return NULL;
}
- return AnalyzeCholeskyWithUserOrdering(A, ordering);
+ return AnalyzeCholeskyWithUserOrdering(A, ordering, status);
}
cholmod_factor* SuiteSparse::AnalyzeCholeskyWithUserOrdering(
cholmod_sparse* A,
- const vector<int>& ordering) {
+ const vector<int>& ordering,
+ string* status) {
CHECK_EQ(ordering.size(), A->nrow);
cc_.nmethods = 1;
@@ -166,39 +167,34 @@
cholmod_factor* factor =
cholmod_analyze_p(A, const_cast<int*>(&ordering[0]), NULL, 0, &cc_);
- if (cc_.status != CHOLMOD_OK) {
- LOG(ERROR) << "cholmod_analyze failed. error code: " << cc_.status;
- return NULL;
- }
-
- CHECK_NOTNULL(factor);
-
if (VLOG_IS_ON(2)) {
cholmod_print_common(const_cast<char*>("Symbolic Analysis"), &cc_);
}
+ if (cc_.status != CHOLMOD_OK) {
+ *status = StringPrintf("cholmod_analyze failed. error code: %d", cc_.status);
+ return NULL;
+ }
- return factor;
+ return CHECK_NOTNULL(factor);
}
cholmod_factor* SuiteSparse::AnalyzeCholeskyWithNaturalOrdering(
- cholmod_sparse* A) {
+ cholmod_sparse* A,
+ string* status) {
cc_.nmethods = 1;
cc_.method[0].ordering = CHOLMOD_NATURAL;
cc_.postorder = 0;
cholmod_factor* factor = cholmod_analyze(A, &cc_);
- if (cc_.status != CHOLMOD_OK) {
- LOG(ERROR) << "cholmod_analyze failed. error code: " << cc_.status;
- return NULL;
- }
-
- CHECK_NOTNULL(factor);
-
if (VLOG_IS_ON(2)) {
cholmod_print_common(const_cast<char*>("Symbolic Analysis"), &cc_);
}
+ if (cc_.status != CHOLMOD_OK) {
+ *status = StringPrintf("cholmod_analyze failed. error code: %d", cc_.status);
+ return NULL;
+ }
- return factor;
+ return CHECK_NOTNULL(factor);
}
bool SuiteSparse::BlockAMDOrdering(const cholmod_sparse* A,
@@ -243,7 +239,9 @@
return true;
}
-LinearSolverTerminationType SuiteSparse::Cholesky(cholmod_sparse* A, cholmod_factor* L) {
+LinearSolverTerminationType SuiteSparse::Cholesky(cholmod_sparse* A,
+ cholmod_factor* L,
+ string* status) {
CHECK_NOTNULL(A);
CHECK_NOTNULL(L);
@@ -255,7 +253,7 @@
cc_.print = 0;
cc_.quick_return_if_not_posdef = 1;
- int status = cholmod_factorize(A, L, &cc_);
+ int cholmod_status = cholmod_factorize(A, L, &cc_);
cc_.print = old_print_level;
// TODO(sameeragarwal): This switch statement is not consistent. It
@@ -267,67 +265,73 @@
// (e.g. out of memory).
switch (cc_.status) {
case CHOLMOD_NOT_INSTALLED:
- LOG(WARNING) << "CHOLMOD failure: Method not installed.";
+ *status = "CHOLMOD failure: Method not installed.";
return FATAL_ERROR;
case CHOLMOD_OUT_OF_MEMORY:
- LOG(WARNING) << "CHOLMOD failure: Out of memory.";
+ *status = "CHOLMOD failure: Out of memory.";
return FATAL_ERROR;
case CHOLMOD_TOO_LARGE:
- LOG(WARNING) << "CHOLMOD failure: Integer overflow occured.";
+ *status = "CHOLMOD failure: Integer overflow occured.";
return FATAL_ERROR;
case CHOLMOD_INVALID:
- LOG(WARNING) << "CHOLMOD failure: Invalid input.";
+ *status = "CHOLMOD failure: Invalid input.";
return FATAL_ERROR;
case CHOLMOD_NOT_POSDEF:
- LOG(WARNING) << "CHOLMOD warning: Matrix not positive definite.";
+ *status = "CHOLMOD warning: Matrix not positive definite.";
return FAILURE;
case CHOLMOD_DSMALL:
- LOG(WARNING) << "CHOLMOD warning: D for LDL' or diag(L) or "
- << "LL' has tiny absolute value.";
+ *status = "CHOLMOD warning: D for LDL' or diag(L) or "
+ "LL' has tiny absolute value.";
return FAILURE;
case CHOLMOD_OK:
- if (status != 0) {
+ if (cholmod_status != 0) {
return TOLERANCE;
}
- LOG(WARNING) << "CHOLMOD failure: cholmod_factorize returned zero "
- << "but cholmod_common::status is CHOLMOD_OK."
- << "Please report this to ceres-solver@googlegroups.com.";
+
+ *status = "CHOLMOD failure: cholmod_factorize returned false "
+ "but cholmod_common::status is CHOLMOD_OK."
+ "Please report this to ceres-solver@googlegroups.com.";
return FATAL_ERROR;
default:
- LOG(WARNING) << "Unknown cholmod return code: " << cc_.status
- << ". Please report this to ceres-solver@googlegroups.com.";
+ *status =
+ StringPrintf("Unknown cholmod return code: %d. "
+ "Please report this to ceres-solver@googlegroups.com.",
+ cc_.status);
return FATAL_ERROR;
}
+
return FATAL_ERROR;
}
cholmod_dense* SuiteSparse::Solve(cholmod_factor* L,
- cholmod_dense* b) {
+ cholmod_dense* b,
+ string* status) {
if (cc_.status != CHOLMOD_OK) {
- LOG(WARNING) << "CHOLMOD status NOT OK";
+ *status = "cholmod_solve failed. CHOLMOD status is not CHOLMOD_OK";
return NULL;
}
return cholmod_solve(CHOLMOD_A, L, b, &cc_);
}
-void SuiteSparse::ApproximateMinimumDegreeOrdering(cholmod_sparse* matrix,
+bool SuiteSparse::ApproximateMinimumDegreeOrdering(cholmod_sparse* matrix,
int* ordering) {
- cholmod_amd(matrix, NULL, 0, ordering, &cc_);
+ return cholmod_amd(matrix, NULL, 0, ordering, &cc_);
}
-void SuiteSparse::ConstrainedApproximateMinimumDegreeOrdering(
+bool SuiteSparse::ConstrainedApproximateMinimumDegreeOrdering(
cholmod_sparse* matrix,
int* constraints,
int* ordering) {
#ifndef CERES_NO_CAMD
- cholmod_camd(matrix, NULL, 0, constraints, ordering, &cc_);
+ return cholmod_camd(matrix, NULL, 0, constraints, ordering, &cc_);
#else
LOG(FATAL) << "Congratulations you have found a bug in Ceres."
<< "Ceres Solver was compiled with SuiteSparse "
<< "version 4.1.0 or less. Calling this function "
<< "in that case is a bug. Please contact the"
<< "the Ceres Solver developers.";
+ return false;
#endif
}
diff --git a/internal/ceres/suitesparse.h b/internal/ceres/suitesparse.h
index c029b6d..0604654 100644
--- a/internal/ceres/suitesparse.h
+++ b/internal/ceres/suitesparse.h
@@ -139,12 +139,15 @@
// A is not modified, only the pattern of non-zeros of A is used,
// the actual numerical values in A are of no consequence.
//
+ // status contains an explanation of the failures if any.
+ //
// Caller owns the result.
- cholmod_factor* AnalyzeCholesky(cholmod_sparse* A);
+ cholmod_factor* AnalyzeCholesky(cholmod_sparse* A, string* status);
cholmod_factor* BlockAnalyzeCholesky(cholmod_sparse* A,
const vector<int>& row_blocks,
- const vector<int>& col_blocks);
+ const vector<int>& col_blocks,
+ string* status);
// If A is symmetric, then compute the symbolic Cholesky
// factorization of A(ordering, ordering). If A is unsymmetric, then
@@ -154,26 +157,38 @@
// A is not modified, only the pattern of non-zeros of A is used,
// the actual numerical values in A are of no consequence.
//
+ // status contains an explanation of the failures if any.
+ //
// Caller owns the result.
cholmod_factor* AnalyzeCholeskyWithUserOrdering(cholmod_sparse* A,
- const vector<int>& ordering);
+ const vector<int>& ordering,
+ string* status);
// Perform a symbolic factorization of A without re-ordering A. No
// postordering of the elimination tree is performed. This ensures
// that the symbolic factor does not introduce an extra permutation
// on the matrix. See the documentation for CHOLMOD for more details.
- cholmod_factor* AnalyzeCholeskyWithNaturalOrdering(cholmod_sparse* A);
+ //
+ // status contains an explanation of the failures if any.
+ cholmod_factor* AnalyzeCholeskyWithNaturalOrdering(cholmod_sparse* A,
+ string* status);
// Use the symbolic factorization in L, to find the numerical
// factorization for the matrix A or AA^T. Return true if
// successful, false otherwise. L contains the numeric factorization
// on return.
- LinearSolverTerminationType Cholesky(cholmod_sparse* A, cholmod_factor* L);
+ //
+ // status contains an explanation of the failures if any.
+ LinearSolverTerminationType Cholesky(cholmod_sparse* A,
+ cholmod_factor* L,
+ string* status);
// Given a Cholesky factorization of a matrix A = LL^T, solve the
// linear system Ax = b, and return the result. If the Solve fails
// NULL is returned. Caller owns the result.
- cholmod_dense* Solve(cholmod_factor* L, cholmod_dense* b);
+ //
+ // status contains an explanation of the failures if any.
+ cholmod_dense* Solve(cholmod_factor* L, cholmod_dense* b, string* solve);
// By virtue of the modeling layer in Ceres being block oriented,
// all the matrices used by Ceres are also block oriented. When
@@ -205,7 +220,7 @@
// Find a fill reducing approximate minimum degree
// ordering. ordering is expected to be large enough to hold the
// ordering.
- void ApproximateMinimumDegreeOrdering(cholmod_sparse* matrix, int* ordering);
+ bool ApproximateMinimumDegreeOrdering(cholmod_sparse* matrix, int* ordering);
// Before SuiteSparse version 4.2.0, cholmod_camd was only enabled
@@ -235,7 +250,7 @@
//
// If CERES_NO_CAMD is defined then calling this function will
// result in a crash.
- void ConstrainedApproximateMinimumDegreeOrdering(cholmod_sparse* matrix,
+ bool ConstrainedApproximateMinimumDegreeOrdering(cholmod_sparse* matrix,
int* constraints,
int* ordering);
diff --git a/internal/ceres/visibility_based_preconditioner.cc b/internal/ceres/visibility_based_preconditioner.cc
index aeaab12..e0737f2 100644
--- a/internal/ceres/visibility_based_preconditioner.cc
+++ b/internal/ceres/visibility_based_preconditioner.cc
@@ -435,14 +435,19 @@
// matrix contains the values.
lhs->stype = 1;
+ // TODO(sameeragarwal): Refactor to pipe this up and out.
+ string status;
+
// Symbolic factorization is computed if we don't already have one handy.
if (factor_ == NULL) {
- factor_ = ss_.BlockAnalyzeCholesky(lhs, block_size_, block_size_);
+ factor_ = ss_.BlockAnalyzeCholesky(lhs, block_size_, block_size_, &status);
}
- LinearSolverTerminationType status = ss_.Cholesky(lhs, factor_);
+ const LinearSolverTerminationType termination_type =
+ (factor_ != NULL) ? ss_.Cholesky(lhs, factor_, &status) : FATAL_ERROR;
+
ss_.Free(lhs);
- return status;
+ return termination_type;
}
void VisibilityBasedPreconditioner::RightMultiply(const double* x,
@@ -453,7 +458,9 @@
const int num_rows = m_->num_rows();
memcpy(CHECK_NOTNULL(tmp_rhs_)->x, x, m_->num_rows() * sizeof(*x));
- cholmod_dense* solution = CHECK_NOTNULL(ss->Solve(factor_, tmp_rhs_));
+ // TODO(sameeragarwal): Better error handling.
+ string status;
+ cholmod_dense* solution = CHECK_NOTNULL(ss->Solve(factor_, tmp_rhs_, &status));
memcpy(y, solution->x, sizeof(*y) * num_rows);
ss->Free(solution);
}
