diff --git a/bazel/ceres.bzl b/bazel/ceres.bzl index a1117aaaa..dc92e0bb8 100644 --- a/bazel/ceres.bzl +++ b/bazel/ceres.bzl @@ -95,7 +95,8 @@ CERES_SRCS = ["internal/ceres/" + filename for filename in [ "manifold.cc", "minimizer.cc", "normal_prior.cc", - "parallel_for_cxx.cc", + "parallel_for.cc", + "parallel_invoke.cc", "parallel_utils.cc", "parameter_block_ordering.cc", "partitioned_matrix_view.cc", @@ -205,7 +206,6 @@ def ceres_library(name, "CERES_NO_EXPORT=", "CERES_NO_LAPACK", "CERES_NO_SUITESPARSE", - "CERES_USE_CXX_THREADS", "CERES_USE_EIGEN_SPARSE", ], includes = [ diff --git a/docs/source/features.rst b/docs/source/features.rst index d905415fb..ba22b9e3b 100644 --- a/docs/source/features.rst +++ b/docs/source/features.rst @@ -47,7 +47,7 @@ Why? linear system. To this end Ceres ships with a variety of linear solvers - dense QR and dense Cholesky factorization (using `Eigen`_, `LAPACK`_ or `CUDA`_) for dense problems, sparse - Cholesky factorization (`SuiteSparse`_, `Accelerate`_, `Eigen`_) + Cholesky factorization (`SuiteSparse`_, `Apple's Accelerate`_, `Eigen`_) for large sparse problems, custom Schur complement based dense, sparse, and iterative linear solvers for `bundle adjustment`_ problems. @@ -59,9 +59,8 @@ Why? of Non-linear Conjugate Gradients, BFGS and LBFGS. * **Speed** - Ceres Solver has been extensively optimized, with C++ - templating, hand written linear algebra routines and OpenMP or - modern C++ threads based multithreading of the Jacobian evaluation - and the linear solvers. + templating, hand written linear algebra routines and modern C++ threads + based multithreading of the Jacobian evaluation and the linear solvers. * **GPU Acceleration** If your system supports `CUDA`_ then Ceres Solver can use the Nvidia GPU on your system to speed up the solver. diff --git a/docs/source/installation.rst b/docs/source/installation.rst index d9e338fdf..3838ccb71 100644 --- a/docs/source/installation.rst +++ b/docs/source/installation.rst @@ -663,11 +663,6 @@ Options controlling Ceres configuration gains in the ``SPARSE_SCHUR`` solver, you can disable some of the template specializations by turning this ``OFF``. -#. ``CERES_THREADING_MODEL [Default: CXX_THREADS > OPENMP > NO_THREADS]``: - Multi-threading backend Ceres should be compiled with. This will - automatically be set to only accept the available subset of threading - options in the CMake GUI. - #. ``BUILD_SHARED_LIBS [Default: OFF]``: By default Ceres is built as a static library, turn this ``ON`` to instead build Ceres as a shared library. @@ -860,11 +855,6 @@ The Ceres components which can be specified are: #. ``SchurSpecializations``: Ceres built with Schur specializations (``SCHUR_SPECIALIZATIONS=ON``). -#. ``OpenMP``: Ceres built with OpenMP (``CERES_THREADING_MODEL=OPENMP``). - -#. ``Multithreading``: Ceres built with *a* multithreading library. - This is equivalent to (``CERES_THREAD != NO_THREADS``). - To specify one/multiple Ceres components use the ``COMPONENTS`` argument to `find_package() `_ like so: diff --git a/docs/source/nnls_solving.rst b/docs/source/nnls_solving.rst index e281b286b..7f7340bea 100644 --- a/docs/source/nnls_solving.rst +++ b/docs/source/nnls_solving.rst @@ -2298,9 +2298,7 @@ The three arrays will be: .. member:: int Solver::Summary::num_threads_used Number of threads actually used by the solver for Jacobian and - residual evaluation. This number is not equal to - :member:`Solver::Summary::num_threads_given` if none of `OpenMP` - or `CXX_THREADS` is available. + residual evaluation. .. member:: LinearSolverType Solver::Summary::linear_solver_type_given diff --git a/include/ceres/solver.h b/include/ceres/solver.h index 09c6525c9..52249c57d 100644 --- a/include/ceres/solver.h +++ b/include/ceres/solver.h @@ -991,8 +991,7 @@ class CERES_EXPORT Solver { int num_threads_given = -1; // Number of threads actually used by the solver for Jacobian and - // residual evaluation. This number is not equal to - // num_threads_given if OpenMP is not available. + // residual evaluation. int num_threads_used = -1; // Type of the linear solver requested by the user. diff --git a/internal/ceres/CMakeLists.txt b/internal/ceres/CMakeLists.txt index 23dbbd686..9729ef609 100644 --- a/internal/ceres/CMakeLists.txt +++ b/internal/ceres/CMakeLists.txt @@ -219,7 +219,8 @@ add_library(ceres_internal OBJECT linear_solver.cc low_rank_inverse_hessian.cc minimizer.cc - parallel_for_cxx.cc + parallel_for.cc + parallel_invoke.cc parallel_utils.cc parameter_block_ordering.cc partitioned_matrix_view.cc diff --git a/internal/ceres/block_random_access_dense_matrix.cc b/internal/ceres/block_random_access_dense_matrix.cc index deb0a80db..30d73f53f 100644 --- a/internal/ceres/block_random_access_dense_matrix.cc +++ b/internal/ceres/block_random_access_dense_matrix.cc @@ -34,7 +34,7 @@ #include #include "ceres/internal/eigen.h" -#include "ceres/parallel_for.h" +#include "ceres/parallel_vector_ops.h" #include "glog/logging.h" namespace ceres::internal { diff --git a/internal/ceres/block_random_access_diagonal_matrix.cc b/internal/ceres/block_random_access_diagonal_matrix.cc index fd109baf0..daabf8be1 100644 --- a/internal/ceres/block_random_access_diagonal_matrix.cc +++ b/internal/ceres/block_random_access_diagonal_matrix.cc @@ -40,6 +40,7 @@ #include "ceres/compressed_row_sparse_matrix.h" #include "ceres/internal/export.h" #include "ceres/parallel_for.h" +#include "ceres/parallel_vector_ops.h" #include "ceres/stl_util.h" #include "ceres/types.h" #include "glog/logging.h" diff --git a/internal/ceres/block_random_access_sparse_matrix.cc b/internal/ceres/block_random_access_sparse_matrix.cc index 36a318886..bd9f1baf3 100644 --- a/internal/ceres/block_random_access_sparse_matrix.cc +++ b/internal/ceres/block_random_access_sparse_matrix.cc @@ -37,7 +37,7 @@ #include #include "ceres/internal/export.h" -#include "ceres/parallel_for.h" +#include "ceres/parallel_vector_ops.h" #include "ceres/triplet_sparse_matrix.h" #include "ceres/types.h" #include "glog/logging.h" diff --git a/internal/ceres/block_sparse_matrix.cc b/internal/ceres/block_sparse_matrix.cc index f49f229ab..a60ee89cc 100644 --- a/internal/ceres/block_sparse_matrix.cc +++ b/internal/ceres/block_sparse_matrix.cc @@ -40,6 +40,7 @@ #include "ceres/crs_matrix.h" #include "ceres/internal/eigen.h" #include "ceres/parallel_for.h" +#include "ceres/parallel_vector_ops.h" #include "ceres/small_blas.h" #include "ceres/triplet_sparse_matrix.h" #include "glog/logging.h" diff --git a/internal/ceres/eigen_vector_ops.h b/internal/ceres/eigen_vector_ops.h index 91f294603..8b519e6d6 100644 --- a/internal/ceres/eigen_vector_ops.h +++ b/internal/ceres/eigen_vector_ops.h @@ -35,6 +35,7 @@ #include "ceres/internal/eigen.h" #include "ceres/parallel_for.h" +#include "ceres/parallel_vector_ops.h" namespace ceres::internal { diff --git a/internal/ceres/implicit_schur_complement.cc b/internal/ceres/implicit_schur_complement.cc index dbba3a759..9c73de689 100644 --- a/internal/ceres/implicit_schur_complement.cc +++ b/internal/ceres/implicit_schur_complement.cc @@ -36,6 +36,7 @@ #include "ceres/internal/eigen.h" #include "ceres/linear_solver.h" #include "ceres/parallel_for.h" +#include "ceres/parallel_vector_ops.h" #include "ceres/types.h" #include "glog/logging.h" diff --git a/internal/ceres/levenberg_marquardt_strategy.cc b/internal/ceres/levenberg_marquardt_strategy.cc index 81ced2cf0..b37f74e7e 100644 --- a/internal/ceres/levenberg_marquardt_strategy.cc +++ b/internal/ceres/levenberg_marquardt_strategy.cc @@ -38,7 +38,7 @@ #include "ceres/internal/eigen.h" #include "ceres/linear_least_squares_problems.h" #include "ceres/linear_solver.h" -#include "ceres/parallel_for.h" +#include "ceres/parallel_vector_ops.h" #include "ceres/sparse_matrix.h" #include "ceres/trust_region_strategy.h" #include "ceres/types.h" diff --git a/internal/ceres/parallel_for.cc b/internal/ceres/parallel_for.cc new file mode 100644 index 000000000..a92276d4f --- /dev/null +++ b/internal/ceres/parallel_for.cc @@ -0,0 +1,63 @@ +// Ceres Solver - A fast non-linear least squares minimizer +// Copyright 2018 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. +// +// Author: vitus@google.com (Michael Vitus) + +#include "ceres/parallel_for.h" + +#include +#include +#include +#include +#include +#include +#include + +#include "ceres/internal/config.h" +#include "ceres/parallel_vector_ops.h" +#include "glog/logging.h" + +namespace ceres::internal { + +int MaxNumThreadsAvailable() { return ThreadPool::MaxNumThreadsAvailable(); } + +void ParallelSetZero(ContextImpl* context, + int num_threads, + double* values, + int num_values) { + ParallelFor(context, + 0, + num_values, + num_threads, + [values](std::tuple range) { + auto [start, end] = range; + std::fill(values + start, values + end, 0.); + }); +} + +} // namespace ceres::internal diff --git a/internal/ceres/parallel_for.h b/internal/ceres/parallel_for.h index bf38603fa..7536c74c3 100644 --- a/internal/ceres/parallel_for.h +++ b/internal/ceres/parallel_for.h @@ -32,14 +32,14 @@ #ifndef CERES_INTERNAL_PARALLEL_FOR_H_ #define CERES_INTERNAL_PARALLEL_FOR_H_ -#include -#include #include +#include #include "ceres/context_impl.h" -#include "ceres/internal/disable_warnings.h" #include "ceres/internal/eigen.h" #include "ceres/internal/export.h" +#include "ceres/parallel_invoke.h" +#include "ceres/partition_range_for_parallel_for.h" #include "glog/logging.h" namespace ceres::internal { @@ -51,237 +51,82 @@ inline decltype(auto) MakeConditionalLock(const int num_threads, : std::unique_lock{m}; } -// Returns the maximum number of threads supported by the threading backend -// Ceres was compiled with. -CERES_NO_EXPORT -int MaxNumThreadsAvailable(); - -// Parallel for implementations share a common set of routines in order -// to enforce inlining of loop bodies, ensuring that single-threaded -// performance is equivalent to a simple for loop -namespace parallel_for_details { -// Get arguments of callable as a tuple -template -std::tuple...> args_of(void (F::*)(Args...) const); - -template -using args_of_t = decltype(args_of(&F::operator())); - -// Parallelizable functions might require passing thread_id as the first -// argument. This class supplies thread_id argument to functions that -// support it and ignores it otherwise. -template -struct InvokeImpl; - -// For parallel for iterations of type [](int i) -> void -template -struct InvokeImpl> { - static void InvokeOnSegment(int thread_id, - int start, - int end, - const F& function) { - (void)thread_id; - for (int i = start; i < end; ++i) { - function(i); - } - } -}; - -// For parallel for iterations of type [](int thread_id, int i) -> void -template -struct InvokeImpl> { - static void InvokeOnSegment(int thread_id, - int start, - int end, - const F& function) { - for (int i = start; i < end; ++i) { - function(thread_id, i); - } - } -}; - -// For parallel for iterations of type [](tuple range) -> void -template -struct InvokeImpl>> { - static void InvokeOnSegment(int thread_id, - int start, - int end, - const F& function) { - (void)thread_id; - function(std::make_tuple(start, end)); - } -}; - -// For parallel for iterations of type [](int thread_id, tuple range) -// -> void -template -struct InvokeImpl>> { - static void InvokeOnSegment(int thread_id, - int start, - int end, - const F& function) { - function(thread_id, std::make_tuple(start, end)); - } -}; - -// Invoke function passing thread_id only if required -template -void InvokeOnSegment(int thread_id, int start, int end, const F& function) { - InvokeImpl>::InvokeOnSegment(thread_id, start, end, function); -} - -// Check if it is possible to split range [start; end) into at most -// max_num_partitions contiguous partitions of cost not greater than -// max_partition_cost. Inclusive integer cumulative costs are provided by -// cumulative_cost_data objects, with cumulative_cost_offset being a total cost -// of all indices (starting from zero) preceding start element. Cumulative costs -// are returned by cumulative_cost_fun called with a reference to -// cumulative_cost_data element with index from range[start; end), and should be -// non-decreasing. Partition of the range is returned via partition argument -template -bool MaxPartitionCostIsFeasible(int start, - int end, - int max_num_partitions, - int max_partition_cost, - int cumulative_cost_offset, - const CumulativeCostData* cumulative_cost_data, - const CumulativeCostFun& cumulative_cost_fun, - std::vector* partition) { - partition->clear(); - partition->push_back(start); - int partition_start = start; - int cost_offset = cumulative_cost_offset; - - while (partition_start < end) { - // Already have max_num_partitions - if (partition->size() > max_num_partitions) { - return false; - } - const int target = max_partition_cost + cost_offset; - const int partition_end = - std::partition_point( - cumulative_cost_data + partition_start, - cumulative_cost_data + end, - [&cumulative_cost_fun, target](const CumulativeCostData& item) { - return cumulative_cost_fun(item) <= target; - }) - - cumulative_cost_data; - // Unable to make a partition from a single element - if (partition_end == partition_start) { - return false; - } - - const int cost_last = - cumulative_cost_fun(cumulative_cost_data[partition_end - 1]); - partition->push_back(partition_end); - partition_start = partition_end; - cost_offset = cost_last; - } - return true; -} - -// Split integer interval [start, end) into at most max_num_partitions -// contiguous intervals, minimizing maximal total cost of a single interval. -// Inclusive integer cumulative costs for each (zero-based) index are provided -// by cumulative_cost_data objects, and are returned by cumulative_cost_fun call -// with a reference to one of the objects from range [start, end) -template -std::vector ComputePartition( - int start, - int end, - int max_num_partitions, - const CumulativeCostData* cumulative_cost_data, - const CumulativeCostFun& cumulative_cost_fun) { - // Given maximal partition cost, it is possible to verify if it is admissible - // and obtain corresponding partition using MaxPartitionCostIsFeasible - // function. In order to find the lowest admissible value, a binary search - // over all potentially optimal cost values is being performed - const int cumulative_cost_last = - cumulative_cost_fun(cumulative_cost_data[end - 1]); - const int cumulative_cost_offset = - start ? cumulative_cost_fun(cumulative_cost_data[start - 1]) : 0; - const int total_cost = cumulative_cost_last - cumulative_cost_offset; - - // Minimal maximal partition cost is not smaller than the average - // We will use non-inclusive lower bound - int partition_cost_lower_bound = total_cost / max_num_partitions - 1; - // Minimal maximal partition cost is not larger than the total cost - // Upper bound is inclusive - int partition_cost_upper_bound = total_cost; - - std::vector partition, partition_upper_bound; - // Binary search over partition cost, returning the lowest admissible cost - while (partition_cost_upper_bound - partition_cost_lower_bound > 1) { - partition.reserve(max_num_partitions + 1); - const int partition_cost = - partition_cost_lower_bound + - (partition_cost_upper_bound - partition_cost_lower_bound) / 2; - bool admissible = MaxPartitionCostIsFeasible(start, - end, - max_num_partitions, - partition_cost, - cumulative_cost_offset, - cumulative_cost_data, - cumulative_cost_fun, - &partition); - if (admissible) { - partition_cost_upper_bound = partition_cost; - std::swap(partition, partition_upper_bound); - } else { - partition_cost_lower_bound = partition_cost; - } - } - - // After binary search over partition cost, interval - // (partition_cost_lower_bound, partition_cost_upper_bound] contains the only - // admissible partition cost value - partition_cost_upper_bound - // - // Partition for this cost value might have been already computed - if (partition_upper_bound.empty() == false) { - return partition_upper_bound; - } - // Partition for upper bound is not computed if and only if upper bound was - // never updated This is a simple case of a single interval containing all - // values, which we were not able to break into pieces - partition = {start, end}; - return partition; -} -} // namespace parallel_for_details - -// Forward declaration of parallel invocation function that is to be -// implemented by each threading backend -template -void ParallelInvoke(ContextImpl* context, - int start, - int end, - int num_threads, - const F& function); +// Returns the maximum supported number of threads +CERES_NO_EXPORT int MaxNumThreadsAvailable(); // Execute the function for every element in the range [start, end) with at most // num_threads. It will execute all the work on the calling thread if // num_threads or (end - start) is equal to 1. +// Depending on function signature, it will be supplied with either loop index +// or a range of loop indicies; function can also be supplied with thread_id. +// The following function signatures are supported: +// - Functions accepting a single loop index: +// - [](int index) { ... } +// - [](int thread_id, int index) { ... } +// - Functions accepting a range of loop index: +// - [](std::tuple index) { ... } +// - [](int thread_id, std::tuple index) { ... } // -// Functions with two arguments will be passed thread_id and loop index on each -// invocation, functions with one argument will be invoked with loop index +// When distributing workload between threads, it is assumed that each loop +// iteration takes approximately equal time to complete. template -void ParallelFor(ContextImpl* context, - int start, - int end, - int num_threads, - const F& function) { - using namespace parallel_for_details; +void ParallelFor( + ContextImpl* context, int start, int end, int num_threads, F&& function) { CHECK_GT(num_threads, 0); if (start >= end) { return; } if (num_threads == 1 || end - start == 1) { - InvokeOnSegment(0, start, end, function); + InvokeOnSegment(0, std::make_tuple(start, end), std::forward(function)); return; } CHECK(context != nullptr); - ParallelInvoke(context, start, end, num_threads, function); + ParallelInvoke(context, start, end, num_threads, std::forward(function)); +} + +// Execute function for every element in the range [start, end) with at most +// num_threads, using user-provided partitions array. +// When distributing workload between threads, it is assumed that each segment +// bounded by adjacent elements of partitions array takes approximately equal +// time to process. +template +void ParallelFor(ContextImpl* context, + int start, + int end, + int num_threads, + F&& function, + const std::vector& partitions) { + CHECK_GT(num_threads, 0); + if (start >= end) { + return; + } + CHECK_EQ(partitions.front(), start); + CHECK_EQ(partitions.back(), end); + if (num_threads == 1 || end - start <= num_threads) { + ParallelFor(context, start, end, num_threads, std::forward(function)); + return; + } + CHECK_GT(partitions.size(), 1); + const int num_partitions = partitions.size() - 1; + ParallelFor(context, + 0, + num_partitions, + num_threads, + [&function, &partitions](int thread_id, + std::tuple partition_ids) { + // partition_ids is a range of partition indices + const auto [partition_start, partition_end] = partition_ids; + // Execution over several adjacent segments is equivalent + // to execution over union of those segments (which is also a + // contiguous segment) + const int range_start = partitions[partition_start]; + const int range_end = partitions[partition_end]; + // Range of original loop indices + const auto range = std::make_tuple(range_start, range_end); + InvokeOnSegment(thread_id, range, function); + }); } // Execute function for every element in the range [start, end) with at most @@ -293,123 +138,43 @@ void ParallelFor(ContextImpl* context, // user-defined object. Only indices from [start, end) will be referenced. This // routine assumes that cumulative_cost_fun is non-decreasing (in other words, // all costs are non-negative); +// When distributing workload between threads, input range of loop indices will +// be partitioned into disjoint contiguous intervals, with the maximal cost +// being minimized. +// For example, with iteration costs of [1, 1, 5, 3, 1, 4] cumulative_cost_fun +// should return [1, 2, 7, 10, 11, 15], and with num_threads = 4 this range +// will be split into segments [0, 2) [2, 3) [3, 5) [5, 6) with costs +// [2, 5, 4, 4]. template void ParallelFor(ContextImpl* context, int start, int end, int num_threads, - const F& function, + F&& function, const CumulativeCostData* cumulative_cost_data, - const CumulativeCostFun& cumulative_cost_fun) { - using namespace parallel_for_details; + CumulativeCostFun&& cumulative_cost_fun) { CHECK_GT(num_threads, 0); if (start >= end) { return; } if (num_threads == 1 || end - start <= num_threads) { - ParallelFor(context, start, end, num_threads, function); + ParallelFor(context, start, end, num_threads, std::forward(function)); return; } // Creating several partitions allows us to tolerate imperfections of // partitioning and user-supplied iteration costs up to a certain extent - const int kNumPartitionsPerThread = 4; + constexpr int kNumPartitionsPerThread = 4; const int kMaxPartitions = num_threads * kNumPartitionsPerThread; - const std::vector partitions = ComputePartition( - start, end, kMaxPartitions, cumulative_cost_data, cumulative_cost_fun); + const auto& partitions = PartitionRangeForParallelFor( + start, + end, + kMaxPartitions, + cumulative_cost_data, + std::forward(cumulative_cost_fun)); CHECK_GT(partitions.size(), 1); - const int num_partitions = partitions.size() - 1; - ParallelFor(context, - 0, - num_partitions, - num_threads, - [&function, &partitions](int thread_id, int partition_id) { - const int partition_start = partitions[partition_id]; - const int partition_end = partitions[partition_id + 1]; - - InvokeOnSegment( - thread_id, partition_start, partition_end, function); - }); + ParallelFor( + context, start, end, num_threads, std::forward(function), partitions); } - -// Execute function for every element in the range [start, end) with at most -// num_threads, using the user provided partitioning. taking into account -// user-provided integer cumulative costs of iterations. -template -void ParallelFor(ContextImpl* context, - int start, - int end, - int num_threads, - const F& function, - const std::vector& partitions) { - using namespace parallel_for_details; - CHECK_GT(num_threads, 0); - if (start >= end) { - return; - } - CHECK_EQ(partitions.front(), start); - CHECK_EQ(partitions.back(), end); - if (num_threads == 1 || end - start <= num_threads) { - ParallelFor(context, start, end, num_threads, function); - return; - } - CHECK_GT(partitions.size(), 1); - const int num_partitions = partitions.size() - 1; - ParallelFor(context, - 0, - num_partitions, - num_threads, - [&function, &partitions](int thread_id, int partition_id) { - const int partition_start = partitions[partition_id]; - const int partition_end = partitions[partition_id + 1]; - - InvokeOnSegment( - thread_id, partition_start, partition_end, function); - }); -} - -// Evaluate vector expression in parallel -// Assuming LhsExpression and RhsExpression are some sort of -// column-vector expression, assignment lhs = rhs -// is eavluated over a set of contiguous blocks in parallel. -// This is expected to work well in the case of vector-based -// expressions (since they typically do not result into -// temporaries). -// This method expects lhs to be size-compatible with rhs -template -void ParallelAssign(ContextImpl* context, - int num_threads, - LhsExpression& lhs, - const RhsExpression& rhs) { - static_assert(LhsExpression::ColsAtCompileTime == 1); - static_assert(RhsExpression::ColsAtCompileTime == 1); - CHECK_EQ(lhs.rows(), rhs.rows()); - const int num_rows = lhs.rows(); - ParallelFor(context, - 0, - num_rows, - num_threads, - [&lhs, &rhs](const std::tuple& range) { - auto [start, end] = range; - lhs.segment(start, end - start) = - rhs.segment(start, end - start); - }); -} - -// Set vector to zero using num_threads -template -void ParallelSetZero(ContextImpl* context, - int num_threads, - VectorType& vector) { - ParallelSetZero(context, num_threads, vector.data(), vector.rows()); -} -void ParallelSetZero(ContextImpl* context, - int num_threads, - double* values, - int num_values); } // namespace ceres::internal -// Backend-specific implementations of ParallelInvoke -#include "ceres/internal/disable_warnings.h" -#include "ceres/parallel_for_cxx.h" - #endif // CERES_INTERNAL_PARALLEL_FOR_H_ diff --git a/internal/ceres/parallel_for_test.cc b/internal/ceres/parallel_for_test.cc index 871bb52d5..b3c210f06 100644 --- a/internal/ceres/parallel_for_test.cc +++ b/internal/ceres/parallel_for_test.cc @@ -41,6 +41,7 @@ #include "ceres/context_impl.h" #include "ceres/internal/config.h" +#include "ceres/parallel_vector_ops.h" #include "glog/logging.h" #include "gmock/gmock.h" #include "gtest/gtest.h" @@ -190,8 +191,6 @@ bool BruteForcePartition( // Basic test if MaxPartitionCostIsFeasible and BruteForcePartition agree on // simple test-cases TEST(GuidedParallelFor, MaxPartitionCostIsFeasible) { - using parallel_for_details::MaxPartitionCostIsFeasible; - std::vector costs, cumulative_costs, partition; costs = {1, 2, 3, 5, 0, 0, 0, 0, 0, 0, 7, 0, 0, 0, 0, 0, 0, 0}; cumulative_costs.resize(costs.size()); @@ -242,8 +241,6 @@ TEST(GuidedParallelFor, MaxPartitionCostIsFeasible) { // Randomized tests for MaxPartitionCostIsFeasible TEST(GuidedParallelFor, MaxPartitionCostIsFeasibleRandomized) { - using parallel_for_details::MaxPartitionCostIsFeasible; - std::vector costs, cumulative_costs, partition; const auto dummy_getter = [](const int v) { return v; }; @@ -315,9 +312,7 @@ TEST(GuidedParallelFor, MaxPartitionCostIsFeasibleRandomized) { } } -TEST(GuidedParallelFor, ComputePartition) { - using parallel_for_details::ComputePartition; - +TEST(GuidedParallelFor, PartitionRangeForParallelFor) { std::vector costs, cumulative_costs, partition; const auto dummy_getter = [](const int v) { return v; }; @@ -359,8 +354,8 @@ TEST(GuidedParallelFor, ComputePartition) { } } EXPECT_TRUE(first_admissible != 0 || total == 0); - partition = - ComputePartition(start, end, M, cumulative_costs.data(), dummy_getter); + partition = PartitionRangeForParallelFor( + start, end, M, cumulative_costs.data(), dummy_getter); ASSERT_GT(partition.size(), 1); EXPECT_EQ(partition.front(), start); EXPECT_EQ(partition.back(), end); diff --git a/internal/ceres/parallel_for_cxx.cc b/internal/ceres/parallel_invoke.cc similarity index 79% rename from internal/ceres/parallel_for_cxx.cc rename to internal/ceres/parallel_invoke.cc index 3c0671ed3..0e387c5ea 100644 --- a/internal/ceres/parallel_for_cxx.cc +++ b/internal/ceres/parallel_invoke.cc @@ -1,5 +1,5 @@ // Ceres Solver - A fast non-linear least squares minimizer -// Copyright 2018 Google Inc. All rights reserved. +// Copyright 2023 Google Inc. All rights reserved. // http://ceres-solver.org/ // // Redistribution and use in source and binary forms, with or without @@ -38,6 +38,7 @@ #include "ceres/internal/config.h" #include "ceres/parallel_for.h" +#include "ceres/parallel_vector_ops.h" #include "glog/logging.h" namespace ceres::internal { @@ -58,10 +59,12 @@ void BlockUntilFinished::Finished(int num_jobs_finished) { void BlockUntilFinished::Block() { std::unique_lock lock(mutex_); condition_.wait( - lock, [&]() { return num_total_jobs_finished_ == num_total_jobs_; }); + lock, [this]() { return num_total_jobs_finished_ == num_total_jobs_; }); } -ThreadPoolState::ThreadPoolState(int start, int end, int num_work_blocks) +ParallelInvokeState::ParallelInvokeState(int start, + int end, + int num_work_blocks) : start(start), end(end), num_work_blocks(num_work_blocks), @@ -71,20 +74,4 @@ ThreadPoolState::ThreadPoolState(int start, int end, int num_work_blocks) thread_id(0), block_until_finished(num_work_blocks) {} -int MaxNumThreadsAvailable() { return ThreadPool::MaxNumThreadsAvailable(); } - -void ParallelSetZero(ContextImpl* context, - int num_threads, - double* values, - int num_values) { - ParallelFor(context, - 0, - num_values, - num_threads, - [values](std::tuple range) { - auto [start, end] = range; - std::fill(values + start, values + end, 0.); - }); -} - } // namespace ceres::internal diff --git a/internal/ceres/parallel_for_cxx.h b/internal/ceres/parallel_invoke.h similarity index 85% rename from internal/ceres/parallel_for_cxx.h rename to internal/ceres/parallel_invoke.h index fb8715607..3e67e37bb 100644 --- a/internal/ceres/parallel_for_cxx.h +++ b/internal/ceres/parallel_invoke.h @@ -1,5 +1,5 @@ // Ceres Solver - A fast non-linear least squares minimizer -// Copyright 2018 Google Inc. All rights reserved. +// Copyright 2023 Google Inc. All rights reserved. // http://ceres-solver.org/ // // Redistribution and use in source and binary forms, with or without @@ -29,19 +29,47 @@ // Authors: vitus@google.com (Michael Vitus), // dmitriy.korchemkin@gmail.com (Dmitriy Korchemkin) -#ifndef CERES_INTERNAL_PARALLEL_FOR_CXX_H_ -#define CERES_INTERNAL_PARALLEL_FOR_CXX_H_ +#ifndef CERES_INTERNAL_PARALLEL_INVOKE_H_ +#define CERES_INTERNAL_PARALLEL_INVOKE_H_ #include -#include #include #include #include - -#include "ceres/internal/config.h" -#include "glog/logging.h" +#include +#include namespace ceres::internal { + +// InvokeWithThreadId handles passing thread_id to the function +template +void InvokeWithThreadId(int thread_id, F&& function, Args&&... args) { + constexpr bool kPassThreadId = std::is_invocable_v; + + if constexpr (kPassThreadId) { + function(thread_id, std::forward(args)...); + } else { + function(std::forward(args)...); + } +} + +// InvokeOnSegment either runs a loop over segment indices or passes it to the +// function +template +void InvokeOnSegment(int thread_id, std::tuple range, F&& function) { + constexpr bool kExplicitLoop = + std::is_invocable_v || std::is_invocable_v; + + if constexpr (kExplicitLoop) { + const auto [start, end] = range; + for (int i = start; i != end; ++i) { + InvokeWithThreadId(thread_id, std::forward(function), i); + } + } else { + InvokeWithThreadId(thread_id, std::forward(function), range); + } +} + // This class creates a thread safe barrier which will block until a // pre-specified number of threads call Finished. This allows us to block the // main thread until all the parallel threads are finished processing all the @@ -62,18 +90,17 @@ class BlockUntilFinished { std::mutex mutex_; std::condition_variable condition_; int num_total_jobs_finished_; - int num_total_jobs_; + const int num_total_jobs_; }; // Shared state between the parallel tasks. Each thread will use this // information to get the next block of work to be performed. -struct ThreadPoolState { +struct ParallelInvokeState { // The entire range [start, end) is split into num_work_blocks contiguous // disjoint intervals (blocks), which are as equal as possible given // total index count and requested number of blocks. // - // Those num_work_blocks blocks are then processed by num_workers - // workers + // Those num_work_blocks blocks are then processed in parallel. // // Total number of integer indices in interval [start, end) is // end - start, and when splitting them into num_work_blocks blocks @@ -87,7 +114,7 @@ struct ThreadPoolState { // Note that this splitting is optimal in the sense of maximal difference // between block sizes, since splitting into equal blocks is possible // if and only if number of indices is divisible by number of blocks. - ThreadPoolState(int start, int end, int num_work_blocks); + ParallelInvokeState(int start, int end, int num_work_blocks); // The start and end index of the for loop. const int start; @@ -135,12 +162,8 @@ struct ThreadPoolState { // A performance analysis has shown this implementation is on par with OpenMP // and TBB. template -void ParallelInvoke(ContextImpl* context, - int start, - int end, - int num_threads, - const F& function) { - using namespace parallel_for_details; +void ParallelInvoke( + ContextImpl* context, int start, int end, int num_threads, F&& function) { CHECK(context != nullptr); // Maximal number of work items scheduled for a single thread @@ -159,8 +182,8 @@ void ParallelInvoke(ContextImpl* context, // We use a std::shared_ptr because the main thread can finish all // the work before the tasks have been popped off the queue. So the // shared state needs to exist for the duration of all the tasks. - std::shared_ptr shared_state( - new ThreadPoolState(start, end, num_work_blocks)); + auto shared_state = + std::make_shared(start, end, num_work_blocks); // A function which tries to perform several chunks of work. auto task = [shared_state, num_threads, &function]() { @@ -213,7 +236,8 @@ void ParallelInvoke(ContextImpl* context, const int curr_end = curr_start + base_block_size + (block_id < num_base_p1_sized_blocks ? 1 : 0); // Perform each task in current block - InvokeOnSegment(thread_id, curr_start, curr_end, function); + const auto range = std::make_tuple(curr_start, curr_end); + InvokeOnSegment(thread_id, range, function); } shared_state->block_until_finished.Finished(num_jobs_finished); }; @@ -239,4 +263,4 @@ void ParallelInvoke(ContextImpl* context, } // namespace ceres::internal -#endif // CERES_INTERNAL_PARALLEL_FOR_CXX_H_ +#endif diff --git a/internal/ceres/parallel_vector_ops.h b/internal/ceres/parallel_vector_ops.h new file mode 100644 index 000000000..20d6dbcf7 --- /dev/null +++ b/internal/ceres/parallel_vector_ops.h @@ -0,0 +1,87 @@ +// Ceres Solver - A fast non-linear least squares minimizer +// Copyright 2023 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: vitus@google.com (Michael Vitus), +// dmitriy.korchemkin@gmail.com (Dmitriy Korchemkin) + +#ifndef CERES_INTERNAL_PARALLEL_VECTOR_OPS_H_ +#define CERES_INTERNAL_PARALLEL_VECTOR_OPS_H_ + +#include +#include + +#include "ceres/context_impl.h" +#include "ceres/internal/eigen.h" +#include "ceres/internal/export.h" +#include "ceres/parallel_for.h" + +namespace ceres::internal { + +// Evaluate vector expression in parallel +// Assuming LhsExpression and RhsExpression are some sort of +// column-vector expression, assignment lhs = rhs +// is eavluated over a set of contiguous blocks in parallel. +// This is expected to work well in the case of vector-based +// expressions (since they typically do not result into +// temporaries). +// This method expects lhs to be size-compatible with rhs +template +void ParallelAssign(ContextImpl* context, + int num_threads, + LhsExpression& lhs, + const RhsExpression& rhs) { + static_assert(LhsExpression::ColsAtCompileTime == 1); + static_assert(RhsExpression::ColsAtCompileTime == 1); + CHECK_EQ(lhs.rows(), rhs.rows()); + const int num_rows = lhs.rows(); + ParallelFor(context, + 0, + num_rows, + num_threads, + [&lhs, &rhs](const std::tuple& range) { + auto [start, end] = range; + lhs.segment(start, end - start) = + rhs.segment(start, end - start); + }); +} + +// Set vector to zero using num_threads +template +void ParallelSetZero(ContextImpl* context, + int num_threads, + VectorType& vector) { + ParallelSetZero(context, num_threads, vector.data(), vector.rows()); +} +void ParallelSetZero(ContextImpl* context, + int num_threads, + double* values, + int num_values); + +} // namespace ceres::internal + +#endif // CERES_INTERNAL_PARALLEL_FOR_H_ diff --git a/internal/ceres/partition_range_for_parallel_for.h b/internal/ceres/partition_range_for_parallel_for.h new file mode 100644 index 000000000..309d7a89d --- /dev/null +++ b/internal/ceres/partition_range_for_parallel_for.h @@ -0,0 +1,150 @@ +// Ceres Solver - A fast non-linear least squares minimizer +// Copyright 2023 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: vitus@google.com (Michael Vitus), +// dmitriy.korchemkin@gmail.com (Dmitriy Korchemkin) + +#ifndef CERES_INTERNAL_PARTITION_RANGE_FOR_PARALLEL_FOR_H_ +#define CERES_INTERNAL_PARTITION_RANGE_FOR_PARALLEL_FOR_H_ + +#include +#include + +namespace ceres::internal { +// Check if it is possible to split range [start; end) into at most +// max_num_partitions contiguous partitions of cost not greater than +// max_partition_cost. Inclusive integer cumulative costs are provided by +// cumulative_cost_data objects, with cumulative_cost_offset being a total cost +// of all indices (starting from zero) preceding start element. Cumulative costs +// are returned by cumulative_cost_fun called with a reference to +// cumulative_cost_data element with index from range[start; end), and should be +// non-decreasing. Partition of the range is returned via partition argument +template +bool MaxPartitionCostIsFeasible(int start, + int end, + int max_num_partitions, + int max_partition_cost, + int cumulative_cost_offset, + const CumulativeCostData* cumulative_cost_data, + CumulativeCostFun&& cumulative_cost_fun, + std::vector* partition) { + partition->clear(); + partition->push_back(start); + int partition_start = start; + int cost_offset = cumulative_cost_offset; + + while (partition_start < end) { + // Already have max_num_partitions + if (partition->size() > max_num_partitions) { + return false; + } + const int target = max_partition_cost + cost_offset; + const int partition_end = + std::partition_point( + cumulative_cost_data + partition_start, + cumulative_cost_data + end, + [&cumulative_cost_fun, target](const CumulativeCostData& item) { + return cumulative_cost_fun(item) <= target; + }) - + cumulative_cost_data; + // Unable to make a partition from a single element + if (partition_end == partition_start) { + return false; + } + + const int cost_last = + cumulative_cost_fun(cumulative_cost_data[partition_end - 1]); + partition->push_back(partition_end); + partition_start = partition_end; + cost_offset = cost_last; + } + return true; +} + +// Split integer interval [start, end) into at most max_num_partitions +// contiguous intervals, minimizing maximal total cost of a single interval. +// Inclusive integer cumulative costs for each (zero-based) index are provided +// by cumulative_cost_data objects, and are returned by cumulative_cost_fun call +// with a reference to one of the objects from range [start, end) +template +std::vector PartitionRangeForParallelFor( + int start, + int end, + int max_num_partitions, + const CumulativeCostData* cumulative_cost_data, + CumulativeCostFun&& cumulative_cost_fun) { + // Given maximal partition cost, it is possible to verify if it is admissible + // and obtain corresponding partition using MaxPartitionCostIsFeasible + // function. In order to find the lowest admissible value, a binary search + // over all potentially optimal cost values is being performed + const int cumulative_cost_last = + cumulative_cost_fun(cumulative_cost_data[end - 1]); + const int cumulative_cost_offset = + start ? cumulative_cost_fun(cumulative_cost_data[start - 1]) : 0; + const int total_cost = cumulative_cost_last - cumulative_cost_offset; + + // Minimal maximal partition cost is not smaller than the average + // We will use non-inclusive lower bound + int partition_cost_lower_bound = total_cost / max_num_partitions - 1; + // Minimal maximal partition cost is not larger than the total cost + // Upper bound is inclusive + int partition_cost_upper_bound = total_cost; + + std::vector partition; + // Range partition corresponding to the latest evaluated upper bound. + // A single segment covering the whole input interval [start, end) corresponds + // to minimal maximal partition cost of total_cost. + std::vector partition_upper_bound = {start, end}; + // Binary search over partition cost, returning the lowest admissible cost + while (partition_cost_upper_bound - partition_cost_lower_bound > 1) { + partition.reserve(max_num_partitions + 1); + const int partition_cost = + partition_cost_lower_bound + + (partition_cost_upper_bound - partition_cost_lower_bound) / 2; + bool admissible = MaxPartitionCostIsFeasible( + start, + end, + max_num_partitions, + partition_cost, + cumulative_cost_offset, + cumulative_cost_data, + std::forward(cumulative_cost_fun), + &partition); + if (admissible) { + partition_cost_upper_bound = partition_cost; + std::swap(partition, partition_upper_bound); + } else { + partition_cost_lower_bound = partition_cost; + } + } + + return partition_upper_bound; +} +} // namespace ceres::internal + +#endif diff --git a/internal/ceres/partitioned_matrix_view_impl.h b/internal/ceres/partitioned_matrix_view_impl.h index e685f5b76..385ceb1f9 100644 --- a/internal/ceres/partitioned_matrix_view_impl.h +++ b/internal/ceres/partitioned_matrix_view_impl.h @@ -37,6 +37,7 @@ #include "ceres/block_structure.h" #include "ceres/internal/eigen.h" #include "ceres/parallel_for.h" +#include "ceres/partition_range_for_parallel_for.h" #include "ceres/partitioned_matrix_view.h" #include "ceres/small_blas.h" #include "glog/logging.h" @@ -87,14 +88,14 @@ PartitionedMatrixView:: const int num_threads = options_.num_threads; if (transpose_bs != nullptr && num_threads > 1) { int kMaxPartitions = num_threads * 4; - e_cols_partition_ = parallel_for_details::ComputePartition( + e_cols_partition_ = PartitionRangeForParallelFor( 0, num_col_blocks_e_, kMaxPartitions, transpose_bs->rows.data(), [](const CompressedRow& row) { return row.cumulative_nnz; }); - f_cols_partition_ = parallel_for_details::ComputePartition( + f_cols_partition_ = PartitionRangeForParallelFor( num_col_blocks_e_, num_col_blocks_e_ + num_col_blocks_f_, kMaxPartitions, diff --git a/internal/ceres/program_evaluator.h b/internal/ceres/program_evaluator.h index 7ccbbd128..372b2602b 100644 --- a/internal/ceres/program_evaluator.h +++ b/internal/ceres/program_evaluator.h @@ -43,7 +43,7 @@ // residual jacobians are written directly into their final position in the // block sparse matrix by the user's CostFunction; there is no copying. // -// The evaluation is threaded with OpenMP or C++ threads. +// The evaluation is threaded with C++ threads. // // The EvaluatePreparer and JacobianWriter interfaces are as follows: // @@ -96,6 +96,7 @@ #include "ceres/execution_summary.h" #include "ceres/internal/eigen.h" #include "ceres/parallel_for.h" +#include "ceres/parallel_vector_ops.h" #include "ceres/parameter_block.h" #include "ceres/program.h" #include "ceres/residual_block.h" diff --git a/internal/ceres/wall_time.h b/internal/ceres/wall_time.h index 36936ee7e..d50349337 100644 --- a/internal/ceres/wall_time.h +++ b/internal/ceres/wall_time.h @@ -41,10 +41,8 @@ namespace ceres::internal { -// Returns time, in seconds, from some arbitrary starting point. If -// OpenMP is available then the high precision openmp_get_wtime() -// function is used. Otherwise on unixes, gettimeofday is used. The -// granularity is in seconds on windows systems. +// Returns time, in seconds, from some arbitrary starting point. On unixes, +// gettimeofday is used. The granularity is microseconds. CERES_NO_EXPORT double WallTimeInSeconds(); // Log a series of events, recording for each event the time elapsed