diff --git a/internal/ceres/eigen_vector_ops.h b/internal/ceres/eigen_vector_ops.h
index 525f2dd..b7f57f0 100644
--- a/internal/ceres/eigen_vector_ops.h
+++ b/internal/ceres/eigen_vector_ops.h
@@ -42,7 +42,10 @@
 // Blas1 operations on Eigen vectors. These functions are needed as an
 // abstraction layer so that we can use different versions of a vector style
 // object in the conjugate gradients linear solver.
-inline double Norm(const Vector& x, ContextImpl* context, int num_threads) {
+template <typename Derived>
+inline double Norm(const Eigen::DenseBase<Derived>& x,
+                   ContextImpl* context,
+                   int num_threads) {
   std::vector<double> norms(num_threads);
   ParallelFor(context,
               0,
diff --git a/internal/ceres/evaluation_benchmark.cc b/internal/ceres/evaluation_benchmark.cc
index 15b0534..c679885 100644
--- a/internal/ceres/evaluation_benchmark.cc
+++ b/internal/ceres/evaluation_benchmark.cc
@@ -43,6 +43,7 @@
 #include "ceres/evaluator.h"
 #include "ceres/implicit_schur_complement.h"
 #include "ceres/partitioned_matrix_view.h"
+#include "ceres/power_series_expansion_preconditioner.h"
 #include "ceres/preprocessor.h"
 #include "ceres/problem.h"
 #include "ceres/problem_impl.h"
@@ -327,6 +328,29 @@
   CHECK_GT(state_plus_delta.squaredNorm(), 0.);
 }
 
+static void PSEPreconditioner(benchmark::State& state,
+                              BALData* data,
+                              ContextImpl* context) {
+  LinearSolver::Options options;
+  options.num_threads = static_cast<int>(state.range(0));
+  options.elimination_groups.push_back(data->bal_problem->num_points());
+  options.context = context;
+
+  auto jacobian = data->ImplicitSchurComplementWithDiagonal(options);
+  Preconditioner::Options preconditioner_options(options);
+
+  PowerSeriesExpansionPreconditioner preconditioner(
+      jacobian, 10, 0, preconditioner_options);
+
+  Vector y = Vector::Zero(jacobian->num_cols());
+  Vector x = Vector::Random(jacobian->num_cols());
+
+  for (auto _ : state) {
+    preconditioner.RightMultiplyAndAccumulate(x.data(), y.data());
+  }
+  CHECK_GT(y.squaredNorm(), 0.);
+}
+
 static void PMVRightMultiplyAndAccumulateF(benchmark::State& state,
                                            BALData* data,
                                            ContextImpl* context) {
@@ -884,6 +908,16 @@
         ->Arg(8)
         ->Arg(16);
 
+    const std::string name_pse =
+        "PSEPreconditionerRightMultiplyAndAccumulate<" + path + ">";
+    ::benchmark::RegisterBenchmark(
+        name_pse.c_str(), ceres::internal::PSEPreconditioner, data, &context)
+        ->Arg(1)
+        ->Arg(2)
+        ->Arg(4)
+        ->Arg(8)
+        ->Arg(16);
+
     const std::string name_isc_no_diag =
         "ISCRightMultiplyAndAccumulate<" + path + ">";
     ::benchmark::RegisterBenchmark(name_isc_no_diag.c_str(),
diff --git a/internal/ceres/iterative_schur_complement_solver.cc b/internal/ceres/iterative_schur_complement_solver.cc
index c2e4b35..bcfb6e4 100644
--- a/internal/ceres/iterative_schur_complement_solver.cc
+++ b/internal/ceres/iterative_schur_complement_solver.cc
@@ -97,10 +97,13 @@
   reduced_linear_system_solution_.resize(schur_complement_->num_rows());
   reduced_linear_system_solution_.setZero();
   if (options_.use_spse_initialization) {
+    Preconditioner::Options preconditioner_options(options_);
+    preconditioner_options.type = SCHUR_POWER_SERIES_EXPANSION;
     PowerSeriesExpansionPreconditioner pse_solver(
         schur_complement_.get(),
         options_.max_num_spse_iterations,
-        options_.spse_tolerance);
+        options_.spse_tolerance,
+        preconditioner_options);
     pse_solver.RightMultiplyAndAccumulate(
         schur_complement_->rhs().data(),
         reduced_linear_system_solution_.data());
@@ -158,19 +161,8 @@
     return;
   }
 
-  Preconditioner::Options preconditioner_options;
-  preconditioner_options.type = options_.preconditioner_type;
-  preconditioner_options.visibility_clustering_type =
-      options_.visibility_clustering_type;
-  preconditioner_options.sparse_linear_algebra_library_type =
-      options_.sparse_linear_algebra_library_type;
-  preconditioner_options.num_threads = options_.num_threads;
-  preconditioner_options.row_block_size = options_.row_block_size;
-  preconditioner_options.e_block_size = options_.e_block_size;
-  preconditioner_options.f_block_size = options_.f_block_size;
-  preconditioner_options.elimination_groups = options_.elimination_groups;
+  Preconditioner::Options preconditioner_options(options_);
   CHECK(options_.context != nullptr);
-  preconditioner_options.context = options_.context;
 
   switch (options_.preconditioner_type) {
     case IDENTITY:
@@ -186,7 +178,10 @@
       // Ignoring the value of spse_tolerance to ensure preconditioner stays
       // fixed during the iterations of cg.
       preconditioner_ = std::make_unique<PowerSeriesExpansionPreconditioner>(
-          schur_complement_.get(), options_.max_num_spse_iterations, 0);
+          schur_complement_.get(),
+          options_.max_num_spse_iterations,
+          0,
+          preconditioner_options);
       break;
     case SCHUR_JACOBI:
       preconditioner_ = std::make_unique<SchurJacobiPreconditioner>(
diff --git a/internal/ceres/power_series_expansion_preconditioner.cc b/internal/ceres/power_series_expansion_preconditioner.cc
index 76de2fc..53fdc27 100644
--- a/internal/ceres/power_series_expansion_preconditioner.cc
+++ b/internal/ceres/power_series_expansion_preconditioner.cc
@@ -30,15 +30,20 @@
 
 #include "ceres/power_series_expansion_preconditioner.h"
 
+#include "ceres/eigen_vector_ops.h"
+#include "ceres/parallel_vector_ops.h"
+
 namespace ceres::internal {
 
 PowerSeriesExpansionPreconditioner::PowerSeriesExpansionPreconditioner(
     const ImplicitSchurComplement* isc,
     const int max_num_spse_iterations,
-    const double spse_tolerance)
+    const double spse_tolerance,
+    const Preconditioner::Options& options)
     : isc_(isc),
       max_num_spse_iterations_(max_num_spse_iterations),
-      spse_tolerance_(spse_tolerance) {}
+      spse_tolerance_(spse_tolerance),
+      options_(options) {}
 
 PowerSeriesExpansionPreconditioner::~PowerSeriesExpansionPreconditioner() =
     default;
@@ -53,17 +58,21 @@
   VectorRef yref(y, num_rows());
   Vector series_term(num_rows());
   Vector previous_series_term(num_rows());
-  yref.setZero();
-  isc_->block_diagonal_FtF_inverse()->RightMultiplyAndAccumulate(x, y);
-  previous_series_term = yref;
+  ParallelSetZero(options_.context, options_.num_threads, yref);
+  isc_->block_diagonal_FtF_inverse()->RightMultiplyAndAccumulate(
+      x, y, options_.context, options_.num_threads);
+  ParallelAssign(
+      options_.context, options_.num_threads, previous_series_term, yref);
 
-  const double norm_threshold = spse_tolerance_ * yref.norm();
+  const double norm_threshold =
+      spse_tolerance_ * Norm(yref, options_.context, options_.num_threads);
 
   for (int i = 1;; i++) {
-    series_term.setZero();
+    ParallelSetZero(options_.context, options_.num_threads, series_term);
     isc_->InversePowerSeriesOperatorRightMultiplyAccumulate(
         previous_series_term.data(), series_term.data());
-    yref += series_term;
+    ParallelAssign(
+        options_.context, options_.num_threads, yref, yref + series_term);
     if (i >= max_num_spse_iterations_ || series_term.norm() < norm_threshold) {
       break;
     }
diff --git a/internal/ceres/power_series_expansion_preconditioner.h b/internal/ceres/power_series_expansion_preconditioner.h
index 84e721e..9a993cf 100644
--- a/internal/ceres/power_series_expansion_preconditioner.h
+++ b/internal/ceres/power_series_expansion_preconditioner.h
@@ -41,20 +41,15 @@
 // This is a preconditioner via power series expansion of Schur
 // complement inverse based on "Weber et al, Power Bundle Adjustment for
 // Large-Scale 3D Reconstruction".
-//
-//
-// TODO(https://github.com/ceres-solver/ceres-solver/issues/934): In
-// PowerSeriesExpansionPreconditioner::RightMultiplyAndAccumulate only
-// operations performed via ImplicitSchurComplement are threaded.
-// PowerSeriesExpansionPreconditioner will benefit from multithreading
-// applied to remaning operations (block-sparse right product and several
-// vector operations)
 class CERES_NO_EXPORT PowerSeriesExpansionPreconditioner
     : public Preconditioner {
  public:
+  // TODO: Consider moving max_num_spse_iterations and spse_tolerance to
+  // Preconditioner::Options
   PowerSeriesExpansionPreconditioner(const ImplicitSchurComplement* isc,
                                      const int max_num_spse_iterations,
-                                     const double spse_tolerance);
+                                     const double spse_tolerance,
+                                     const Preconditioner::Options& options);
   PowerSeriesExpansionPreconditioner(
       const PowerSeriesExpansionPreconditioner&) = delete;
   void operator=(const PowerSeriesExpansionPreconditioner&) = delete;
@@ -68,6 +63,7 @@
   const ImplicitSchurComplement* isc_;
   const int max_num_spse_iterations_;
   const double spse_tolerance_;
+  const Preconditioner::Options options_;
 };
 
 }  // namespace ceres::internal
diff --git a/internal/ceres/power_series_expansion_preconditioner_test.cc b/internal/ceres/power_series_expansion_preconditioner_test.cc
index 6afca0b..1c04162 100644
--- a/internal/ceres/power_series_expansion_preconditioner_test.cc
+++ b/internal/ceres/power_series_expansion_preconditioner_test.cc
@@ -49,6 +49,7 @@
 
     options_.elimination_groups.push_back(problem_->num_eliminate_blocks);
     options_.preconditioner_type = SCHUR_POWER_SERIES_EXPANSION;
+    preconditioner_options_ = Preconditioner::Options(options_);
     isc_ = std::make_unique<ImplicitSchurComplement>(options_);
     isc_->Init(*A, D, problem_->b.get());
     num_f_cols_ = isc_->rhs().rows();
@@ -78,6 +79,7 @@
   int num_f_cols_;
   Matrix sc_inverse_expected_;
   LinearSolver::Options options_;
+  Preconditioner::Options preconditioner_options_;
 };
 
 TEST_F(PowerSeriesExpansionPreconditionerTest,
@@ -85,7 +87,7 @@
   const double spse_tolerance = kEpsilon;
   const int max_num_iterations = 50;
   PowerSeriesExpansionPreconditioner preconditioner(
-      isc_.get(), max_num_iterations, spse_tolerance);
+      isc_.get(), max_num_iterations, spse_tolerance, preconditioner_options_);
 
   Vector x(num_f_cols_), y(num_f_cols_);
   for (int i = 0; i < num_f_cols_; i++) {
@@ -108,7 +110,7 @@
   const double spse_tolerance = 0;
   const int max_num_iterations = 50;
   PowerSeriesExpansionPreconditioner preconditioner_fixed_n_iterations(
-      isc_.get(), max_num_iterations, spse_tolerance);
+      isc_.get(), max_num_iterations, spse_tolerance, preconditioner_options_);
 
   Vector x(num_f_cols_), y(num_f_cols_);
   for (int i = 0; i < num_f_cols_; i++) {
@@ -132,7 +134,7 @@
   const double spse_tolerance = 1 / kEpsilon;
   const int max_num_iterations = 50;
   PowerSeriesExpansionPreconditioner preconditioner_bad_tolerance(
-      isc_.get(), max_num_iterations, spse_tolerance);
+      isc_.get(), max_num_iterations, spse_tolerance, preconditioner_options_);
 
   Vector x(num_f_cols_), y(num_f_cols_);
   for (int i = 0; i < num_f_cols_; i++) {
@@ -153,7 +155,7 @@
   const double spse_tolerance = kEpsilon;
   const int max_num_iterations = 1;
   PowerSeriesExpansionPreconditioner preconditioner_bad_iterations_limit(
-      isc_.get(), max_num_iterations, spse_tolerance);
+      isc_.get(), max_num_iterations, spse_tolerance, preconditioner_options_);
 
   Vector x(num_f_cols_), y(num_f_cols_);
   for (int i = 0; i < num_f_cols_; i++) {
diff --git a/internal/ceres/preconditioner.h b/internal/ceres/preconditioner.h
index 2903d3b..07a7b1f 100644
--- a/internal/ceres/preconditioner.h
+++ b/internal/ceres/preconditioner.h
@@ -51,6 +51,20 @@
 class CERES_NO_EXPORT Preconditioner : public LinearOperator {
  public:
   struct Options {
+    Options() = default;
+    Options(const LinearSolver::Options& linear_solver_options)
+        : type(linear_solver_options.preconditioner_type),
+          visibility_clustering_type(
+              linear_solver_options.visibility_clustering_type),
+          sparse_linear_algebra_library_type(
+              linear_solver_options.sparse_linear_algebra_library_type),
+          num_threads(linear_solver_options.num_threads),
+          row_block_size(linear_solver_options.row_block_size),
+          e_block_size(linear_solver_options.e_block_size),
+          f_block_size(linear_solver_options.f_block_size),
+          elimination_groups(linear_solver_options.elimination_groups),
+          context(linear_solver_options.context) {}
+
     PreconditionerType type = JACOBI;
     VisibilityClusteringType visibility_clustering_type = CANONICAL_VIEWS;
     SparseLinearAlgebraLibraryType sparse_linear_algebra_library_type =