diff --git a/include/ceres/problem.h b/include/ceres/problem.h
index 7219699..909d1fc 100644
--- a/include/ceres/problem.h
+++ b/include/ceres/problem.h
@@ -543,9 +543,10 @@
       double* residuals,
       double** jacobians) const;
 
+  // Returns pointer to Problem implementation
+  internal::ProblemImpl* mutable_impl();
+
  private:
-  friend class Solver;
-  friend class Covariance;
   std::unique_ptr<internal::ProblemImpl> impl_;
 };
 
diff --git a/internal/ceres/CMakeLists.txt b/internal/ceres/CMakeLists.txt
index cde3895..fef1c2d 100644
--- a/internal/ceres/CMakeLists.txt
+++ b/internal/ceres/CMakeLists.txt
@@ -574,6 +574,9 @@
 endmacro()
 
 if (BUILD_BENCHMARKS)
+  add_executable(evaluation_benchmark evaluation_benchmark.cc)
+  add_dependencies_to_benchmark(evaluation_benchmark)
+
   add_executable(small_blas_gemv_benchmark small_blas_gemv_benchmark.cc)
   add_dependencies_to_benchmark(small_blas_gemv_benchmark)
 
diff --git a/internal/ceres/bundle_adjustment_test_util.h b/internal/ceres/bundle_adjustment_test_util.h
index 44281b6..8e6197d 100644
--- a/internal/ceres/bundle_adjustment_test_util.h
+++ b/internal/ceres/bundle_adjustment_test_util.h
@@ -65,6 +65,10 @@
 // problem is hard coded in the constructor.
 class BundleAdjustmentProblem {
  public:
+  BundleAdjustmentProblem(const std::string input_file) {
+    ReadData(input_file);
+    BuildProblem();
+  }
   BundleAdjustmentProblem() {
     const string input_file = TestFileAbsolutePath("problem-16-22106-pre.txt");
     ReadData(input_file);
diff --git a/internal/ceres/covariance.cc b/internal/ceres/covariance.cc
index d63dd37..029f3f5 100644
--- a/internal/ceres/covariance.cc
+++ b/internal/ceres/covariance.cc
@@ -51,12 +51,12 @@
 bool Covariance::Compute(
     const vector<pair<const double*, const double*>>& covariance_blocks,
     Problem* problem) {
-  return impl_->Compute(covariance_blocks, problem->impl_.get());
+  return impl_->Compute(covariance_blocks, problem->mutable_impl());
 }
 
 bool Covariance::Compute(const vector<const double*>& parameter_blocks,
                          Problem* problem) {
-  return impl_->Compute(parameter_blocks, problem->impl_.get());
+  return impl_->Compute(parameter_blocks, problem->mutable_impl());
 }
 
 bool Covariance::GetCovarianceBlock(const double* parameter_block1,
diff --git a/internal/ceres/evaluation_benchmark.cc b/internal/ceres/evaluation_benchmark.cc
new file mode 100644
index 0000000..64bb9a5
--- /dev/null
+++ b/internal/ceres/evaluation_benchmark.cc
@@ -0,0 +1,170 @@
+// Ceres Solver - A fast non-linear least squares minimizer
+// Copyright 2022 Google Inc. All rights reserved.
+// http://ceres-solver.org/
+//
+// Redistribution and use in source and binary forms, with or without
+// modification, are permitted provided that the following conditions are met:
+//
+// * Redistributions of source code must retain the above copyright notice,
+//   this list of conditions and the following disclaimer.
+// * Redistributions in binary form must reproduce the above copyright notice,
+//   this list of conditions and the following disclaimer in the documentation
+//   and/or other materials provided with the distribution.
+// * Neither the name of Google Inc. nor the names of its contributors may be
+//   used to endorse or promote products derived from this software without
+//   specific prior written permission.
+//
+// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
+// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+// POSSIBILITY OF SUCH DAMAGE.
+//
+// Authors: dmitriy.korchemkin@gmail.com (Dmitriy Korchemkin)
+
+#include <benchmark/benchmark.h>
+#include <gflags/gflags.h>
+
+#include "ceres/bundle_adjustment_test_util.h"
+#include "ceres/evaluator.h"
+#include "ceres/problem.h"
+#include "ceres/problem_impl.h"
+#include "ceres/program.h"
+#include "ceres/sparse_matrix.h"
+
+namespace ceres::internal {
+
+// Benchmark library might invoke benchmark function multiple times.
+// In order to save time required to parse BAL data, we ensure that
+// each dataset is being loaded at most once.
+struct BALData {
+  BALData(const std::string& path) {
+    bal_problem = std::make_unique<BundleAdjustmentProblem>(path);
+    CHECK(bal_problem != nullptr);
+
+    auto problem_impl = bal_problem->mutable_problem()->mutable_impl();
+    auto program = problem_impl->mutable_program();
+    parameters.resize(program->NumParameters());
+    program->ParameterBlocksToStateVector(parameters.data());
+  }
+
+  Vector parameters;
+  std::unique_ptr<BundleAdjustmentProblem> bal_problem;
+};
+
+static void Residuals(benchmark::State& state,
+                      BALData* data,
+                      ContextImpl* context) {
+  auto problem = data->bal_problem->mutable_problem();
+  auto problem_impl = problem->mutable_impl();
+  CHECK(problem_impl != nullptr);
+  const int num_threads = state.range(0);
+
+  Evaluator::Options options;
+  options.linear_solver_type = SPARSE_NORMAL_CHOLESKY;
+  options.num_threads = num_threads;
+  options.context = context;
+  options.num_eliminate_blocks = 0;
+
+  std::string error;
+  auto program = problem_impl->mutable_program();
+  auto evaluator = Evaluator::Create(options, program, &error);
+  CHECK(evaluator != nullptr);
+
+  double cost = 0.;
+  Vector residuals = Vector::Zero(program->NumResiduals());
+
+  Evaluator::EvaluateOptions eval_options;
+  for (auto _ : state) {
+    CHECK(evaluator->Evaluate(eval_options,
+                              data->parameters.data(),
+                              &cost,
+                              residuals.data(),
+                              nullptr,
+                              nullptr));
+  }
+}
+
+static void ResidualsAndJacobian(benchmark::State& state,
+                                 BALData* data,
+                                 ContextImpl* context) {
+  auto problem = data->bal_problem->mutable_problem();
+  auto problem_impl = problem->mutable_impl();
+  CHECK(problem_impl != nullptr);
+  const int num_threads = state.range(0);
+
+  Evaluator::Options options;
+  options.linear_solver_type = SPARSE_NORMAL_CHOLESKY;
+  options.num_threads = num_threads;
+  options.context = context;
+  options.num_eliminate_blocks = 0;
+
+  std::string error;
+  auto program = problem_impl->mutable_program();
+  auto evaluator = Evaluator::Create(options, program, &error);
+  CHECK(evaluator != nullptr);
+
+  double cost = 0.;
+  Vector residuals = Vector::Zero(program->NumResiduals());
+  auto jacobian = evaluator->CreateJacobian();
+
+  Evaluator::EvaluateOptions eval_options;
+  for (auto _ : state) {
+    CHECK(evaluator->Evaluate(eval_options,
+                              data->parameters.data(),
+                              &cost,
+                              residuals.data(),
+                              nullptr,
+                              jacobian.get()));
+  }
+}
+
+}  // namespace ceres::internal
+
+int main(int argc, char** argv) {
+  using namespace ceres::internal;
+  ::benchmark::Initialize(&argc, argv);
+
+  std::vector<std::unique_ptr<BALData>> benchmark_data;
+  if (argc == 1) {
+    LOG(FATAL) << "No input datasets specified. Usage: " << argv[0]
+               << " [benchmark flags] path_to_BAL_data_1.txt ... "
+                  "path_to_BAL_data_N.txt";
+    return -1;
+  }
+
+  ContextImpl context;
+  context.EnsureMinimumThreads(16);
+
+  for (int i = 1; i < argc; ++i) {
+    const std::string path(argv[i]);
+    const std::string name_residuals = "Residuals<" + path + ">";
+    benchmark_data.emplace_back(std::make_unique<BALData>(path));
+    auto data = benchmark_data.back().get();
+    ::benchmark::RegisterBenchmark(
+        name_residuals.c_str(), Residuals, data, &context)
+        ->Arg(1)
+        ->Arg(2)
+        ->Arg(4)
+        ->Arg(8)
+        ->Arg(16);
+    const std::string name_jacobians = "ResidualsAndJacobian<" + path + ">";
+    ::benchmark::RegisterBenchmark(
+        name_jacobians.c_str(), ResidualsAndJacobian, data, &context)
+        ->Arg(1)
+        ->Arg(2)
+        ->Arg(4)
+        ->Arg(8)
+        ->Arg(16);
+  }
+
+  ::benchmark::RunSpecifiedBenchmarks();
+  ::benchmark::Shutdown();
+  return 0;
+}
diff --git a/internal/ceres/problem.cc b/internal/ceres/problem.cc
index f4c2656..8c6d504 100644
--- a/internal/ceres/problem.cc
+++ b/internal/ceres/problem.cc
@@ -212,4 +212,6 @@
   impl_->GetResidualBlocksForParameterBlock(values, residual_blocks);
 }
 
+internal::ProblemImpl* Problem::mutable_impl() { return impl_.get(); }
+
 }  // namespace ceres
diff --git a/internal/ceres/solver.cc b/internal/ceres/solver.cc
index 07e7216..12dfc59 100644
--- a/internal/ceres/solver.cc
+++ b/internal/ceres/solver.cc
@@ -723,7 +723,7 @@
     return;
   }
 
-  ProblemImpl* problem_impl = problem->impl_.get();
+  ProblemImpl* problem_impl = problem->mutable_impl();
   Program* program = problem_impl->mutable_program();
   PreSolveSummarize(options, problem_impl, summary);
 
@@ -805,7 +805,7 @@
   }
 
   const double postprocessor_start_time = WallTimeInSeconds();
-  problem_impl = problem->impl_.get();
+  problem_impl = problem->mutable_impl();
   program = problem_impl->mutable_program();
   // On exit, ensure that the parameter blocks again point at the user
   // provided values and the parameter blocks are numbered according