Skip to content

Multithreaded sky map operations #219

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
arahlin merged 12 commits into from
Mar 19, 2026
Merged

Multithreaded sky map operations #219

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
12 commits
Select commit Hold shift + click to select a range
f11bada
[feat] multithread `RemoveWeights`
LiYunyang Jan 31, 2026
02db487
multithread `ApplyWeights`
LiYunyang Feb 2, 2026
f838b24
multithread `ApplyWeights`
LiYunyang Feb 2, 2026
52b8164
multithread `RemoveWeightsT`
LiYunyang Feb 4, 2026
82dc643
Merge remote-tracking branch 'origin/multithread-remove-weights' into…
LiYunyang Feb 4, 2026
9b46e84
Merge branch 'master' into multithread-remove-weights
arahlin Mar 18, 2026
2037709
more openmp support
LiYunyang Mar 19, 2026
2d8cb53
[fix] remove omp for `RemoveWeights*` to avoid race condition in spar…
LiYunyang Mar 19, 2026
eebc329
Fix tabs, schedule(static) is default
arahlin Mar 19, 2026
7941b79
Parallelize common dense map operations
arahlin Mar 19, 2026
670ae56
Merge branch 'master' into multithread-remove-weights
arahlin Mar 19, 2026
ae5466a
Merge branch 'master' into multithread-remove-weights
arahlin Mar 19, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
54 changes: 54 additions & 0 deletions maps/src/mapdata.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,4 +1,7 @@
#include "mapdata.h"
#ifdef _OPENMP
#include <omp.h>
#endif

template <typename T>
typename SparseMapData<T>::const_iterator
Expand Down Expand Up @@ -46,6 +49,9 @@ DenseMapData *
SparseMapData<T>::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];
Expand Down Expand Up @@ -190,6 +196,9 @@ SparseMapData<double>::operator*=(const SparseMapData<double> &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];
Expand All @@ -209,6 +218,9 @@ SparseMapData<double>::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];
Expand All @@ -225,6 +237,9 @@ template <>
SparseMapData<double> &
SparseMapData<double>::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++) {
Expand Down Expand Up @@ -283,6 +298,9 @@ SparseMapData<double>::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++) {
Expand All @@ -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);
Expand All @@ -314,6 +335,9 @@ DenseMapData::operator+=(const SparseMapData<double> &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);
Expand All @@ -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;
Expand All @@ -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);
Expand All @@ -359,6 +389,9 @@ DenseMapData::operator-=(const SparseMapData<double> &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);
Expand All @@ -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;
Expand All @@ -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);
Expand All @@ -404,6 +443,9 @@ DenseMapData::operator*=(const SparseMapData<double> &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);
Expand All @@ -416,6 +458,9 @@ DenseMapData::operator*=(const SparseMapData<double> &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;
Expand All @@ -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);
Expand All @@ -446,6 +494,9 @@ DenseMapData::operator/=(const SparseMapData<double> &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);
Expand All @@ -458,6 +509,9 @@ DenseMapData::operator/=(const SparseMapData<double> &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;
Expand Down
23 changes: 18 additions & 5 deletions maps/src/maputils.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -40,17 +40,15 @@ 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);
}
} 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);
Expand Down Expand Up @@ -130,6 +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
#endif
for (size_t pix = 0; pix < T.size(); pix++) {
if (T.at(pix) == 0 && Q.at(pix) == 0 && U.at(pix) == 0)
continue;
Expand Down Expand Up @@ -163,6 +164,9 @@ py::tuple GetRaDecMap(const G3SkyMap &m)
ra->ConvertToDense();
dec->ConvertToDense();

#ifdef _OPENMP
#pragma omp parallel for
#endif
for (size_t i = 0; i < m.size(); i++) {
std::vector<double> radec = m.PixelToAngle(i);
(*ra)[i] = radec[0];
Expand Down Expand Up @@ -197,6 +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
#endif
for (size_t i = 0; i < m.size(); i++) {
std::vector<double> radec = m.PixelToAngle(i);
double ra = wrap_ra(radec[0]);
Expand All @@ -223,6 +230,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
#endif
for (size_t i = 0; i < m.size(); i++) {
// compute just latitude part of each coordinate
auto q = q_rot * m.PixelToQuat(i);
Expand All @@ -232,6 +242,9 @@ G3SkyMapMaskPtr GetGalacticPlaneMask(const G3SkyMap &m, double lat)
(*mask)[i] = true;
}
} else if (m.coord_ref == MapCoordReference::Galactic) {
#ifdef _OPENMP
#pragma omp parallel for
#endif
for (size_t i = 0; i < m.size(); i++) {
auto q = m.PixelToQuat(i);
if (fabs(q.d()) <= slat)
Expand Down
Loading