From 17cffc5006df99c277f02cc8301108e7a87cbb25 Mon Sep 17 00:00:00 2001 From: Sandro Wenzel Date: Thu, 5 Nov 2026 21:32:02 +0100 Subject: [PATCH] Make the magnetic field evaluation reentrant This removes the shared scratch state from the Chebyshev field evaluation, so that one field can be evaluated from several threads. - Chebyshev3D keeps the mapped coordinates in a local array instead of a member. - Chebyshev3DCalc evaluates into scratch on the stack, or into a thread-local buffer for parameterisations larger than 64 rows or columns. The ALICE maps use at most 14. - MagneticWrapperChebyshev no longer uses static coordinates in the TPC integrals. - The arithmetic is unchanged, so the field values are bit-identical. - A unit test evaluates the field from 8 threads and compares with a sequential pass. This prepares the simulation for multithreaded Geant4, where all workers share one field. https://its.cern.ch/jira/browse/O2-7305 Co-Authored-By: Claude Opus 5.5 --- Common/Field/src/MagneticWrapperChebyshev.cxx | 4 +- Common/Field/test/testMagneticField.cxx | 41 ++++++++++++++ .../MathUtils/include/MathUtils/Chebyshev3D.h | 50 ++++++++++------- .../include/MathUtils/Chebyshev3DCalc.h | 55 ++++++++++++++++--- Common/MathUtils/src/Chebyshev3DCalc.cxx | 42 ++++++++------ 5 files changed, 144 insertions(+), 48 deletions(-) diff --git a/Common/Field/src/MagneticWrapperChebyshev.cxx b/Common/Field/src/MagneticWrapperChebyshev.cxx index 98f1e191129b4..4d3a3a43c1990 100644 --- a/Common/Field/src/MagneticWrapperChebyshev.cxx +++ b/Common/Field/src/MagneticWrapperChebyshev.cxx @@ -690,7 +690,7 @@ Int_t MagneticWrapperChebyshev::findTPCRatSegment(const Double_t* rpz) const void MagneticWrapperChebyshev::getTPCIntegral(const Double_t* xyz, Double_t* b) const { - static Double_t rphiz[3]; + Double_t rphiz[3]; // TPCInt region // convert coordinates to cyl system @@ -712,7 +712,7 @@ void MagneticWrapperChebyshev::getTPCIntegral(const Double_t* xyz, Double_t* b) void MagneticWrapperChebyshev::getTPCRatIntegral(const Double_t* xyz, Double_t* b) const { - static Double_t rphiz[3]; + Double_t rphiz[3]; // TPCRatIntegral region // convert coordinates to cylindrical system diff --git a/Common/Field/test/testMagneticField.cxx b/Common/Field/test/testMagneticField.cxx index df300a4ac7194..2c41ce0ec7f67 100644 --- a/Common/Field/test/testMagneticField.cxx +++ b/Common/Field/test/testMagneticField.cxx @@ -21,6 +21,10 @@ #include #include #include +#include +#include +#include +#include using namespace o2::field; @@ -144,3 +148,40 @@ BOOST_AUTO_TEST_CASE(MagneticField_reinitialization_test) } } } + +BOOST_AUTO_TEST_CASE(MagneticField_threads) +{ + // evaluating one field from several threads at once gives bit-identical values to a sequential pass + MagneticField fld("Maps", "Maps", 1., 1., o2::field::MagFieldParam::k5kG); + constexpr int npoints = 20000; + constexpr int nthreads = 8; + std::vector> xyz(npoints); + for (auto& p : xyz) { + p = {gRandom->Uniform(-450., 450.), gRandom->Uniform(-450., 450.), gRandom->Uniform(-1500., 1500.)}; + } + // field, Bz, TPC integral and TPC ratio integral for every point + constexpr int nval = 10; + auto evaluate = [&fld, &xyz](std::vector& out) { + out.resize(npoints * nval); + for (int i = 0; i < npoints; ++i) { + fld.Field(xyz[i].data(), &out[i * nval]); + out[i * nval + 3] = fld.getBz(xyz[i].data()); + fld.getTPCIntegral(xyz[i].data(), &out[i * nval + 4]); + fld.getTPCRatIntegral(xyz[i].data(), &out[i * nval + 7]); + } + }; + std::vector reference; + evaluate(reference); + std::vector> results(nthreads); + std::vector threads; + for (int t = 0; t < nthreads; ++t) { + threads.emplace_back(evaluate, std::ref(results[t])); + } + for (auto& th : threads) { + th.join(); + } + for (int t = 0; t < nthreads; ++t) { + auto diff = std::mismatch(reference.begin(), reference.end(), results[t].begin()); + BOOST_CHECK_MESSAGE(diff.first == reference.end(), "thread " << t << " differs at value " << (diff.first - reference.begin())); + } +} diff --git a/Common/MathUtils/include/MathUtils/Chebyshev3D.h b/Common/MathUtils/include/MathUtils/Chebyshev3D.h index c7addd170b61b..0aa18de6d2ec3 100644 --- a/Common/MathUtils/include/MathUtils/Chebyshev3D.h +++ b/Common/MathUtils/include/MathUtils/Chebyshev3D.h @@ -292,52 +292,57 @@ inline Bool_t Chebyshev3D::isInside(const Double_t* par) const /// Evaluates Chebyshev parameterization for 3d->DimOut function inline void Chebyshev3D::Eval(const Float_t* par, Float_t* res) { + Float_t coef[3]; // mapped arguments, local so that evaluation is reentrant for (int i = 3; i--;) { - mTemporaryCoefficient[i] = mapToInternal(par[i], i); + coef[i] = mapToInternal(par[i], i); } for (int i = mOutputArrayDimension; i--;) { - res[i] = getChebyshevCalc(i)->Eval(mTemporaryCoefficient); + res[i] = getChebyshevCalc(i)->Eval(coef); } } /// Evaluates Chebyshev parameterization for 3d->DimOut function inline void Chebyshev3D::Eval(const Double_t* par, Double_t* res) { + Float_t coef[3]; // mapped arguments, local so that evaluation is reentrant for (int i = 3; i--;) { - mTemporaryCoefficient[i] = mapToInternal(par[i], i); + coef[i] = mapToInternal(par[i], i); } for (int i = mOutputArrayDimension; i--;) { - res[i] = getChebyshevCalc(i)->Eval(mTemporaryCoefficient); + res[i] = getChebyshevCalc(i)->Eval(coef); } } /// Evaluates Chebyshev parameterization for idim-th output dimension of 3d->DimOut function inline Double_t Chebyshev3D::Eval(const Double_t* par, int idim) { + Float_t coef[3]; // mapped arguments, local so that evaluation is reentrant for (int i = 3; i--;) { - mTemporaryCoefficient[i] = mapToInternal(par[i], i); + coef[i] = mapToInternal(par[i], i); } - return getChebyshevCalc(idim)->Eval(mTemporaryCoefficient); + return getChebyshevCalc(idim)->Eval(coef); } /// Evaluates Chebyshev parameterization for idim-th output dimension of 3d->DimOut function inline Float_t Chebyshev3D::Eval(const Float_t* par, int idim) { + Float_t coef[3]; // mapped arguments, local so that evaluation is reentrant for (int i = 3; i--;) { - mTemporaryCoefficient[i] = mapToInternal(par[i], i); + coef[i] = mapToInternal(par[i], i); } - return getChebyshevCalc(idim)->Eval(mTemporaryCoefficient); + return getChebyshevCalc(idim)->Eval(coef); } /// Returns the gradient matrix inline void Chebyshev3D::evaluateDerivative3D(const Float_t* par, Float_t dbdr[3][3]) { + Float_t coef[3]; // mapped arguments, local so that evaluation is reentrant for (int i = 3; i--;) { - mTemporaryCoefficient[i] = mapToInternal(par[i], i); + coef[i] = mapToInternal(par[i], i); } for (int ib = 3; ib--;) { for (int id = 3; id--;) { - dbdr[ib][id] = getChebyshevCalc(ib)->evaluateDerivative(id, mTemporaryCoefficient) * mBoundaryMappingScale[id]; + dbdr[ib][id] = getChebyshevCalc(ib)->evaluateDerivative(id, coef) * mBoundaryMappingScale[id]; } } } @@ -345,13 +350,14 @@ inline void Chebyshev3D::evaluateDerivative3D(const Float_t* par, Float_t dbdr[3 /// Returns the gradient matrix inline void Chebyshev3D::evaluateDerivative3D2(const Float_t* par, Float_t dbdrdr[3][3][3]) { + Float_t coef[3]; // mapped arguments, local so that evaluation is reentrant for (int i = 3; i--;) { - mTemporaryCoefficient[i] = mapToInternal(par[i], i); + coef[i] = mapToInternal(par[i], i); } for (int ib = 3; ib--;) { for (int id = 3; id--;) { for (int id1 = 3; id1--;) { - dbdrdr[ib][id][id1] = getChebyshevCalc(ib)->evaluateDerivative2(id, id1, mTemporaryCoefficient) * + dbdrdr[ib][id][id1] = getChebyshevCalc(ib)->evaluateDerivative2(id, id1, coef) * mBoundaryMappingScale[id] * mBoundaryMappingScale[id1]; } } @@ -361,22 +367,24 @@ inline void Chebyshev3D::evaluateDerivative3D2(const Float_t* par, Float_t dbdrd // Evaluates Chebyshev parameterization derivative for 3d->DimOut function inline void Chebyshev3D::evaluateDerivative(int dimd, const Float_t* par, Float_t* res) { + Float_t coef[3]; // mapped arguments, local so that evaluation is reentrant for (int i = 3; i--;) { - mTemporaryCoefficient[i] = mapToInternal(par[i], i); + coef[i] = mapToInternal(par[i], i); } for (int i = mOutputArrayDimension; i--;) { - res[i] = getChebyshevCalc(i)->evaluateDerivative(dimd, mTemporaryCoefficient) * mBoundaryMappingScale[dimd]; + res[i] = getChebyshevCalc(i)->evaluateDerivative(dimd, coef) * mBoundaryMappingScale[dimd]; }; } // Evaluates Chebyshev parameterization 2nd derivative over dimd1 and dimd2 dimensions for 3d->DimOut function inline void Chebyshev3D::evaluateDerivative2(int dimd1, int dimd2, const Float_t* par, Float_t* res) { + Float_t coef[3]; // mapped arguments, local so that evaluation is reentrant for (int i = 3; i--;) { - mTemporaryCoefficient[i] = mapToInternal(par[i], i); + coef[i] = mapToInternal(par[i], i); } for (int i = mOutputArrayDimension; i--;) { - res[i] = getChebyshevCalc(i)->evaluateDerivative2(dimd1, dimd2, mTemporaryCoefficient) * + res[i] = getChebyshevCalc(i)->evaluateDerivative2(dimd1, dimd2, coef) * mBoundaryMappingScale[dimd1] * mBoundaryMappingScale[dimd2]; } } @@ -385,20 +393,22 @@ inline void Chebyshev3D::evaluateDerivative2(int dimd1, int dimd2, const Float_t /// function inline Float_t Chebyshev3D::evaluateDerivative(int dimd, const Float_t* par, int idim) { + Float_t coef[3]; // mapped arguments, local so that evaluation is reentrant for (int i = 3; i--;) { - mTemporaryCoefficient[i] = mapToInternal(par[i], i); + coef[i] = mapToInternal(par[i], i); } - return getChebyshevCalc(idim)->evaluateDerivative(dimd, mTemporaryCoefficient) * mBoundaryMappingScale[dimd]; + return getChebyshevCalc(idim)->evaluateDerivative(dimd, coef) * mBoundaryMappingScale[dimd]; } /// Evaluates Chebyshev parameterization 2ns derivative over dimd1 and dimd2 dimensions for idim-th output dimension of /// 3d->DimOut function inline Float_t Chebyshev3D::evaluateDerivative2(int dimd1, int dimd2, const Float_t* par, int idim) { + Float_t coef[3]; // mapped arguments, local so that evaluation is reentrant for (int i = 3; i--;) { - mTemporaryCoefficient[i] = mapToInternal(par[i], i); + coef[i] = mapToInternal(par[i], i); } - return getChebyshevCalc(idim)->evaluateDerivative2(dimd1, dimd2, mTemporaryCoefficient) * + return getChebyshevCalc(idim)->evaluateDerivative2(dimd1, dimd2, coef) * mBoundaryMappingScale[dimd1] * mBoundaryMappingScale[dimd2]; } diff --git a/Common/MathUtils/include/MathUtils/Chebyshev3DCalc.h b/Common/MathUtils/include/MathUtils/Chebyshev3DCalc.h index 5f611172928eb..9ccd2daef42b1 100644 --- a/Common/MathUtils/include/MathUtils/Chebyshev3DCalc.h +++ b/Common/MathUtils/include/MathUtils/Chebyshev3DCalc.h @@ -19,6 +19,8 @@ #include // for TNamed #include // for FILE, stdout #include // for std::fma +#include +#include #include "Rtypes.h" // for Float_t, UShort_t, Int_t, Double_t, etc class TString; @@ -171,6 +173,37 @@ class Chebyshev3DCalc : public TNamed Double_t Eval(const Double_t* par) const; private: + /// Scratch arrays for one evaluation: on the stack for the usual sizes, thread-local beyond them, + /// so that evaluation is reentrant. The ALICE field maps use at most 14 rows and columns. + class EvalScratch + { + public: + EvalScratch(int nrows, int ncols) + { + if (nrows <= kStackSize && ncols <= kStackSize) { + mRows = mStackRows; + mCols = mStackCols; + } else { + thread_local std::vector rows, cols; + rows.resize(std::max(nrows, rows.size())); + cols.resize(std::max(ncols, cols.size())); + mRows = rows.data(); + mCols = cols.data(); + } + } + EvalScratch(const EvalScratch&) = delete; + EvalScratch& operator=(const EvalScratch&) = delete; + Float_t* rows() { return mRows; } + Float_t* cols() { return mCols; } + + private: + static constexpr int kStackSize = 64; + Float_t mStackRows[kStackSize]; + Float_t mStackCols[kStackSize]; + Float_t* mRows; + Float_t* mCols; + }; + Int_t mNumberOfCoefficients; ///< total number of coeeficients Int_t mNumberOfRows; ///< number of significant rows in the 3D coeffs matrix Int_t mNumberOfColumns; ///< max number of significant cols in the 3D coeffs matrix @@ -188,8 +221,8 @@ class Chebyshev3DCalc : public TNamed // coeffs for col/row Float_t* mCoefficients; //[mNumberOfCoefficients] array of Chebyshev coefficients - Float_t* mTemporaryCoefficients2D; //[mNumberOfColumns] temp. coeffs for 2d summation - Float_t* mTemporaryCoefficients1D; //[mNumberOfRows] temp. coeffs for 1d summation + Float_t* mTemporaryCoefficients2D; //[mNumberOfColumns] no longer used for evaluation + Float_t* mTemporaryCoefficients1D; //[mNumberOfRows] no longer used for evaluation ClassDefOverride(o2::math_utils::Chebyshev3DCalc, 2) // Class for interpolation of 3D->1 function by Chebyshev parametrization @@ -223,32 +256,38 @@ inline Float_t Chebyshev3DCalc::chebyshevEvaluation1D(Float_t x, const Float_t* /// VERY IMPORTANT: par must contain the function arguments ALREADY MAPPED to [-1:1] interval inline Float_t Chebyshev3DCalc::Eval(const Float_t* par) const { + EvalScratch scratch(mNumberOfRows, mNumberOfColumns); + Float_t* coefs2D = scratch.cols(); + Float_t* coefs1D = scratch.rows(); for (int id0 = mNumberOfRows; id0--;) { int nCLoc = mNumberOfColumnsAtRow[id0]; // number of significant coefs on this row int col0 = mColumnAtRowBeginning[id0]; // beginning of local column in the 2D boundary matrix for (int id1 = nCLoc; id1--;) { int id = id1 + col0; - mTemporaryCoefficients2D[id1] = chebyshevEvaluation1D(par[2], mCoefficients + mCoefficientBound2D1[id], mCoefficientBound2D0[id]); + coefs2D[id1] = chebyshevEvaluation1D(par[2], mCoefficients + mCoefficientBound2D1[id], mCoefficientBound2D0[id]); } - mTemporaryCoefficients1D[id0] = chebyshevEvaluation1D(par[1], mTemporaryCoefficients2D, nCLoc); + coefs1D[id0] = chebyshevEvaluation1D(par[1], coefs2D, nCLoc); } - return chebyshevEvaluation1D(par[0], mTemporaryCoefficients1D, mNumberOfRows); + return chebyshevEvaluation1D(par[0], coefs1D, mNumberOfRows); } /// Evaluates Chebyshev parameterization for 3D function. /// VERY IMPORTANT: par must contain the function arguments ALREADY MAPPED to [-1:1] interval inline Double_t Chebyshev3DCalc::Eval(const Double_t* par) const { + EvalScratch scratch(mNumberOfRows, mNumberOfColumns); + Float_t* coefs2D = scratch.cols(); + Float_t* coefs1D = scratch.rows(); for (int id0 = mNumberOfRows; id0--;) { int nCLoc = mNumberOfColumnsAtRow[id0]; // number of significant coefs on this row int col0 = mColumnAtRowBeginning[id0]; // beginning of local column in the 2D boundary matrix for (int id1 = nCLoc; id1--;) { int id = id1 + col0; - mTemporaryCoefficients2D[id1] = chebyshevEvaluation1D(par[2], mCoefficients + mCoefficientBound2D1[id], mCoefficientBound2D0[id]); + coefs2D[id1] = chebyshevEvaluation1D(par[2], mCoefficients + mCoefficientBound2D1[id], mCoefficientBound2D0[id]); } - mTemporaryCoefficients1D[id0] = chebyshevEvaluation1D(par[1], mTemporaryCoefficients2D, nCLoc); + coefs1D[id0] = chebyshevEvaluation1D(par[1], coefs2D, nCLoc); } - return chebyshevEvaluation1D(par[0], mTemporaryCoefficients1D, mNumberOfRows); + return chebyshevEvaluation1D(par[0], coefs1D, mNumberOfRows); } } // namespace math_utils } // namespace o2 diff --git a/Common/MathUtils/src/Chebyshev3DCalc.cxx b/Common/MathUtils/src/Chebyshev3DCalc.cxx index 7a402215a1cf3..6772ab887ff4e 100644 --- a/Common/MathUtils/src/Chebyshev3DCalc.cxx +++ b/Common/MathUtils/src/Chebyshev3DCalc.cxx @@ -205,11 +205,14 @@ void Chebyshev3DCalc::Print(const Option_t*) const Float_t Chebyshev3DCalc::evaluateDerivative(int dim, const Float_t* par) const { + EvalScratch scratch(mNumberOfRows, mNumberOfColumns); + Float_t* coefs2D = scratch.cols(); + Float_t* coefs1D = scratch.rows(); int ncfRC; for (int id0 = mNumberOfRows; id0--;) { int nCLoc = mNumberOfColumnsAtRow[id0]; // number of significant coefs on this row if (!nCLoc) { - mTemporaryCoefficients1D[id0] = 0; + coefs1D[id0] = 0; continue; } // @@ -217,62 +220,65 @@ Float_t Chebyshev3DCalc::evaluateDerivative(int dim, const Float_t* par) const for (int id1 = nCLoc; id1--;) { int id = id1 + col0; if (!(ncfRC = mCoefficientBound2D0[id])) { - mTemporaryCoefficients2D[id1] = 0; + coefs2D[id1] = 0; continue; } if (dim == 2) { - mTemporaryCoefficients2D[id1] = + coefs2D[id1] = chebyshevEvaluation1Derivative(par[2], mCoefficients + mCoefficientBound2D1[id], ncfRC); } else { - mTemporaryCoefficients2D[id1] = chebyshevEvaluation1D(par[2], mCoefficients + mCoefficientBound2D1[id], ncfRC); + coefs2D[id1] = chebyshevEvaluation1D(par[2], mCoefficients + mCoefficientBound2D1[id], ncfRC); } } if (dim == 1) { - mTemporaryCoefficients1D[id0] = chebyshevEvaluation1Derivative(par[1], mTemporaryCoefficients2D, nCLoc); + coefs1D[id0] = chebyshevEvaluation1Derivative(par[1], coefs2D, nCLoc); } else { - mTemporaryCoefficients1D[id0] = chebyshevEvaluation1D(par[1], mTemporaryCoefficients2D, nCLoc); + coefs1D[id0] = chebyshevEvaluation1D(par[1], coefs2D, nCLoc); } } - return (dim == 0) ? chebyshevEvaluation1Derivative(par[0], mTemporaryCoefficients1D, mNumberOfRows) - : chebyshevEvaluation1D(par[0], mTemporaryCoefficients1D, mNumberOfRows); + return (dim == 0) ? chebyshevEvaluation1Derivative(par[0], coefs1D, mNumberOfRows) + : chebyshevEvaluation1D(par[0], coefs1D, mNumberOfRows); } Float_t Chebyshev3DCalc::evaluateDerivative2(int dim1, int dim2, const Float_t* par) const { + EvalScratch scratch(mNumberOfRows, mNumberOfColumns); + Float_t* coefs2D = scratch.cols(); + Float_t* coefs1D = scratch.rows(); Bool_t same = dim1 == dim2; int ncfRC; for (int id0 = mNumberOfRows; id0--;) { int nCLoc = mNumberOfColumnsAtRow[id0]; // number of significant coefs on this row if (!nCLoc) { - mTemporaryCoefficients1D[id0] = 0; + coefs1D[id0] = 0; continue; } int col0 = mColumnAtRowBeginning[id0]; // beginning of local column in the 2D boundary matrix for (int id1 = nCLoc; id1--;) { int id = id1 + col0; if (!(ncfRC = mCoefficientBound2D0[id])) { - mTemporaryCoefficients2D[id1] = 0; + coefs2D[id1] = 0; continue; } if (dim1 == 2 || dim2 == 2) { - mTemporaryCoefficients2D[id1] = + coefs2D[id1] = same ? chebyshevEvaluation1Derivative2(par[2], mCoefficients + mCoefficientBound2D1[id], ncfRC) : chebyshevEvaluation1Derivative(par[2], mCoefficients + mCoefficientBound2D1[id], ncfRC); } else { - mTemporaryCoefficients2D[id1] = chebyshevEvaluation1D(par[2], mCoefficients + mCoefficientBound2D1[id], ncfRC); + coefs2D[id1] = chebyshevEvaluation1D(par[2], mCoefficients + mCoefficientBound2D1[id], ncfRC); } } if (dim1 == 1 || dim2 == 1) { - mTemporaryCoefficients1D[id0] = same ? chebyshevEvaluation1Derivative2(par[1], mTemporaryCoefficients2D, nCLoc) - : chebyshevEvaluation1Derivative(par[1], mTemporaryCoefficients2D, nCLoc); + coefs1D[id0] = same ? chebyshevEvaluation1Derivative2(par[1], coefs2D, nCLoc) + : chebyshevEvaluation1Derivative(par[1], coefs2D, nCLoc); } else { - mTemporaryCoefficients1D[id0] = chebyshevEvaluation1D(par[1], mTemporaryCoefficients2D, nCLoc); + coefs1D[id0] = chebyshevEvaluation1D(par[1], coefs2D, nCLoc); } } return (dim1 == 0 || dim2 == 0) - ? (same ? chebyshevEvaluation1Derivative2(par[0], mTemporaryCoefficients1D, mNumberOfRows) - : chebyshevEvaluation1Derivative(par[0], mTemporaryCoefficients1D, mNumberOfRows)) - : chebyshevEvaluation1D(par[0], mTemporaryCoefficients1D, mNumberOfRows); + ? (same ? chebyshevEvaluation1Derivative2(par[0], coefs1D, mNumberOfRows) + : chebyshevEvaluation1Derivative(par[0], coefs1D, mNumberOfRows)) + : chebyshevEvaluation1D(par[0], coefs1D, mNumberOfRows); } #ifdef _INC_CREATION_Chebyshev3D_