From f11bada966b49e5a31b26ce3f388ddc9de44a6f2 Mon Sep 17 00:00:00 2001 From: YL <20972148+LiYunyang@users.noreply.github.com> Date: Sat, 31 Jan 2026 15:41:16 -0600 Subject: [PATCH 1/8] [feat] multithread `RemoveWeights` --- maps/src/maputils.cxx | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/maps/src/maputils.cxx b/maps/src/maputils.cxx index ee46ec1d..15b81b05 100644 --- a/maps/src/maputils.cxx +++ b/maps/src/maputils.cxx @@ -5,6 +5,9 @@ #include #include #include +#ifdef _OPENMP +#include +#endif #include #include @@ -37,11 +40,17 @@ void RemoveWeights(G3SkyMap &T, G3SkyMap &Q, G3SkyMap &U, const G3SkyMapWeights T.ConvertToDense(); Q.ConvertToDense(); U.ConvertToDense(); + #ifdef _OPENMP + #pragma omp parallel for schedule(static) + #endif for (size_t pix = 0; pix < T.size(); pix++) { StokesVector v(T[pix], Q[pix], U[pix]); v /= W.at(pix); } } else { + #ifdef _OPENMP + #pragma omp parallel for schedule(static) + #endif for (size_t pix = 0; pix < W.size(); pix++) { double t = T.at(pix); MuellerMatrix m = W.at(pix); From 02db4879615aa56630d78c96d7982c23aeb5aebb Mon Sep 17 00:00:00 2001 From: YL <20972148+LiYunyang@users.noreply.github.com> Date: Sun, 1 Feb 2026 18:12:30 -0600 Subject: [PATCH 2/8] multithread `ApplyWeights` --- maps/src/maputils.cxx | 3 +++ 1 file changed, 3 insertions(+) diff --git a/maps/src/maputils.cxx b/maps/src/maputils.cxx index 15b81b05..435ab498 100644 --- a/maps/src/maputils.cxx +++ b/maps/src/maputils.cxx @@ -130,6 +130,9 @@ void ApplyWeights(G3SkyMap &T, G3SkyMap &Q, G3SkyMap &U, const G3SkyMapWeights & g3_assert(!Q.weighted); g3_assert(!U.weighted); + #ifdef _OPENMP + #pragma omp parallel for schedule(static) + #endif for (size_t pix = 0; pix < T.size(); pix++) { if (T.at(pix) == 0 && Q.at(pix) == 0 && U.at(pix) == 0) continue; From f838b2458cb46494e8478f4b831828b8ce5bf3a5 Mon Sep 17 00:00:00 2001 From: YL <20972148+LiYunyang@users.noreply.github.com> Date: Sun, 1 Feb 2026 18:12:30 -0600 Subject: [PATCH 3/8] multithread `ApplyWeights` --- maps/src/maputils.cxx | 3 +++ 1 file changed, 3 insertions(+) diff --git a/maps/src/maputils.cxx b/maps/src/maputils.cxx index 15b81b05..435ab498 100644 --- a/maps/src/maputils.cxx +++ b/maps/src/maputils.cxx @@ -130,6 +130,9 @@ void ApplyWeights(G3SkyMap &T, G3SkyMap &Q, G3SkyMap &U, const G3SkyMapWeights & g3_assert(!Q.weighted); g3_assert(!U.weighted); + #ifdef _OPENMP + #pragma omp parallel for schedule(static) + #endif for (size_t pix = 0; pix < T.size(); pix++) { if (T.at(pix) == 0 && Q.at(pix) == 0 && U.at(pix) == 0) continue; From 52b81648b0d0eced624e994b42c51e739a4a421e Mon Sep 17 00:00:00 2001 From: YL <20972148+LiYunyang@users.noreply.github.com> Date: Wed, 4 Feb 2026 12:32:27 -0600 Subject: [PATCH 4/8] multithread `RemoveWeightsT` --- maps/src/maputils.cxx | 3 +++ 1 file changed, 3 insertions(+) diff --git a/maps/src/maputils.cxx b/maps/src/maputils.cxx index 435ab498..24501e92 100644 --- a/maps/src/maputils.cxx +++ b/maps/src/maputils.cxx @@ -94,6 +94,9 @@ void RemoveWeightsT(G3SkyMap &T, const G3SkyMapWeights &W, bool zero_nans) T.ConvertToDense(); T /= *(W.TT); } else { + #ifdef _OPENMP + #pragma omp parallel for schedule(static) + #endif for (size_t pix = 0; pix < W.size(); pix++) { double t = T.at(pix); MuellerMatrix m = W.at(pix); From 20377096bae07d4a4a8124f7fd5433c6d91af45d Mon Sep 17 00:00:00 2001 From: YL <20972148+LiYunyang@users.noreply.github.com> Date: Wed, 18 Mar 2026 20:37:30 -0500 Subject: [PATCH 5/8] more openmp support --- maps/src/maputils.cxx | 13 ++++++++++++- 1 file changed, 12 insertions(+), 1 deletion(-) diff --git a/maps/src/maputils.cxx b/maps/src/maputils.cxx index 24501e92..77b0d796 100644 --- a/maps/src/maputils.cxx +++ b/maps/src/maputils.cxx @@ -169,6 +169,9 @@ py::tuple GetRaDecMap(const G3SkyMap &m) ra->ConvertToDense(); dec->ConvertToDense(); + #ifdef _OPENMP + #pragma omp parallel for schedule(static) + #endif for (size_t i = 0; i < m.size(); i++) { std::vector radec = m.PixelToAngle(i); (*ra)[i] = radec[0]; @@ -203,6 +206,9 @@ G3SkyMapMaskPtr GetRaDecMask(const G3SkyMap &m, double ra_left, double ra_right, ra_left = wrap_ra(ra_left); ra_right = wrap_ra(ra_right); + #ifdef _OPENMP + #pragma omp parallel for schedule(static) + #endif for (size_t i = 0; i < m.size(); i++) { std::vector radec = m.PixelToAngle(i); double ra = wrap_ra(radec[0]); @@ -228,7 +234,9 @@ G3SkyMapMaskPtr GetGalacticPlaneMask(const G3SkyMap &m, double lat) if (m.coord_ref == MapCoordReference::Equatorial) { auto q_rot = get_fk5_j2000_to_gal_quat(); - + #ifdef _OPENMP + #pragma omp parallel for schedule(static) + #endif for (size_t i = 0; i < m.size(); i++) { // compute just latitude part of each coordinate auto q = q_rot * m.PixelToQuat(i); @@ -238,6 +246,9 @@ G3SkyMapMaskPtr GetGalacticPlaneMask(const G3SkyMap &m, double lat) (*mask)[i] = true; } } else if (m.coord_ref == MapCoordReference::Galactic) { + #ifdef _OPENMP + #pragma omp parallel for schedule(static) + #endif for (size_t i = 0; i < m.size(); i++) { auto q = m.PixelToQuat(i); if (fabs(q.d()) <= slat) From 2d8cb537ad5d1510195fafecca970aec1aae0f60 Mon Sep 17 00:00:00 2001 From: YL <20972148+LiYunyang@users.noreply.github.com> Date: Thu, 19 Mar 2026 11:57:18 -0500 Subject: [PATCH 6/8] [fix] remove omp for `RemoveWeights*` to avoid race condition in sparse format --- maps/src/maputils.cxx | 6 ------ 1 file changed, 6 deletions(-) diff --git a/maps/src/maputils.cxx b/maps/src/maputils.cxx index 77b0d796..9ec61d24 100644 --- a/maps/src/maputils.cxx +++ b/maps/src/maputils.cxx @@ -48,9 +48,6 @@ void RemoveWeights(G3SkyMap &T, G3SkyMap &Q, G3SkyMap &U, const G3SkyMapWeights v /= W.at(pix); } } else { - #ifdef _OPENMP - #pragma omp parallel for schedule(static) - #endif for (size_t pix = 0; pix < W.size(); pix++) { double t = T.at(pix); MuellerMatrix m = W.at(pix); @@ -94,9 +91,6 @@ void RemoveWeightsT(G3SkyMap &T, const G3SkyMapWeights &W, bool zero_nans) T.ConvertToDense(); T /= *(W.TT); } else { - #ifdef _OPENMP - #pragma omp parallel for schedule(static) - #endif for (size_t pix = 0; pix < W.size(); pix++) { double t = T.at(pix); MuellerMatrix m = W.at(pix); From eebc3293facee432f75f7b627dff2e9b8057337c Mon Sep 17 00:00:00 2001 From: Sasha Rahlin Date: Thu, 19 Mar 2026 13:09:19 -0500 Subject: [PATCH 7/8] Fix tabs, schedule(static) is default --- maps/src/maputils.cxx | 34 ++++++++++++++++++---------------- 1 file changed, 18 insertions(+), 16 deletions(-) diff --git a/maps/src/maputils.cxx b/maps/src/maputils.cxx index 9ec61d24..ac77159d 100644 --- a/maps/src/maputils.cxx +++ b/maps/src/maputils.cxx @@ -40,9 +40,10 @@ void RemoveWeights(G3SkyMap &T, G3SkyMap &Q, G3SkyMap &U, const G3SkyMapWeights T.ConvertToDense(); Q.ConvertToDense(); U.ConvertToDense(); + #ifdef _OPENMP - #pragma omp parallel for schedule(static) - #endif + #pragma omp parallel for + #endif for (size_t pix = 0; pix < T.size(); pix++) { StokesVector v(T[pix], Q[pix], U[pix]); v /= W.at(pix); @@ -127,9 +128,9 @@ void ApplyWeights(G3SkyMap &T, G3SkyMap &Q, G3SkyMap &U, const G3SkyMapWeights & g3_assert(!Q.weighted); g3_assert(!U.weighted); - #ifdef _OPENMP - #pragma omp parallel for schedule(static) - #endif + #ifdef _OPENMP + #pragma omp parallel for + #endif for (size_t pix = 0; pix < T.size(); pix++) { if (T.at(pix) == 0 && Q.at(pix) == 0 && U.at(pix) == 0) continue; @@ -163,9 +164,9 @@ py::tuple GetRaDecMap(const G3SkyMap &m) ra->ConvertToDense(); dec->ConvertToDense(); - #ifdef _OPENMP - #pragma omp parallel for schedule(static) - #endif + #ifdef _OPENMP + #pragma omp parallel for + #endif for (size_t i = 0; i < m.size(); i++) { std::vector radec = m.PixelToAngle(i); (*ra)[i] = radec[0]; @@ -200,9 +201,9 @@ G3SkyMapMaskPtr GetRaDecMask(const G3SkyMap &m, double ra_left, double ra_right, ra_left = wrap_ra(ra_left); ra_right = wrap_ra(ra_right); - #ifdef _OPENMP - #pragma omp parallel for schedule(static) - #endif + #ifdef _OPENMP + #pragma omp parallel for + #endif for (size_t i = 0; i < m.size(); i++) { std::vector radec = m.PixelToAngle(i); double ra = wrap_ra(radec[0]); @@ -228,9 +229,10 @@ G3SkyMapMaskPtr GetGalacticPlaneMask(const G3SkyMap &m, double lat) if (m.coord_ref == MapCoordReference::Equatorial) { auto q_rot = get_fk5_j2000_to_gal_quat(); - #ifdef _OPENMP - #pragma omp parallel for schedule(static) - #endif + + #ifdef _OPENMP + #pragma omp parallel for + #endif for (size_t i = 0; i < m.size(); i++) { // compute just latitude part of each coordinate auto q = q_rot * m.PixelToQuat(i); @@ -241,8 +243,8 @@ G3SkyMapMaskPtr GetGalacticPlaneMask(const G3SkyMap &m, double lat) } } else if (m.coord_ref == MapCoordReference::Galactic) { #ifdef _OPENMP - #pragma omp parallel for schedule(static) - #endif + #pragma omp parallel for + #endif for (size_t i = 0; i < m.size(); i++) { auto q = m.PixelToQuat(i); if (fabs(q.d()) <= slat) From 7941b7988e395f07a0ea9e09e4ab06bc35ce40b7 Mon Sep 17 00:00:00 2001 From: Sasha Rahlin Date: Thu, 19 Mar 2026 13:10:49 -0500 Subject: [PATCH 8/8] Parallelize common dense map operations --- maps/src/mapdata.cxx | 54 ++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 54 insertions(+) diff --git a/maps/src/mapdata.cxx b/maps/src/mapdata.cxx index 91dabe9b..e6459e95 100644 --- a/maps/src/mapdata.cxx +++ b/maps/src/mapdata.cxx @@ -1,4 +1,7 @@ #include "mapdata.h" +#ifdef _OPENMP +#include +#endif template typename SparseMapData::const_iterator @@ -46,6 +49,9 @@ DenseMapData * SparseMapData::to_dense() const { DenseMapData *rv = new DenseMapData(xlen_, ylen_); + #ifdef _OPENMP + #pragma omp parallel for + #endif for (size_t ix = 0; ix < data_.size(); ix++) { size_t x = offset_ + ix; const data_element &column = data_[ix]; @@ -190,6 +196,9 @@ SparseMapData::operator*=(const SparseMapData &r) assert(xlen_ == r.xlen_); assert(ylen_ == r.ylen_); + #ifdef _OPENMP + #pragma omp parallel for + #endif for (size_t ix = 0; ix < data_.size(); ix++) { size_t x = offset_ + ix; data_element &column = data_[ix]; @@ -209,6 +218,9 @@ SparseMapData::operator*=(const DenseMapData &r) assert(xlen_ == r.xdim()); assert(ylen_ == r.ydim()); + #ifdef _OPENMP + #pragma omp parallel for + #endif for (size_t ix = 0; ix < data_.size(); ix++) { size_t x = offset_ + ix; data_element &column = data_[ix]; @@ -225,6 +237,9 @@ template <> SparseMapData & SparseMapData::operator*=(double r) { + #ifdef _OPENMP + #pragma omp parallel for + #endif for (size_t ix = 0; ix < data_.size(); ix++) { data_element &column = data_[ix]; for (size_t iy = 0; iy < column.second.size(); iy++) { @@ -283,6 +298,9 @@ SparseMapData::operator/=(double r) { assert(r != 0); + #ifdef _OPENMP + #pragma omp parallel for + #endif for (size_t ix = 0; ix < data_.size(); ix++) { data_element &column = data_[ix]; for (size_t iy = 0; iy < column.second.size(); iy++) { @@ -299,6 +317,9 @@ DenseMapData::operator+=(const DenseMapData &r) assert(xlen_ == r.xlen_); assert(ylen_ == r.ylen_); + #ifdef _OPENMP + #pragma omp parallel for + #endif for (size_t x = 0; x < xlen_; x++) { for (size_t y = 0; y < ylen_; y++) { (*this)(x, y) += r.at(x, y); @@ -314,6 +335,9 @@ DenseMapData::operator+=(const SparseMapData &r) assert(xlen_ == r.xdim()); assert(ylen_ == r.ydim()); + #ifdef _OPENMP + #pragma omp parallel for + #endif for (size_t x = 0; x < xlen_; x++) { for (size_t y = 0; y < ylen_; y++) { (*this)(x, y) += r.at(x, y); @@ -329,6 +353,9 @@ DenseMapData::operator+=(double r) if (r == 0) return *this; + #ifdef _OPENMP + #pragma omp parallel for + #endif for (size_t x = 0; x < xlen_; x++) { for (size_t y = 0; y < ylen_; y++) { (*this)(x, y) += r; @@ -344,6 +371,9 @@ DenseMapData::operator-=(const DenseMapData &r) assert(xlen_ == r.xlen_); assert(ylen_ == r.ylen_); + #ifdef _OPENMP + #pragma omp parallel for + #endif for (size_t x = 0; x < xlen_; x++) { for (size_t y = 0; y < ylen_; y++) { (*this)(x, y) -= r.at(x, y); @@ -359,6 +389,9 @@ DenseMapData::operator-=(const SparseMapData &r) assert(xlen_ == r.xdim()); assert(ylen_ == r.ydim()); + #ifdef _OPENMP + #pragma omp parallel for + #endif for (size_t x = 0; x < xlen_; x++) { for (size_t y = 0; y < ylen_; y++) { (*this)(x, y) -= r.at(x, y); @@ -374,6 +407,9 @@ DenseMapData::operator-=(double r) if (r == 0) return *this; + #ifdef _OPENMP + #pragma omp parallel for + #endif for (size_t x = 0; x < xlen_; x++) { for (size_t y = 0; y < ylen_; y++) { (*this)(x, y) -= r; @@ -389,6 +425,9 @@ DenseMapData::operator*=(const DenseMapData &r) assert(xlen_ == r.xlen_); assert(ylen_ == r.ylen_); + #ifdef _OPENMP + #pragma omp parallel for + #endif for (size_t x = 0; x < xlen_; x++) { for (size_t y = 0; y < ylen_; y++) { (*this)(x, y) *= r.at(x, y); @@ -404,6 +443,9 @@ DenseMapData::operator*=(const SparseMapData &r) assert(xlen_ == r.xdim()); assert(ylen_ == r.ydim()); + #ifdef _OPENMP + #pragma omp parallel for + #endif for (size_t x = 0; x < xlen_; x++) { for (size_t y = 0; y < ylen_; y++) { (*this)(x, y) *= r.at(x, y); @@ -416,6 +458,9 @@ DenseMapData::operator*=(const SparseMapData &r) DenseMapData & DenseMapData::operator*=(double r) { + #ifdef _OPENMP + #pragma omp parallel for + #endif for (size_t x = 0; x < xlen_; x++) { for (size_t y = 0; y < ylen_; y++) { (*this)(x, y) *= r; @@ -431,6 +476,9 @@ DenseMapData::operator/=(const DenseMapData &r) assert(xlen_ == r.xlen_); assert(ylen_ == r.ylen_); + #ifdef _OPENMP + #pragma omp parallel for + #endif for (size_t x = 0; x < xlen_; x++) { for (size_t y = 0; y < ylen_; y++) { (*this)(x, y) /= r.at(x, y); @@ -446,6 +494,9 @@ DenseMapData::operator/=(const SparseMapData &r) assert(xlen_ == r.xdim()); assert(ylen_ == r.ydim()); + #ifdef _OPENMP + #pragma omp parallel for + #endif for (size_t x = 0; x < xlen_; x++) { for (size_t y = 0; y < ylen_; y++) { (*this)(x, y) /= r.at(x, y); @@ -458,6 +509,9 @@ DenseMapData::operator/=(const SparseMapData &r) DenseMapData & DenseMapData::operator/=(double r) { + #ifdef _OPENMP + #pragma omp parallel for + #endif for (size_t x = 0; x < xlen_; x++) { for (size_t y = 0; y < ylen_; y++) { (*this)(x, y) /= r;