diff --git a/examples/bundle_adjuster.cc b/examples/bundle_adjuster.cc
index 718fde9..1216a82 100644
--- a/examples/bundle_adjuster.cc
+++ b/examples/bundle_adjuster.cc
@@ -70,17 +70,20 @@
 DEFINE_string(preconditioner_type, "jacobi", "Options are: "
               "identity, jacobi, schur_jacobi, cluster_jacobi, "
               "cluster_tridiagonal");
+DEFINE_string(sparse_linear_algebra_library, "suitesparse",
+              "Options are: suitesparse and cxsparse");
 DEFINE_int32(num_iterations, 5, "Number of iterations");
 DEFINE_int32(num_threads, 1, "Number of threads");
 DEFINE_double(eta, 1e-2, "Default value for eta. Eta determines the "
              "accuracy of each linear solve of the truncated newton step. "
              "Changing this parameter can affect solve performance ");
-DEFINE_bool(use_schur_ordering, false, "Use automatic Schur ordering.");
+DEFINE_string(ordering_type, "schur", "Options are: schur, user, natural");
 DEFINE_bool(use_quaternions, false, "If true, uses quaternions to represent "
             "rotations. If false, angle axis is used");
 DEFINE_bool(use_local_parameterization, false, "For quaternions, use a local "
             "parameterization.");
 DEFINE_bool(robustify, false, "Use a robust loss function");
+DEFINE_bool(use_block_amd, true, "Use a block oriented fill reducing ordering.");
 
 namespace ceres {
 namespace examples {
@@ -138,10 +141,31 @@
     }
   }
 
+  if (FLAGS_sparse_linear_algebra_library == "suitesparse") {
+    options->sparse_linear_algebra_library = SUITE_SPARSE;
+  } else if (FLAGS_sparse_linear_algebra_library == "cxsparse") {
+    options->sparse_linear_algebra_library = CX_SPARSE;
+  } else {
+    LOG(FATAL) << "Unknown sparse linear algebra library type.";
+  }
+
   options->num_linear_solver_threads = FLAGS_num_threads;
 }
 
 void SetOrdering(BALProblem* bal_problem, Solver::Options* options) {
+  options->use_block_amd = FLAGS_use_block_amd;
+
+  // Only non-Schur solvers support the natural ordering for this
+  // problem.
+  if (FLAGS_ordering_type == "natural") {
+    if (options->linear_solver_type == SPARSE_SCHUR ||
+        options->linear_solver_type == DENSE_SCHUR ||
+        options->linear_solver_type == ITERATIVE_SCHUR) {
+      LOG(FATAL) << "Natural ordering with Schur type solver does not work.";
+    }
+    return;
+  }
+
   // Bundle adjustment problems have a sparsity structure that makes
   // them amenable to more specialized and much more efficient
   // solution strategies. The SPARSE_SCHUR, DENSE_SCHUR and
@@ -156,7 +180,7 @@
   // the right ParameterBlock ordering, or by manually specifying a
   // suitable ordering vector and defining
   // Options::num_eliminate_blocks.
-  if (FLAGS_use_schur_ordering) {
+  if (FLAGS_ordering_type == "schur") {
     options->ordering_type = ceres::SCHUR;
     return;
   }
@@ -272,7 +296,6 @@
   BuildProblem(&bal_problem, &problem);
   Solver::Options options;
   SetSolverOptionsFromFlags(&bal_problem, &options);
-
   Solver::Summary summary;
   Solve(options, &problem, &summary);
   std::cout << summary.FullReport() << "\n";
diff --git a/include/ceres/solver.h b/include/ceres/solver.h
index ed4c9b8..12351ab 100644
--- a/include/ceres/solver.h
+++ b/include/ceres/solver.h
@@ -83,6 +83,12 @@
       sparse_linear_algebra_library = CX_SPARSE;
 #endif
 
+#if defined(CERES_NO_SUITESPARSE)
+      use_block_amd = false;
+#else
+      use_block_amd = true;
+#endif
+
       preconditioner_type = JACOBI;
       num_linear_solver_threads = 1;
       num_eliminate_blocks = 0;
@@ -207,6 +213,19 @@
     // non-empty.
     vector<double*> ordering;
 
+    // By virtue of the modeling layer in Ceres being block oriented,
+    // all the matrices used by Ceres are also block oriented. When
+    // doing sparse direct factorization of these matrices (for
+    // SPARSE_NORMAL_CHOLESKY, SPARSE_SCHUR and ITERATIVE in
+    // conjunction with CLUSTER_TRIDIAGONAL AND CLUSTER_JACOBI
+    // preconditioners), the fill-reducing ordering algorithms can
+    // either be run on the block or the scalar form of these matrices.
+    // Running it on the block form exposes more of the super-nodal
+    // structure of the matrix to the factorization routines. Setting
+    // this parameter to true runs the ordering algorithms in block
+    // form. Currently this option only makes sense with
+    // sparse_linear_algebra_library = SUITE_SPARSE.
+    bool use_block_amd;
 
     // Minimum number of iterations for which the linear solver should
     // run, even if the convergence criterion is satisfied.
diff --git a/internal/ceres/CMakeLists.txt b/internal/ceres/CMakeLists.txt
index 2a0a6af..432346f 100644
--- a/internal/ceres/CMakeLists.txt
+++ b/internal/ceres/CMakeLists.txt
@@ -196,6 +196,7 @@
   ENDIF (GFLAGS)
   CERES_TEST(schur_ordering)
   CERES_TEST(solver_impl)
+  CERES_TEST(suitesparse)
   CERES_TEST(symmetric_linear_solver)
   CERES_TEST(triplet_sparse_matrix)
   CERES_TEST(trust_region_minimizer)
diff --git a/internal/ceres/compressed_row_jacobian_writer.cc b/internal/ceres/compressed_row_jacobian_writer.cc
index aa883b7..912c484 100644
--- a/internal/ceres/compressed_row_jacobian_writer.cc
+++ b/internal/ceres/compressed_row_jacobian_writer.cc
@@ -130,8 +130,24 @@
     }
     row_pos += num_residuals;
   }
-
   CHECK_EQ(num_jacobian_nonzeros, rows[total_num_residuals]);
+
+  // Populate the row and column block vectors for use by block
+  // oriented ordering algorithms. This is useful when
+  // Solver::Options::use_block_amd = true.
+  const vector<ParameterBlock*>& parameter_blocks = program_->parameter_blocks();
+  vector<int>& col_blocks = *(jacobian->mutable_col_blocks());
+  col_blocks.resize(parameter_blocks.size());
+  for (int i = 0; i <  parameter_blocks.size(); ++i) {
+    col_blocks[i] = parameter_blocks[i]->LocalSize();
+  }
+
+  vector<int>& row_blocks = *(jacobian->mutable_row_blocks());
+  row_blocks.resize(residual_blocks.size());
+  for (int i = 0; i <  residual_blocks.size(); ++i) {
+    row_blocks[i] = residual_blocks[i]->NumResiduals();
+  }
+
   return jacobian;
 }
 
diff --git a/internal/ceres/compressed_row_sparse_matrix.h b/internal/ceres/compressed_row_sparse_matrix.h
index c44f1db..7fb460a 100644
--- a/internal/ceres/compressed_row_sparse_matrix.h
+++ b/internal/ceres/compressed_row_sparse_matrix.h
@@ -31,11 +31,13 @@
 #ifndef CERES_INTERNAL_COMPRESSED_ROW_SPARSE_MATRIX_H_
 #define CERES_INTERNAL_COMPRESSED_ROW_SPARSE_MATRIX_H_
 
+#include <vector>
 #include <glog/logging.h>
 #include "ceres/sparse_matrix.h"
 #include "ceres/triplet_sparse_matrix.h"
 #include "ceres/internal/eigen.h"
 #include "ceres/internal/macros.h"
+#include "ceres/internal/port.h"
 #include "ceres/types.h"
 
 namespace ceres {
@@ -110,6 +112,12 @@
   const int* rows() const { return rows_.get(); }
   int* mutable_rows() { return rows_.get(); }
 
+  const vector<int>& row_blocks() const { return row_blocks_; };
+  vector<int>* mutable_row_blocks() { return &row_blocks_; };
+
+  const vector<int>& col_blocks() const { return col_blocks_; };
+  vector<int>* mutable_col_blocks() { return &col_blocks_; };
+
  private:
   scoped_array<int> cols_;
   scoped_array<int> rows_;
@@ -117,9 +125,16 @@
 
   int num_rows_;
   int num_cols_;
-
   int max_num_nonzeros_;
 
+  // If the matrix has an underlying block structure, then it can also
+  // carry with it row and column block sizes. This is auxilliary and
+  // optional information for use by algorithms operating on the
+  // matrix. The class itself does not make use of this information in
+  // any way.
+  vector<int> row_blocks_;
+  vector<int> col_blocks_;
+
   CERES_DISALLOW_COPY_AND_ASSIGN(CompressedRowSparseMatrix);
 };
 
diff --git a/internal/ceres/linear_solver.h b/internal/ceres/linear_solver.h
index 853abcb..31f8874 100644
--- a/internal/ceres/linear_solver.h
+++ b/internal/ceres/linear_solver.h
@@ -72,6 +72,7 @@
         : type(SPARSE_NORMAL_CHOLESKY),
           preconditioner_type(JACOBI),
           sparse_linear_algebra_library(SUITE_SPARSE),
+          use_block_amd(true),
           min_num_iterations(1),
           max_num_iterations(1),
           num_threads(1),
@@ -88,6 +89,9 @@
 
     SparseLinearAlgebraLibraryType sparse_linear_algebra_library;
 
+    // See solver.h for explanation of this option.
+    bool use_block_amd;
+
     // Number of internal iterations that the solver uses. This
     // parameter only makes sense for iterative solvers like CG.
     int min_num_iterations;
diff --git a/internal/ceres/schur_complement_solver.cc b/internal/ceres/schur_complement_solver.cc
index b30d6ed..5679a09 100644
--- a/internal/ceres/schur_complement_solver.cc
+++ b/internal/ceres/schur_complement_solver.cc
@@ -94,11 +94,11 @@
   const time_t backsubstitute_time = time(NULL);
   summary.termination_type = TOLERANCE;
 
-  VLOG(2) << "time (sec) total: " << backsubstitute_time - start_time
-          << " init: " << init_time - start_time
-          << " eliminate: " << eliminate_time - init_time
-          << " solve: " << solve_time - eliminate_time
-          << " backsubstitute: " << backsubstitute_time - solve_time;
+  VLOG(2) << "time (sec) total: " << (backsubstitute_time - start_time)
+          << " init: " << (init_time - start_time)
+          << " eliminate: " << (eliminate_time - init_time)
+          << " solve: " << (solve_time - eliminate_time)
+          << " backsubstitute: " << (backsubstitute_time - solve_time);
   return summary;
 }
 
@@ -150,15 +150,15 @@
     const LinearSolver::Options& options)
     : SchurComplementSolver(options) {
 #ifndef CERES_NO_SUITESPARSE
-  symbolic_factor_ = NULL;
+  factor_ = NULL;
 #endif  // CERES_NO_SUITESPARSE
 }
 
 SparseSchurComplementSolver::~SparseSchurComplementSolver() {
 #ifndef CERES_NO_SUITESPARSE
-  if (symbolic_factor_ != NULL) {
-    ss_.Free(symbolic_factor_);
-    symbolic_factor_ = NULL;
+  if (factor_ != NULL) {
+    ss_.Free(factor_);
+    factor_ = NULL;
   }
 #endif  // CERES_NO_SUITESPARSE
 }
@@ -171,13 +171,13 @@
   const int num_col_blocks = bs->cols.size();
   const int num_row_blocks = bs->rows.size();
 
-  vector<int> blocks(num_col_blocks - num_eliminate_blocks, 0);
+  blocks_.resize(num_col_blocks - num_eliminate_blocks, 0);
   for (int i = num_eliminate_blocks; i < num_col_blocks; ++i) {
-    blocks[i - num_eliminate_blocks] = bs->cols[i].size;
+    blocks_[i - num_eliminate_blocks] = bs->cols[i].size;
   }
 
   set<pair<int, int> > block_pairs;
-  for (int i = 0; i < blocks.size(); ++i) {
+  for (int i = 0; i < blocks_.size(); ++i) {
     block_pairs.insert(make_pair(i, i));
   }
 
@@ -230,7 +230,7 @@
     }
   }
 
-  set_lhs(new BlockRandomAccessSparseMatrix(blocks, block_pairs));
+  set_lhs(new BlockRandomAccessSparseMatrix(blocks_, block_pairs));
   set_rhs(new double[lhs()->num_rows()]);
 }
 
@@ -256,7 +256,8 @@
 // CHOLMOD's sparse cholesky factorization routines.
 bool SparseSchurComplementSolver::SolveReducedLinearSystemUsingSuiteSparse(
     double* solution) {
-  // Extract the TripletSparseMatrix that is used for actually storing S.
+  const time_t start_time = time(NULL);
+
   TripletSparseMatrix* tsm =
       const_cast<TripletSparseMatrix*>(
           down_cast<const BlockRandomAccessSparseMatrix*>(lhs())->matrix());
@@ -273,17 +274,32 @@
   // The matrix is symmetric, and the upper triangular part of the
   // matrix contains the values.
   cholmod_lhs->stype = 1;
+  const time_t lhs_time = time(NULL);
 
   cholmod_dense*  cholmod_rhs =
       ss_.CreateDenseVector(const_cast<double*>(rhs()), num_rows, num_rows);
+  const time_t rhs_time = time(NULL);
 
   // Symbolic factorization is computed if we don't already have one handy.
-  if (symbolic_factor_ == NULL) {
-    symbolic_factor_ = ss_.AnalyzeCholesky(cholmod_lhs);
+  if (factor_ == NULL) {
+    if (options().use_block_amd) {
+      factor_ = ss_.BlockAnalyzeCholesky(cholmod_lhs, blocks_, blocks_);
+    } else {
+      factor_ = ss_.AnalyzeCholesky(cholmod_lhs);
+    }
   }
 
+  if (VLOG_IS_ON(2)) {
+    cholmod_print_common("Symbolic Analysis", ss_.mutable_cc());
+  }
+
+  CHECK_NOTNULL(factor_);
+
+  const time_t symbolic_time = time(NULL);
   cholmod_dense* cholmod_solution =
-      ss_.SolveCholesky(cholmod_lhs, symbolic_factor_, cholmod_rhs);
+      ss_.SolveCholesky(cholmod_lhs, factor_, cholmod_rhs);
+
+  const time_t solve_time = time(NULL);
 
   ss_.Free(cholmod_lhs);
   cholmod_lhs = NULL;
@@ -298,6 +314,13 @@
   VectorRef(solution, num_rows)
       = VectorRef(static_cast<double*>(cholmod_solution->x), num_rows);
   ss_.Free(cholmod_solution);
+  const time_t final_time = time(NULL);
+  VLOG(2) << "time: " << (final_time - start_time)
+          << " lhs : " << (lhs_time - start_time)
+          << " rhs:  " << (rhs_time - lhs_time)
+          << " analyze: " <<  (symbolic_time - rhs_time)
+          << " factor_and_solve: " << (solve_time - symbolic_time)
+          << " cleanup: " << (final_time - solve_time);
   return true;
 }
 #else
diff --git a/internal/ceres/schur_complement_solver.h b/internal/ceres/schur_complement_solver.h
index 0a25ce1..a8f0852 100644
--- a/internal/ceres/schur_complement_solver.h
+++ b/internal/ceres/schur_complement_solver.h
@@ -31,6 +31,8 @@
 #ifndef CERES_INTERNAL_SCHUR_COMPLEMENT_SOLVER_H_
 #define CERES_INTERNAL_SCHUR_COMPLEMENT_SOLVER_H_
 
+#include <set>
+#include <utility>
 #include "ceres/block_random_access_matrix.h"
 #include "ceres/block_sparse_matrix.h"
 #include "ceres/block_structure.h"
@@ -158,11 +160,14 @@
   bool SolveReducedLinearSystemUsingSuiteSparse(double* solution);
   bool SolveReducedLinearSystemUsingCXSparse(double* solution);
 
+  // Size of the blocks in the Schur complement.
+  vector<int> blocks_;
+
 #ifndef CERES_NO_SUITESPARSE
   SuiteSparse ss_;
   // Symbolic factorization of the reduced linear system. Precomputed
   // once and reused in subsequent calls.
-  cholmod_factor* symbolic_factor_;
+  cholmod_factor* factor_;
 #endif  // CERES_NO_SUITESPARSE
   CERES_DISALLOW_COPY_AND_ASSIGN(SparseSchurComplementSolver);
 };
diff --git a/internal/ceres/solver_impl.cc b/internal/ceres/solver_impl.cc
index 4eb3d1f..fde5811 100644
--- a/internal/ceres/solver_impl.cc
+++ b/internal/ceres/solver_impl.cc
@@ -481,6 +481,7 @@
   linear_solver_options.preconditioner_type = options->preconditioner_type;
   linear_solver_options.sparse_linear_algebra_library =
       options->sparse_linear_algebra_library;
+  linear_solver_options.use_block_amd = options->use_block_amd;
 
 #ifdef CERES_NO_SUITESPARSE
   if (linear_solver_options.preconditioner_type == SCHUR_JACOBI) {
diff --git a/internal/ceres/sparse_normal_cholesky_solver.cc b/internal/ceres/sparse_normal_cholesky_solver.cc
index 5a082de..4191f3c 100644
--- a/internal/ceres/sparse_normal_cholesky_solver.cc
+++ b/internal/ceres/sparse_normal_cholesky_solver.cc
@@ -46,8 +46,6 @@
 #include "ceres/internal/scoped_ptr.h"
 #include "ceres/types.h"
 
-
-
 namespace ceres {
 namespace internal {
 
@@ -55,15 +53,15 @@
     const LinearSolver::Options& options)
     : options_(options) {
 #ifndef CERES_NO_SUITESPARSE
-  symbolic_factor_ = NULL;
+  factor_ = NULL;
 #endif
 }
 
 SparseNormalCholeskySolver::~SparseNormalCholeskySolver() {
 #ifndef CERES_NO_SUITESPARSE
-  if (symbolic_factor_ != NULL) {
-    ss_.Free(symbolic_factor_);
-    symbolic_factor_ = NULL;
+  if (factor_ != NULL) {
+    ss_.Free(factor_);
+    factor_ = NULL;
   }
 #endif
 }
@@ -183,13 +181,25 @@
   cholmod_dense* rhs = ss_.CreateDenseVector(Atb.data(), num_cols, num_cols);
   const time_t init_time = time(NULL);
 
-  if (symbolic_factor_ == NULL) {
-    symbolic_factor_ = CHECK_NOTNULL(ss_.AnalyzeCholesky(lhs.get()));
+  if (factor_ == NULL) {
+    if (options_.use_block_amd) {
+      factor_ = ss_.BlockAnalyzeCholesky(lhs.get(),
+                                         A->col_blocks(),
+                                         A->row_blocks());
+    } else {
+      factor_ = ss_.AnalyzeCholesky(lhs.get());
+    }
   }
 
+  if (VLOG_IS_ON(2)) {
+    cholmod_print_common("Symbolic Analysis", ss_.mutable_cc());
+  }
+
+  CHECK_NOTNULL(factor_);
+
   const time_t symbolic_time = time(NULL);
 
-  cholmod_dense* sol = ss_.SolveCholesky(lhs.get(), symbolic_factor_, rhs);
+  cholmod_dense* sol = ss_.SolveCholesky(lhs.get(), factor_, rhs);
   const time_t solve_time = time(NULL);
 
   ss_.Free(rhs);
@@ -209,11 +219,11 @@
   }
 
   const time_t cleanup_time = time(NULL);
-  VLOG(2) << "time (sec) total: " << cleanup_time - start_time
-          << " init: " << init_time - start_time
-          << " symbolic: " << symbolic_time - init_time
-          << " solve: " << solve_time - symbolic_time
-          << " cleanup: " << cleanup_time - solve_time;
+  VLOG(2) << "time (sec) total: " << (cleanup_time - start_time)
+          << " init: " << (init_time - start_time)
+          << " symbolic: " << (symbolic_time - init_time)
+          << " solve: " << (solve_time - symbolic_time)
+          << " cleanup: " << (cleanup_time - solve_time);
   return summary;
 }
 #else
diff --git a/internal/ceres/sparse_normal_cholesky_solver.h b/internal/ceres/sparse_normal_cholesky_solver.h
index e38f39b..6366f86 100644
--- a/internal/ceres/sparse_normal_cholesky_solver.h
+++ b/internal/ceres/sparse_normal_cholesky_solver.h
@@ -77,7 +77,7 @@
 #ifndef CERES_NO_SUITESPARSE
   SuiteSparse ss_;
   // Cached factorization
-  cholmod_factor* symbolic_factor_;
+  cholmod_factor* factor_;
 #endif  // CERES_NO_SUITESPARSE
 
 
diff --git a/internal/ceres/suitesparse.cc b/internal/ceres/suitesparse.cc
index 1cf6a74..c02d305 100644
--- a/internal/ceres/suitesparse.cc
+++ b/internal/ceres/suitesparse.cc
@@ -29,15 +29,14 @@
 // Author: sameeragarwal@google.com (Sameer Agarwal)
 
 #ifndef CERES_NO_SUITESPARSE
-
 #include "ceres/suitesparse.h"
 
+#include <vector>
 #include "cholmod.h"
 #include "ceres/compressed_row_sparse_matrix.h"
 #include "ceres/triplet_sparse_matrix.h"
 namespace ceres {
 namespace internal {
-
 cholmod_sparse* SuiteSparse::CreateSparseMatrix(TripletSparseMatrix* A) {
   cholmod_triplet triplet;
 
@@ -111,6 +110,13 @@
 }
 
 cholmod_factor* SuiteSparse::AnalyzeCholesky(cholmod_sparse* A) {
+  // 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
+  // factorization.
+  cc_.nmethods = 1;
+  cc_.method[0].ordering = CHOLMOD_AMD;
+  cc_.supernodal = CHOLMOD_AUTO;
   cholmod_factor* factor = cholmod_analyze(A, &cc_);
   CHECK_EQ(cc_.status, CHOLMOD_OK)
       << "Cholmod symbolic analysis failed " << cc_.status;
@@ -118,6 +124,146 @@
   return factor;
 }
 
+cholmod_factor* SuiteSparse::BlockAnalyzeCholesky(
+    cholmod_sparse* A,
+    const vector<int>& row_blocks,
+    const vector<int>& col_blocks) {
+  vector<int> ordering;
+  if (!BlockAMDOrdering(A, row_blocks, col_blocks, &ordering)) {
+    return NULL;
+  }
+  return AnalyzeCholeskyWithUserOrdering(A, ordering);
+}
+
+cholmod_factor* SuiteSparse::AnalyzeCholeskyWithUserOrdering(cholmod_sparse* A,
+                                                             const vector<int>& ordering) {
+  CHECK_EQ(ordering.size(), A->nrow);
+  cc_.nmethods = 1 ;
+  cc_.method[0].ordering = CHOLMOD_GIVEN;
+  cholmod_factor* factor  =
+      cholmod_analyze_p(A, const_cast<int*>(&ordering[0]), NULL, 0, &cc_);
+  CHECK_EQ(cc_.status, CHOLMOD_OK)
+      << "Cholmod symbolic analysis failed " << cc_.status;
+  CHECK_NOTNULL(factor);
+  return factor;
+}
+
+bool SuiteSparse::BlockAMDOrdering(const cholmod_sparse* A,
+                                   const vector<int>& row_blocks,
+                                   const vector<int>& col_blocks,
+                                   vector<int>* ordering) {
+  const int num_row_blocks = row_blocks.size();
+  const int num_col_blocks = col_blocks.size();
+
+  // Arrays storing the compressed column structure of the matrix
+  // incoding the block sparsity of A.
+  vector<int> block_cols;
+  vector<int> block_rows;
+
+  ScalarMatrixToBlockMatrix(A,
+                            row_blocks,
+                            col_blocks,
+                            &block_rows,
+                            &block_cols);
+
+  cholmod_sparse_struct block_matrix;
+  block_matrix.nrow = num_row_blocks;
+  block_matrix.ncol = num_col_blocks;
+  block_matrix.nzmax = block_rows.size();
+  block_matrix.p = reinterpret_cast<void*>(&block_cols[0]);
+  block_matrix.i = reinterpret_cast<void*>(&block_rows[0]);
+  block_matrix.x = NULL;
+  block_matrix.stype = A->stype;
+  block_matrix.itype = CHOLMOD_INT;
+  block_matrix.xtype = CHOLMOD_PATTERN;
+  block_matrix.dtype = CHOLMOD_DOUBLE;
+  block_matrix.sorted = 1;
+  block_matrix.packed = 1;
+
+  vector<int> block_ordering(num_row_blocks);
+  if (!cholmod_amd(&block_matrix, NULL, 0, &block_ordering[0], &cc_)) {
+    return false;
+  }
+
+  BlockOrderingToScalarOrdering(row_blocks, block_ordering, ordering);
+  return true;
+}
+
+void SuiteSparse::ScalarMatrixToBlockMatrix(const cholmod_sparse* A,
+                                            const vector<int>& row_blocks,
+                                            const vector<int>& col_blocks,
+                                            vector<int>* block_rows,
+                                            vector<int>* block_cols) {
+  CHECK_NOTNULL(block_rows)->clear();
+  CHECK_NOTNULL(block_cols)->clear();
+  const int num_row_blocks = row_blocks.size();
+  const int num_col_blocks = col_blocks.size();
+
+  vector<int> row_block_starts(num_row_blocks);
+  for (int i = 0, cursor = 0; i < num_row_blocks; ++i) {
+    row_block_starts[i] = cursor;
+    cursor += row_blocks[i];
+  }
+
+  // The reinterpret_cast is needed here because CHOLMOD stores arrays
+  // as void*.
+  const int* scalar_cols =  reinterpret_cast<const int*>(A->p);
+  const int* scalar_rows =  reinterpret_cast<const int*>(A->i);
+
+  // This loop extracts the block sparsity of the scalar sparse matrix
+  // A. It does so by iterating over the columns, but only considering
+  // the columns corresponding to the first element of each column
+  // block. Within each column, the inner loop iterates over the rows,
+  // and detects the presence of a row block by checking for the
+  // presence of a non-zero entry corresponding to its first element.
+  block_cols->push_back(0);
+  int c = 0;
+  for (int col_block = 0; col_block < num_col_blocks; ++col_block) {
+    int column_size = 0;
+    for (int idx = scalar_cols[c]; idx < scalar_cols[c + 1]; ++idx) {
+      vector<int>::const_iterator it = lower_bound(row_block_starts.begin(),
+                                                   row_block_starts.end(),
+                                                   scalar_rows[idx]);
+      DCHECK(it != row_block_starts.end());
+      // Only consider the first row of each row block.
+      if (*it != scalar_rows[idx]) {
+        continue;
+      }
+
+      block_rows->push_back(it - row_block_starts.begin());
+      ++column_size;
+    }
+    block_cols->push_back(block_cols->back() + column_size);
+    c += col_blocks[col_block];
+  }
+}
+
+void SuiteSparse::BlockOrderingToScalarOrdering(
+    const vector<int>& blocks,
+    const vector<int>& block_ordering,
+    vector<int>* scalar_ordering) {
+  CHECK_EQ(blocks.size(), block_ordering.size());
+  const int num_blocks = blocks.size();
+
+  // block_starts = [0, block1, block1 + block2 ..]
+  vector<int> block_starts(num_blocks);
+  for (int i = 0, cursor = 0; i < num_blocks ; ++i) {
+    block_starts[i] = cursor;
+    cursor += blocks[i];
+  }
+
+  scalar_ordering->resize(block_starts.back() + blocks.back());
+  int cursor = 0;
+  for (int i = 0; i < num_blocks; ++i) {
+    const int block_id = block_ordering[i];
+    const int block_size = blocks[block_id];
+    int block_position = block_starts[block_id];
+    for (int j = 0; j < block_size; ++j) {
+      (*scalar_ordering)[cursor++] = block_position++;
+    }
+  }
+}
+
 bool SuiteSparse::Cholesky(cholmod_sparse* A, cholmod_factor* L) {
   CHECK_NOTNULL(A);
   CHECK_NOTNULL(L);
diff --git a/internal/ceres/suitesparse.h b/internal/ceres/suitesparse.h
index 091e67a..eb691c0 100644
--- a/internal/ceres/suitesparse.h
+++ b/internal/ceres/suitesparse.h
@@ -37,6 +37,7 @@
 
 #include <cstring>
 #include <string>
+#include <vector>
 
 #include <glog/logging.h>
 #include "cholmod.h"
@@ -105,12 +106,35 @@
     cholmod_sdmult(A, 0, alpha_, beta_, x, y, &cc_);
   }
 
-  // Analyze the sparsity structure of the matrix A compute the
-  // symbolic factorization of A. A is not modified, only the pattern
-  // of non-zeros of A is used, the actual numerical values in A are
-  // of no consequence. Caller owns the result.
+  // Find an ordering of A or AA' (if A is unsymmetric) that minimizes
+  // the fill-in in the Cholesky factorization of the corresponding
+  // matrix. This is done by using the AMD algorithm.
+  //
+  // Using this ordering, the symbolic Cholesky factorization of A (or
+  // AA') is computed and returned.
+  //
+  // A is not modified, only the pattern of non-zeros of A is used,
+  // the actual numerical values in A are of no consequence.
+  //
+  // Caller owns the result.
   cholmod_factor* AnalyzeCholesky(cholmod_sparse* A);
 
+  cholmod_factor* BlockAnalyzeCholesky(cholmod_sparse* A,
+                                       const vector<int>& row_blocks,
+                                       const vector<int>& col_blocks);
+
+  // If A is symmetric, then compute the symbolic Cholesky
+  // factorization of A(ordering, ordering). If A is unsymmetric, then
+  // compute the symbolic factorization of
+  // A(ordering,:) A(ordering,:)'.
+  //
+  // A is not modified, only the pattern of non-zeros of A is used,
+  // the actual numerical values in A are of no consequence.
+  //
+  // Caller owns the result.
+  cholmod_factor* AnalyzeCholeskyWithUserOrdering(cholmod_sparse* A,
+                                                  const vector<int>& ordering);
+
   // 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
@@ -129,6 +153,56 @@
                                cholmod_factor* L,
                                cholmod_dense* b);
 
+  // By virtue of the modeling layer in Ceres being block oriented,
+  // all the matrices used by Ceres are also block oriented. When
+  // doing sparse direct factorization of these matrices the
+  // fill-reducing ordering algorithms (in particular AMD) can either
+  // be run on the block or the scalar form of these matrices. The two
+  // SuiteSparse::AnalyzeCholesky methods allows the the client to
+  // compute the symbolic factorization of a matrix by either using
+  // AMD on the matrix or a user provided ordering of the rows.
+  //
+  // But since the underlying matrices are block oriented, it is worth
+  // running AMD on just the block structre of these matrices and then
+  // lifting these block orderings to a full scalar ordering. This
+  // preserves the block structure of the permuted matrix, and exposes
+  // more of the super-nodal structure of the matrix to the numerical
+  // factorization routines.
+  //
+  // Find the block oriented AMD ordering of a matrix A, whose row and
+  // column blocks are given by row_blocks, and col_blocks
+  // respectively. The matrix may or may not be symmetric. The entries
+  // of col_blocks do not need to sum to the number of columns in
+  // A. If this is the case, only the first sum(col_blocks) are used
+  // to compute the ordering.
+  bool BlockAMDOrdering(const cholmod_sparse* A,
+                        const vector<int>& row_blocks,
+                        const vector<int>& col_blocks,
+                        vector<int>* ordering);
+
+  // Given a set of blocks and a permutation of these blocks, compute
+  // the corresponding "scalar" ordering, where the scalar ordering of
+  // size sum(blocks).
+  static void BlockOrderingToScalarOrdering(const vector<int>& blocks,
+                                            const vector<int>& block_ordering,
+                                            vector<int>* scalar_ordering);
+
+  // Extract the block sparsity pattern of the scalar sparse matrix
+  // A and return it in compressed column form. The compressed column
+  // form is stored in two vectors block_rows, and block_cols, which
+  // correspond to the row and column arrays in a compressed column sparse
+  // matrix.
+  //
+  // If c_ij is the block in the matrix A corresponding to row block i
+  // and column block j, then it is expected that A contains at least
+  // one non-zero entry corresponding to the top left entry of c_ij,
+  // as that entry is used to detect the presence of a non-zero c_ij.
+  static void ScalarMatrixToBlockMatrix(const cholmod_sparse* A,
+                                        const vector<int>& row_blocks,
+                                        const vector<int>& col_blocks,
+                                        vector<int>* block_rows,
+                                        vector<int>* block_cols);
+
   void Free(cholmod_sparse* m) { cholmod_free_sparse(&m, &cc_); }
   void Free(cholmod_dense* m)  { cholmod_free_dense(&m, &cc_);  }
   void Free(cholmod_factor* m) { cholmod_free_factor(&m, &cc_); }
diff --git a/internal/ceres/suitesparse_test.cc b/internal/ceres/suitesparse_test.cc
new file mode 100644
index 0000000..9211788
--- /dev/null
+++ b/internal/ceres/suitesparse_test.cc
@@ -0,0 +1,194 @@
+// Ceres Solver - A fast non-linear least squares minimizer
+// Copyright 2012 Google Inc. All rights reserved.
+// http://code.google.com/p/ceres-solver/
+//
+// 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 <algorithm>
+#include <glog/logging.h>
+#include "gtest/gtest.h"
+#include "ceres/suitesparse.h"
+#include "ceres/triplet_sparse_matrix.h"
+#include "ceres/internal/port.h"
+
+namespace ceres {
+namespace internal {
+
+TEST(SuiteSparse, BlockPermutationToScalarPermutation) {
+  vector<int> blocks;
+  //  Block structure
+  //  0  --1-  ---2---  ---3---  4
+  // [0, 1, 2, 3, 4, 5, 6, 7, 8, 9]
+  blocks.push_back(1);
+  blocks.push_back(2);
+  blocks.push_back(3);
+  blocks.push_back(3);
+  blocks.push_back(1);
+
+  // Block ordering
+  // [1, 0, 2, 4, 5]
+  vector<int> block_ordering;
+  block_ordering.push_back(1);
+  block_ordering.push_back(0);
+  block_ordering.push_back(2);
+  block_ordering.push_back(4);
+  block_ordering.push_back(3);
+
+  // Expected ordering
+  // [1, 2, 0, 3, 4, 5, 9, 6, 7, 8]
+  vector<int> expected_scalar_ordering;
+  expected_scalar_ordering.push_back(1);
+  expected_scalar_ordering.push_back(2);
+  expected_scalar_ordering.push_back(0);
+  expected_scalar_ordering.push_back(3);
+  expected_scalar_ordering.push_back(4);
+  expected_scalar_ordering.push_back(5);
+  expected_scalar_ordering.push_back(9);
+  expected_scalar_ordering.push_back(6);
+  expected_scalar_ordering.push_back(7);
+  expected_scalar_ordering.push_back(8);
+
+  vector<int> scalar_ordering;
+  SuiteSparse::BlockOrderingToScalarOrdering(blocks,
+                                             block_ordering,
+                                             &scalar_ordering);
+  EXPECT_EQ(scalar_ordering.size(), expected_scalar_ordering.size());
+  for (int i = 0; i < expected_scalar_ordering.size(); ++i) {
+    EXPECT_EQ(scalar_ordering[i], expected_scalar_ordering[i]);
+  }
+}
+
+// Helper function to fill the sparsity pattern of a TripletSparseMatrix.
+int FillBlock(const vector<int>& row_blocks,
+              const vector<int>& col_blocks,
+              const int row_block_id,
+              const int col_block_id,
+              int* rows,
+              int* cols) {
+  int row_pos = 0;
+  for (int i = 0; i < row_block_id; ++i) {
+    row_pos += row_blocks[i];
+  }
+
+  int col_pos = 0;
+  for (int i = 0; i < col_block_id; ++i) {
+    col_pos += col_blocks[i];
+  }
+
+  int offset = 0;
+  for (int r = 0; r < row_blocks[row_block_id]; ++r) {
+    for (int c = 0; c < col_blocks[col_block_id]; ++c, ++offset) {
+      rows[offset] = row_pos + r;
+      cols[offset] = col_pos + c;
+    }
+  }
+  return offset;
+}
+
+TEST(SuiteSparse, ScalarMatrixToBlockMatrix) {
+  // Block sparsity.
+  //
+  //     [1 2 3 2]
+  // [1]  x   x
+  // [2]    x   x
+  // [1]  x x
+  // num_nonzeros = 1 + 3 + 4 + 4 + 1 + 2 = 15
+
+  vector<int> col_blocks;
+  col_blocks.push_back(1);
+  col_blocks.push_back(2);
+  col_blocks.push_back(3);
+  col_blocks.push_back(2);
+
+  vector<int> row_blocks;
+  row_blocks.push_back(1);
+  row_blocks.push_back(2);
+  row_blocks.push_back(1);
+
+  TripletSparseMatrix tsm(4, 8, 15);
+  int* rows = tsm.mutable_rows();
+  int* cols = tsm.mutable_cols();
+  fill(tsm.mutable_values(), tsm.mutable_values() + 15, 1.0);
+  int offset = 0;
+
+#define CERES_TEST_FILL_BLOCK(r, c) \
+  offset += FillBlock(row_blocks, col_blocks, \
+                      row_block_id, col_block_id, \
+                      rows + offset, cols + offset);
+
+  CERES_TEST_FILL_BLOCK(0, 0);
+  CERES_TEST_FILL_BLOCK(2, 0);
+  CERES_TEST_FILL_BLOCK(1, 1);
+  CERES_TEST_FILL_BLOCK(2, 1);
+  CERES_TEST_FILL_BLOCK(0, 2);
+  CERES_TEST_FILL_BLOCK(1, 3);
+#undef CERES_TEST_FILL_BLOCK
+
+  tsm.set_num_nonzeros(offset);
+
+  SuiteSparse ss;
+  scoped_ptr<cholmod_sparse> ccsm(ss.CreateSparseMatrix(&tsm));
+
+  vector<int> expected_block_rows;
+  expected_block_rows.push_back(0);
+  expected_block_rows.push_back(2);
+  expected_block_rows.push_back(1);
+  expected_block_rows.push_back(2);
+  expected_block_rows.push_back(0);
+  expected_block_rows.push_back(1);
+
+  vector<int> expected_block_cols;
+  expected_block_cols.push_back(0);
+  expected_block_cols.push_back(2);
+  expected_block_cols.push_back(4);
+  expected_block_cols.push_back(5);
+  expected_block_cols.push_back(6);
+
+  vector<int> block_rows;
+  vector<int> block_cols;
+  SuiteSparse::ScalarMatrixToBlockMatrix(ccsm.get(),
+                                         row_blocks,
+                                         col_blocks,
+                                         &block_rows,
+                                         &block_cols);
+
+  EXPECT_EQ(block_cols.size(), expected_block_cols.size());
+  EXPECT_EQ(block_rows.size(), expected_block_rows.size());
+
+  for (int i = 0; i < expected_block_cols.size(); ++i) {
+    EXPECT_EQ(block_cols[i], expected_block_cols[i]);
+  }
+
+  for (int i = 0; i < expected_block_rows.size(); ++i) {
+    EXPECT_EQ(block_rows[i], expected_block_rows[i]);
+  }
+
+  ss.Free(ccsm.release());
+}
+
+}  // namespace internal
+}  // namespace ceres
diff --git a/internal/ceres/unsymmetric_linear_solver_test.cc b/internal/ceres/unsymmetric_linear_solver_test.cc
index a7ee23e..35cec80 100644
--- a/internal/ceres/unsymmetric_linear_solver_test.cc
+++ b/internal/ceres/unsymmetric_linear_solver_test.cc
@@ -62,6 +62,7 @@
     LinearSolver::Options options;
     options.type = linear_solver_type;
     options.sparse_linear_algebra_library = sparse_linear_algebra_library;
+    options.use_block_amd = false;
     scoped_ptr<LinearSolver> solver(LinearSolver::Create(options));
 
     LinearSolver::PerSolveOptions per_solve_options;
diff --git a/internal/ceres/visibility_based_preconditioner.cc b/internal/ceres/visibility_based_preconditioner.cc
index a930595..5a01006 100644
--- a/internal/ceres/visibility_based_preconditioner.cc
+++ b/internal/ceres/visibility_based_preconditioner.cc
@@ -33,7 +33,6 @@
 #include <algorithm>
 #include <functional>
 #include <iterator>
-#include <numeric>
 #include <set>
 #include <utility>
 #include <vector>
@@ -447,12 +446,21 @@
   // matrix contains the values.
   lhs->stype = 1;
 
-  // Symbolic factorization is computed if we don't already have one
-  // handy.
+  // Symbolic factorization is computed if we don't already have one handy.
   if (factor_ == NULL) {
-    factor_ = ss_.AnalyzeCholesky(lhs);
+    if (options_.use_block_amd) {
+      factor_ = ss_.BlockAnalyzeCholesky(lhs, block_size_, block_size_);
+    } else {
+      factor_ = ss_.AnalyzeCholesky(lhs);
+    }
   }
 
+  if (VLOG_IS_ON(2)) {
+    cholmod_print_common("Symbolic Analysis", ss_.mutable_cc());
+  }
+
+  CHECK_NOTNULL(factor_);
+
   bool status = ss_.Cholesky(lhs, factor_);
   ss_.Free(lhs);
   return status;