diff --git a/docs/source/building.rst b/docs/source/building.rst
index 17e17d4..b91027b 100644
--- a/docs/source/building.rst
+++ b/docs/source/building.rst
@@ -72,7 +72,7 @@
 .. code-block:: bash
 
      # CMake
-     sudo apt-hey install cmake
+     sudo apt-get install cmake
      # gflags
      tar -xvzf gflags-2.0.tar.gz
      cd gflags-2.0
diff --git a/include/ceres/covariance.h b/include/ceres/covariance.h
index f65d9eb..b093026 100644
--- a/include/ceres/covariance.h
+++ b/include/ceres/covariance.h
@@ -51,7 +51,7 @@
 // 1. This is experimental code and the API WILL CHANGE before
 //    release.
 //
-// 2. WARNING: It is very easy to use this class incorrectly without
+// 2. It is very easy to use this class incorrectly without
 //    understanding the underlying mathematics. Please read and
 //    understand the documentation completely before attempting to use
 //    this class.
@@ -166,6 +166,21 @@
 // with indeterminacy. IEEE Transactions on Information Theory 47(5):
 // 2017-2028 (2001)
 //
+// Speed
+// -----
+//
+// When use_dense_linear_algebra = true, Eigen's JacobiSVD algorithm
+// is used to perform the computations. It is an accurate but slow
+// method and should only be used for small to moderate sized
+// problems.
+//
+// When use_dense_linear_algebra = false, SuiteSparse/CHOLMOD is used
+// to perform the computation. Recent versions of SuiteSparse (>=
+// 4.2.0) provide a much more efficient method for solving for rows of
+// the covariance matrix. Therefore, if you are doing large scale
+// covariance estimation, we strongly recommend using a recent version
+// of SuiteSparse.
+//
 // Example Usage
 // =============
 //
@@ -248,6 +263,27 @@
     // and by default Covariance::Compute will return false if it
     // encounters such a matrix.
     //
+    // use_dense_linear_algebra = false
+    // --------------------------------
+    //
+    // When performing large scale sparse covariance estimation,
+    // computing the exact value of the reciprocal condition number is
+    // not possible as it would require computing the eigenvalues of
+    // J'J.
+    //
+    // In this case we use cholmod_rcond, which uses the ratio of the
+    // smallest to the largest diagonal entries of the Cholesky
+    // factorization as an approximation to the reciprocal condition
+    // number.
+    //
+    // However, care must be taken as this is a heuristic and can
+    // sometimes be a very crude estimate. The default value of
+    // min_reciprocal_condition_number has been set to a conservative
+    // value, and sometimes the Covariance::Compute may return false
+    // even if it is possible to estimate the covariance reliably. In
+    // such cases, the user should exercise their judgement before
+    // lowering the value of min_reciprocal_condition_number.
+    //
     // use_dense_linear_algebra = true
     // -------------------------------
     //
diff --git a/internal/ceres/covariance_impl.cc b/internal/ceres/covariance_impl.cc
index 966ba35..978acc1 100644
--- a/internal/ceres/covariance_impl.cc
+++ b/internal/ceres/covariance_impl.cc
@@ -30,6 +30,10 @@
 
 #include "ceres/covariance_impl.h"
 
+#ifdef CERES_USE_OPENMP
+#include <omp.h>
+#endif
+
 #include <algorithm>
 #include <utility>
 #include <vector>
@@ -47,6 +51,36 @@
 
 namespace ceres {
 namespace internal {
+namespace {
+
+// Per thread storage for SuiteSparse.
+struct PerThreadContext {
+  PerThreadContext(int num_rows)
+      : solution(NULL),
+        solution_set(NULL),
+        y_workspace(NULL),
+        e_workspace(NULL),
+        rhs(NULL) {
+    rhs = ss.CreateDenseVector(NULL, num_rows, num_rows);
+  }
+
+  ~PerThreadContext() {
+    ss.Free(solution);
+    ss.Free(solution_set);
+    ss.Free(y_workspace);
+    ss.Free(e_workspace);
+    ss.Free(rhs);
+  }
+
+  cholmod_dense* solution;
+  cholmod_sparse* solution_set;
+  cholmod_dense* y_workspace;
+  cholmod_dense* e_workspace;
+  cholmod_dense* rhs;
+  SuiteSparse ss;
+};
+
+} // namespace
 
 typedef vector<pair<const double*, const double*> > CovarianceBlocks;
 
@@ -424,9 +458,6 @@
   const int* cols = covariance_matrix_->cols();
   double* values = covariance_matrix_->mutable_values();
 
-  cholmod_dense* rhs = ss_.CreateDenseVector(NULL, num_rows, num_rows);
-  double* rhs_x = reinterpret_cast<double*>(rhs->x);
-
   // The following loop exploits the fact that the i^th column of A^{-1}
   // is given by the solution to the linear system
   //
@@ -441,6 +472,9 @@
   // versions of SuiteSparse have the cholmod_solve2 function which
   // re-uses memory across calls.
 #if (SUITESPARSE_VERSION < 4002)
+  cholmod_dense* rhs = ss_.CreateDenseVector(NULL, num_rows, num_rows);
+  double* rhs_x = reinterpret_cast<double*>(rhs->x);
+
   for (int r = 0; r < num_rows; ++r) {
     int row_begin = rows[r];
     int row_end = rows[r + 1];
@@ -458,15 +492,35 @@
     ss_.Free(solution);
     rhs_x[r] = 0.0;
   }
-#else
-  // TODO(sameeragarwal) There should be a more efficient way
-  // involving the use of Bset but I am unable to make it work right
-  // now.
-  cholmod_dense* solution = NULL;
-  cholmod_sparse* solution_set = NULL;
-  cholmod_dense* y_workspace = NULL;
-  cholmod_dense* e_workspace = NULL;
 
+  ss_.Free(rhs);
+#else  // SUITESPARSE_VERSION < 4002
+
+  const int num_threads = options_.num_threads;
+  vector<PerThreadContext*> contexts(num_threads);
+  for (int i = 0; i < num_threads; ++i) {
+    contexts[i] = new PerThreadContext(num_rows);
+  }
+
+  // The first call to cholmod_solve2 is not thread safe, since it
+  // changes the factorization from supernodal to simplicial etc.
+  {
+    PerThreadContext* context = contexts[0];
+    double* context_rhs_x =  reinterpret_cast<double*>(context->rhs->x);
+    context_rhs_x[0] = 1.0;
+    cholmod_solve2(CHOLMOD_A,
+                   factor,
+                   context->rhs,
+                   NULL,
+                   &context->solution,
+                   &context->solution_set,
+                   &context->y_workspace,
+                   &context->e_workspace,
+                   context->ss.mutable_cc());
+    context_rhs_x[0] = 0.0;
+  }
+
+#pragma omp parallel for num_threads(num_threads) schedule(dynamic)
   for (int r = 0; r < num_rows; ++r) {
     int row_begin = rows[r];
     int row_end = rows[r + 1];
@@ -474,40 +528,52 @@
       continue;
     }
 
-    rhs_x[r] = 1.0;
+#  ifdef CERES_USE_OPENMP
+    int thread_id = omp_get_thread_num();
+#  else
+    int thread_id = 0;
+#  endif
 
+    PerThreadContext* context = contexts[thread_id];
+    double* context_rhs_x =  reinterpret_cast<double*>(context->rhs->x);
+    context_rhs_x[r] = 1.0;
+
+    // TODO(sameeragarwal) There should be a more efficient way
+    // involving the use of Bset but I am unable to make it work right
+    // now.
     cholmod_solve2(CHOLMOD_A,
                    factor,
-                   rhs,
+                   context->rhs,
                    NULL,
-                   &solution,
-                   &solution_set,
-                   &y_workspace,
-                   &e_workspace,
-                   ss_.mutable_cc());
+                   &context->solution,
+                   &context->solution_set,
+                   &context->y_workspace,
+                   &context->e_workspace,
+                   context->ss.mutable_cc());
 
-    double* solution_x = reinterpret_cast<double*>(solution->x);
+    double* solution_x = reinterpret_cast<double*>(context->solution->x);
     for (int idx = row_begin; idx < row_end; ++idx) {
       const int c = cols[idx];
       values[idx] = solution_x[c];
     }
-    rhs_x[r] = 0.0;
+    context_rhs_x[r] = 0.0;
   }
 
-  ss_.Free(solution);
-  ss_.Free(solution_set);
-  ss_.Free(y_workspace);
-  ss_.Free(e_workspace);
+  for (int i = 0; i < num_threads; ++i) {
+    delete contexts[i];
+  }
 
-#endif
+#endif  // SUITESPARSE_VERSION < 4002
 
-  ss_.Free(rhs);
   ss_.Free(factor);
   event_logger.AddEvent("Inversion");
   return true;
-#else
+
+#else  // CERES_NO_SUITESPARSE
+
   return false;
-#endif
+
+#endif  // CERES_NO_SUITESPARSE
 };
 
 bool CovarianceImpl::ComputeCovarianceValuesUsingEigen() {
@@ -549,19 +615,32 @@
   const int max_rank = min(num_singular_values,
                            num_singular_values - options_.null_space_rank);
 
+  // Compute the squared inverse of the singular values. Truncate the
+  // computation based on min_singular_value_ratio and
+  // null_space_rank. When either of these two quantities are active,
+  // the resulting covariance matrix is a Moore-Penrose inverse
+  // instead of a regular inverse.
   for (int i = 0; i < max_rank; ++i) {
     const double singular_value_ratio = singular_values[i] / max_singular_value;
-    if (singular_value_ratio >= min_singular_value_ratio) {
-      inverse_squared_singular_values[i] =
-          1.0 / (singular_values[i] * singular_values[i]);
-    } else if (!automatic_truncation) {
-      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;
-      return false;
+    if (singular_value_ratio < min_singular_value_ratio) {
+      // Since the singular values are in decreasing order, if
+      // automatic truncation is enabled, then from this point on
+      // all values will fail the ratio test and there is nothing to
+      // do in this loop.
+      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;
+        return false;
+      }
     }
+
+    inverse_squared_singular_values[i] =
+        1.0 / (singular_values[i] * singular_values[i]);
   }
 
   Matrix dense_covariance =
@@ -585,7 +664,5 @@
   return true;
 };
 
-
-
 }  // namespace internal
 }  // namespace ceres
diff --git a/internal/ceres/covariance_test.cc b/internal/ceres/covariance_test.cc
index 05fb234..bd0a0e2 100644
--- a/internal/ceres/covariance_test.cc
+++ b/internal/ceres/covariance_test.cc
@@ -405,6 +405,50 @@
   ComputeAndCompareCovarianceBlocks(options, expected_covariance);
 }
 
+#ifdef CERES_USE_OPENMP
+
+TEST_F(CovarianceTest, ThreadedNormalBehavior) {
+  // J
+  //
+  //   1  0  0  0  0  0
+  //   0  1  0  0  0  0
+  //   0  0  2  0  0  0
+  //   0  0  0  2  0  0
+  //   0  0  0  0  2  0
+  //   0  0  0  0  0  5
+  //  -5 -6  1  2  3  0
+  //   3 -2  0  0  0  2
+
+  // J'J
+  //
+  //   35  24 -5 -10 -15  6
+  //   24  41 -6 -12 -18 -4
+  //   -5  -6  5   2   3  0
+  //  -10 -12  2   8   6  0
+  //  -15 -18  3   6  13  0
+  //    6  -4  0   0   0 29
+
+  // inv(J'J) computed using octave.
+  double expected_covariance[] = {
+     7.0747e-02,  -8.4923e-03,   1.6821e-02,   3.3643e-02,   5.0464e-02,  -1.5809e-02,  // NOLINT
+    -8.4923e-03,   8.1352e-02,   2.4758e-02,   4.9517e-02,   7.4275e-02,   1.2978e-02,  // NOLINT
+     1.6821e-02,   2.4758e-02,   2.4904e-01,  -1.9271e-03,  -2.8906e-03,  -6.5325e-05,  // NOLINT
+     3.3643e-02,   4.9517e-02,  -1.9271e-03,   2.4615e-01,  -5.7813e-03,  -1.3065e-04,  // NOLINT
+     5.0464e-02,   7.4275e-02,  -2.8906e-03,  -5.7813e-03,   2.4133e-01,  -1.9598e-04,  // NOLINT
+    -1.5809e-02,   1.2978e-02,  -6.5325e-05,  -1.3065e-04,  -1.9598e-04,   3.9544e-02,  // NOLINT
+  };
+
+  Covariance::Options options;
+  options.num_threads = 4;
+
+  options.use_dense_linear_algebra = false;
+  ComputeAndCompareCovarianceBlocks(options, expected_covariance);
+
+  options.use_dense_linear_algebra = true;
+  ComputeAndCompareCovarianceBlocks(options, expected_covariance);
+}
+
+#endif  // CERES_USE_OPENMP
 
 TEST_F(CovarianceTest, ConstantParameterBlock) {
   problem_.SetParameterBlockConstant(parameters_);
@@ -644,5 +688,80 @@
   ComputeAndCompareCovarianceBlocks(options, expected_covariance);
 }
 
+class LargeScaleCovarianceTest : public ::testing::Test {
+ protected:
+  virtual void SetUp() {
+    num_parameter_blocks_ = 2000;
+    parameter_block_size_ = 5;
+    parameters_.reset(new double[parameter_block_size_ * num_parameter_blocks_]);
+
+    Matrix jacobian(parameter_block_size_, parameter_block_size_);
+    for (int i = 0; i < num_parameter_blocks_; ++i) {
+      jacobian.setIdentity();
+      jacobian *= (i + 1);
+
+      double* block_i = parameters_.get() + i * parameter_block_size_;
+      problem_.AddResidualBlock(new UnaryCostFunction(parameter_block_size_,
+                                                      parameter_block_size_,
+                                                      jacobian.data()),
+                                NULL,
+                                block_i );
+      for (int j = i; j < num_parameter_blocks_; ++j) {
+        double* block_j = parameters_.get() + j * parameter_block_size_;
+        all_covariance_blocks_.push_back(make_pair(block_i, block_j));
+      }
+    }
+  }
+
+  void ComputeAndCompare(int num_threads) {
+    Covariance::Options options;
+    options.use_dense_linear_algebra = false;
+    options.num_threads = num_threads;
+    Covariance covariance(options);
+    EXPECT_TRUE(covariance.Compute(all_covariance_blocks_, &problem_));
+
+    Matrix expected(parameter_block_size_, parameter_block_size_);
+    Matrix actual(parameter_block_size_, parameter_block_size_);
+    const double kTolerance = 1e-16;
+
+    for (int i = 0; i < num_parameter_blocks_; ++i) {
+      expected.setIdentity();
+      expected /= (i + 1.0) * (i + 1.0);
+
+      double* block_i = parameters_.get() + i * parameter_block_size_;
+      covariance.GetCovarianceBlock(block_i, block_i, actual.data());
+      EXPECT_NEAR((expected - actual).norm(), 0.0, kTolerance)
+          << "block: " << i << ", " << i << "\n"
+          << "expected: \n" << expected << "\n"
+          << "actual: \n" << actual;
+
+      expected.setZero();
+      for (int j = i + 1; j < num_parameter_blocks_; ++j) {
+        double* block_j = parameters_.get() + j * parameter_block_size_;
+        covariance.GetCovarianceBlock(block_i, block_j, actual.data());
+        EXPECT_NEAR((expected - actual).norm(), 0.0, kTolerance)
+            << "block: " << i << ", " << j << "\n"
+            << "expected: \n" << expected << "\n"
+            << "actual: \n" << actual;
+      }
+    }
+  }
+
+  scoped_array<double> parameters_;
+  int parameter_block_size_;
+  int num_parameter_blocks_;
+
+  Problem problem_;
+  vector<pair<const double*, const double*> > all_covariance_blocks_;
+};
+
+#if !defined(CERES_NO_SUITESPARSE) && defined(CERES_USE_OPENMP)
+
+TEST_F(LargeScaleCovarianceTest, Parallel) {
+  ComputeAndCompare(4);
+}
+
+#endif  // !defined(CERES_NO_SUITESPARSE) && defined(CERES_USE_OPENMP)
+
 }  // namespace internal
 }  // namespace ceres