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 <omp.h>
+#endif
 
 template <typename T>
 typename SparseMapData<T>::const_iterator
@@ -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];
@@ -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];
@@ -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];
@@ -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++) {
@@ -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++) {
@@ -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<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);
@@ -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<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);
@@ -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<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);
@@ -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;
@@ -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<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);
@@ -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;
diff --git a/maps/src/maputils.cxx b/maps/src/maputils.cxx
index d4ddf0b5..936a0390 100644
--- a/maps/src/maputils.cxx
+++ b/maps/src/maputils.cxx
@@ -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);
@@ -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;
@@ -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];
@@ -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]);
@@ -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);
@@ -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)