commit b16e118b96c55451c0d8556f3c5b52ad36b69cac
author Sameer Agarwal <sameeragarwal@google.com> Mon Nov 25 05:47:43 2013 -0800
committer Sameer Agarwal <sameeragarwal@google.com> Tue Nov 26 10:00:03 2013 -0800
tree 2d5322d2f9fafde0a298e95356f70962acc57957
parent 5794d41be2d8d6a67dcdfe607e66050f0ac04c55

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

internal/ceres/sparse_normal_cholesky_solver.cc

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");