Repository: hpc4cmb/toast (language: python)

Issue:
More OpenMP target offload work
Work in progress, opening a draft PR for comparison.


Discussion / hints:

Fix patch:
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 6d99d618f..50f36be96 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -94,7 +94,8 @@ if(OpenMP_CXX_FOUND)
             # and need to override the flags. NOTE:  The use of "nordc"
             # is painful and needed since nvc++ cannot (as of 22.5) build
             # standalone shared loadable modules.
-            list(APPEND omp_opts -mp=gpu -gpu=nordc -cuda -cudalibs)
+            #list(APPEND omp_opts -mp=gpu -gpu=nordc -cuda -cudalibs)
+            list(APPEND omp_opts -mp=gpu)
         else()
             list(APPEND omp_opts ${OpenMP_CXX_FLAGS})
         endif()
diff --git a/cmake/FindFFTW.cmake b/cmake/FindFFTW.cmake
index ca07d5df0..6c675a8c1 100644
--- a/cmake/FindFFTW.cmake
+++ b/cmake/FindFFTW.cmake
@@ -65,7 +65,7 @@ if( FFTW_ROOT )
     FFTW_DOUBLE_LIB
     NAMES "fftw3" libfftw3-3
     PATHS ${FFTW_ROOT}
-    PATH_SUFFIXES "lib" "lib64"
+    PATH_SUFFIXES "lib" "lib64" "lib/x86_64-linux-gnu"
     NO_DEFAULT_PATH
   )
 
@@ -73,7 +73,7 @@ if( FFTW_ROOT )
     FFTW_DOUBLE_THREADS_LIB
     NAMES "fftw3_threads"
     PATHS ${FFTW_ROOT}
-    PATH_SUFFIXES "lib" "lib64"
+    PATH_SUFFIXES "lib" "lib64" "lib/x86_64-linux-gnu"
     NO_DEFAULT_PATH
   )
 
@@ -81,7 +81,7 @@ if( FFTW_ROOT )
     FFTW_DOUBLE_OPENMP_LIB
     NAMES "fftw3_omp"
     PATHS ${FFTW_ROOT}
-    PATH_SUFFIXES "lib" "lib64"
+    PATH_SUFFIXES "lib" "lib64" "lib/x86_64-linux-gnu"
     NO_DEFAULT_PATH
   )
 
@@ -89,7 +89,7 @@ if( FFTW_ROOT )
     FFTW_FLOAT_LIB
     NAMES "fftw3f" libfftw3f-3
     PATHS ${FFTW_ROOT}
-    PATH_SUFFIXES "lib" "lib64"
+    PATH_SUFFIXES "lib" "lib64" "lib/x86_64-linux-gnu"
     NO_DEFAULT_PATH
   )
 
@@ -97,7 +97,7 @@ if( FFTW_ROOT )
     FFTW_FLOAT_THREADS_LIB
     NAMES "fftw3f_threads"
     PATHS ${FFTW_ROOT}
-    PATH_SUFFIXES "lib" "lib64"
+    PATH_SUFFIXES "lib" "lib64" "lib/x86_64-linux-gnu"
     NO_DEFAULT_PATH
   )
 
@@ -105,7 +105,7 @@ if( FFTW_ROOT )
     FFTW_FLOAT_OPENMP_LIB
     NAMES "fftw3f_omp"
     PATHS ${FFTW_ROOT}
-    PATH_SUFFIXES "lib" "lib64"
+    PATH_SUFFIXES "lib" "lib64" "lib/x86_64-linux-gnu"
     NO_DEFAULT_PATH
   )
 
@@ -113,7 +113,7 @@ if( FFTW_ROOT )
     FFTW_LONGDOUBLE_LIB
     NAMES "fftw3l" libfftw3l-3
     PATHS ${FFTW_ROOT}
-    PATH_SUFFIXES "lib" "lib64"
+    PATH_SUFFIXES "lib" "lib64" "lib/x86_64-linux-gnu"
     NO_DEFAULT_PATH
   )
 
@@ -121,7 +121,7 @@ if( FFTW_ROOT )
     FFTW_LONGDOUBLE_THREADS_LIB
     NAMES "fftw3l_threads"
     PATHS ${FFTW_ROOT}
-    PATH_SUFFIXES "lib" "lib64"
+    PATH_SUFFIXES "lib" "lib64" "lib/x86_64-linux-gnu"
     NO_DEFAULT_PATH
   )
 
@@ -129,7 +129,7 @@ if( FFTW_ROOT )
     FFTW_LONGDOUBLE_OPENMP_LIB
     NAMES "fftw3l_omp"
     PATHS ${FFTW_ROOT}
-    PATH_SUFFIXES "lib" "lib64"
+    PATH_SUFFIXES "lib" "lib64" "lib/x86_64-linux-gnu"
     NO_DEFAULT_PATH
   )
 
diff --git a/etc/pyomptarget/Makefile b/etc/pyomptarget/Makefile
index 5c4924c46..c7498aceb 100644
--- a/etc/pyomptarget/Makefile
+++ b/etc/pyomptarget/Makefile
@@ -1,27 +1,40 @@
 
 # GCC 12.1.0 with offload support
-CXX := g++-12
-CXXFLAGS := -O0 -g -fPIC -pthread -std=c++11 -fcf-protection=none -fno-stack-protector
-CXX_OMP_FLAGS := -fopenmp -foffload=nvptx-none='-Wa,-m,sm_80 -misa=sm_80 -fPIC -lm -latomic'
+#CXX := g++-12
+#CXXFLAGS := -O0 -g -fPIC -pthread -std=c++11 -fcf-protection=none -fno-stack-protector
+#CXX_OMP_FLAGS := -fopenmp -foffload=nvptx-none='-Wa,-m,sm_80 -misa=sm_80 -fPIC -lm -latomic'
 
 # NVHPC
 # CXX = nvc++
 # CXXFLAGS = -O0 -g -fPIC -pthread -std=c++11
 # CXX_OMP_FLAGS = -mp=gpu -Minfo=mp -gpu=cc86
 
+# AMD LLVM
+CXX := amdclang++
+CXXFLAGS := -O3 -g -fPIC -std=c++11
+CXX_OMP_FLAGS := -fPIC -fopenmp -fopenmp-target-debug=3 -fopenmp-targets=amdgcn-amd-amdhsa -Xopenmp-target=amdgcn-amd-amdhsa -march=gfx1030
+CXX_LD_FLAGS := -lomp -lomptarget -lomptarget.rtl.amdgpu
 
 pybind := ../../src/toast/pybind11
 modext := $(shell python3-config --extension-suffix)
 pyincl := $(shell python3-config --includes)
 
 
+all : device_info test pyomptarget$(modext)
+
+
 pyomptarget$(modext) : module.o
-	$(CXX) -shared -o pyomptarget$(modext) module.o $(CXX_OMP_FLAGS)
+	$(CXX) -shared $(CXX_OMP_FLAGS) -o pyomptarget$(modext) module.o $(CXX_LD_FLAGS)
 
 module.o : module.cpp
 	$(CXX) $(CXXFLAGS) $(CXX_OMP_FLAGS) -I. -I$(pybind)/include $(pyincl) -c -o module.o module.cpp
 
+device_info : device_info.cpp
+	$(CXX) $(CXXFLAGS) $(CXX_OMP_FLAGS) -o device_info device_info.cpp $(CXX_LD_FLAGS)
+
+test : test.cpp
+	$(CXX) $(CXXFLAGS) $(CXX_OMP_FLAGS) -o test test.cpp $(CXX_LD_FLAGS)
 
 clean :
-	rm -f *.so *.o
+	rm -f *.so *.o device_info test
 
diff --git a/etc/pyomptarget/device_info.cpp b/etc/pyomptarget/device_info.cpp
new file mode 100644
index 000000000..a8d57d1f2
--- /dev/null
+++ b/etc/pyomptarget/device_info.cpp
@@ -0,0 +1,17 @@
+
+#include <omp.h>
+//#include "omptarget.h"
+#include <iostream>
+  
+int main(int argc, char **argv) {
+  int ndev = omp_get_num_devices();
+  std::cout << "OMP found " << ndev << " available target offload devices" << std::endl;
+  int target = ndev - 1;
+  int host = omp_get_initial_device();
+  int defdev = omp_get_default_device();
+  std::cout << "OMP initial host device = " << host << std::endl;
+  std::cout << "OMP target device = " << target << std::endl;
+  std::cout << "OMP default device = " << defdev << std::endl;
+  return 0;
+}
+
diff --git a/etc/pyomptarget/module.cpp b/etc/pyomptarget/module.cpp
index daacdf9a5..4c9f39e1b 100644
--- a/etc/pyomptarget/module.cpp
+++ b/etc/pyomptarget/module.cpp
@@ -23,9 +23,23 @@ namespace py = pybind11;
 
 using size_container = py::detail::any_container <ssize_t>;
 
+// 2/PI
+#define TWOINVPI 0.63661977236758134308
+
+// 2/3
+#define TWOTHIRDS 0.66666666666666666667
+
+// Helper table initialization
+void hpix_init_utab(uint64_t * utab) {
+    for (uint64_t m = 0; m < 256; ++m) {
+        utab[m] = (m & 0x1) | ((m & 0x2) << 1) | ((m & 0x4) << 2) |
+            ((m & 0x8) << 3) | ((m & 0x10) << 4) | ((m & 0x20) << 5) |
+            ((m & 0x40) << 6) | ((m & 0x80) << 7);
+    }
+    return;
+}
 
 // Device code
-
 #pragma omp declare target
 
 struct Interval {
@@ -47,6 +61,31 @@ void qa_mult(double const * p, double const * q, double * r) {
     return;
 }
 
+void qa_rotate(double const * q_in, double const * v_in, double * v_out) {
+    // The input quaternion has already been normalized on the host.
+
+    double xw =  q_in[3] * q_in[0];
+    double yw =  q_in[3] * q_in[1];
+    double zw =  q_in[3] * q_in[2];
+    double x2 = -q_in[0] * q_in[0];
+    double xy =  q_in[0] * q_in[1];
+    double xz =  q_in[0] * q_in[2];
+    double y2 = -q_in[1] * q_in[1];
+    double yz =  q_in[1] * q_in[2];
+    double z2 = -q_in[2] * q_in[2];
+
+    v_out[0] = 2 * ((y2 + z2) * v_in[0] + (xy - zw) * v_in[1] +
+                    (yw + xz) * v_in[2]) + v_in[0];
+
+    v_out[1] = 2 * ((zw + xy) * v_in[0] + (x2 + z2) * v_in[1] +
+                    (yz - xw) * v_in[2]) + v_in[1];
+
+    v_out[2] = 2 * ((xz - yw) * v_in[0] + (xw + yz) * v_in[1] +
+                    (x2 + y2) * v_in[2]) + v_in[2];
+
+    return;
+}
+
 void pointing_detector_inner(
     int32_t const * q_index,
     uint8_t const * flags,
@@ -79,6 +118,148 @@ void pointing_detector_inner(
     return;
 }
 
+uint64_t hpix_xy2pix(uint64_t * utab, uint64_t x, uint64_t y) {
+    return utab[x & 0xff] | (utab[(x >> 8) & 0xff] << 16) |
+           (utab[(x >> 16) & 0xff] << 32) |
+           (utab[(x >> 24) & 0xff] << 48) |
+           (utab[y & 0xff] << 1) | (utab[(y >> 8) & 0xff] << 17) |
+           (utab[(y >> 16) & 0xff] << 33) |
+           (utab[(y >> 24) & 0xff] << 49);
+}
+
+void hpix_vec2zphi(double const * vec,
+                   double * phi, int * region, double * z,
+                   double * rtz) {
+    // region encodes BOTH the sign of Z and whether its
+    // absolute value is greater than 2/3.
+    (*z) = vec[2];
+    double za = fabs(*z);
+    int itemp = ((*z) > 0.0) ? 1 : -1;
+    (*region) = (za <= TWOTHIRDS) ? itemp : itemp + itemp;
+    (*rtz) = sqrt(3.0 * (1.0 - za));
+    (*phi) = atan2(vec[1], vec[0]);
+    return;
+}
+
+void hpix_zphi2nest(int64_t nside, int64_t factor, uint64_t * utab, 
+                    double phi, int region, double z,
+                    double rtz, int64_t * pix) {
+    double tt = (phi >= 0.0) ? phi * TWOINVPI : phi * TWOINVPI + 4.0;
+    int64_t x;
+    int64_t y;
+    double temp1;
+    double temp2;
+    int64_t jp;
+    int64_t jm;
+    int64_t ifp;
+    int64_t ifm;
+    int64_t face;
+    int64_t ntt;
+    double tp;
+
+    double dnside = static_cast <double> (nside);
+    int64_t twonside = 2 * nside;
+    double halfnside = 0.5 * dnside;
+    double tqnside = 0.75 * dnside;
+    int64_t nsideminusone = nside - 1;
+
+    if ((region == 1) || (region == -1)) {
+        temp1 = halfnside + dnside * tt;
+        temp2 = tqnside * z;
+
+        jp = (int64_t)(temp1 - temp2);
+        jm = (int64_t)(temp1 + temp2);
+
+        ifp = jp >> factor;
+        ifm = jm >> factor;
+
+        if (ifp == ifm) {
+            face = (ifp == 4) ? (int64_t)4 : ifp + 4;
+        } else if (ifp < ifm) {
+            face = ifp;
+        } else {
+            face = ifm + 8;
+        }
+
+        x = jm & nsideminusone;
+        y = nsideminusone - (jp & nsideminusone);
+    } else {
+        ntt = (int64_t)tt;
+
+        tp = tt - (double)ntt;
+
+        temp1 = dnside * rtz;
+
+        jp = (int64_t)(tp * temp1);
+        jm = (int64_t)((1.0 - tp) * temp1);
+
+        if (jp >= nside) {
+            jp = nsideminusone;
+        }
+        if (jm >= nside) {
+            jm = nsideminusone;
+        }
+
+        if (z >= 0) {
+            face = ntt;
+            x = nsideminusone - jm;
+            y = nsideminusone - jp;
+        } else {
+            face = ntt + 8;
+            x = jp;
+            y = jm;
+        }
+    }
+
+    uint64_t sipf = hpix_xy2pix(utab, (uint64_t)x, (uint64_t)y);
+
+    (*pix) = (int64_t)sipf + (face << (2 * factor));
+
+    return;
+}
+
+void pixels_healpix_nest_inner(
+    int64_t nside,
+    int64_t factor,
+    uint64_t * utab,
+    int32_t const * quat_index,
+    int32_t const * pixel_index,
+    double const * quats,
+    uint8_t const * flags,
+    uint8_t * hsub,
+    int64_t * pixels,
+    int64_t n_pix_submap,
+    int64_t isamp,
+    int64_t n_samp,
+    int64_t idet
+) {
+    const double zaxis[3] = {0.0, 0.0, 1.0};
+    int32_t p_indx = pixel_index[idet];
+    int32_t q_indx = quat_index[idet];
+    double dir[3];
+    double z;
+    double rtz;
+    double phi;
+    int region;
+    size_t qoff = (q_indx * 4 * n_samp) + 4 * isamp;
+    size_t poff = p_indx * n_samp + isamp;
+    int64_t sub_map;
+
+    uint8_t check = flags[isamp] & 255;
+
+    if (check != 0) {
+        pixels[poff] = -1;
+    } else {
+        qa_rotate(&(quats[qoff]), zaxis, dir);
+        hpix_vec2zphi(dir, &phi, &region, &z, &rtz);
+        hpix_zphi2nest(nside, factor, utab, phi, region, z, rtz, &(pixels[poff]));
+        sub_map = (int64_t)(pixels[poff] / n_pix_submap);
+        hsub[sub_map] = 1;
+    }
+
+    return;
+}
+
 #pragma omp end declare target
 
 template <typename T>
@@ -296,7 +477,8 @@ PYBIND11_MODULE(pyomptarget, m) {
             py::buffer boresight,
             py::buffer quats,
             py::buffer intervals,
-            py::buffer shared_flags
+            py::buffer shared_flags,
+            py::buffer pixels
         ) {
             int ndev = omp_get_num_devices();
             std::cout << "OMP found " << ndev << " available target offload devices" << std::endl;
@@ -326,6 +508,10 @@ PYBIND11_MODULE(pyomptarget, m) {
             );
             int64_t n_det = temp_shape[0];
 
+            int64_t * raw_pixels = extract_buffer <int64_t> (
+                pixels, "pixels", 2, temp_shape, {n_det, n_samp}
+            );
+
             Interval * raw_intervals = extract_buffer <Interval> (
                 intervals, "intervals", 1, temp_shape, {-1}
             );
@@ -335,6 +521,8 @@ PYBIND11_MODULE(pyomptarget, m) {
 
             host_to_device(mem, n_samp, raw_flags, "flags");
 
+            host_to_device(mem, n_samp * n_det, raw_pixels, "pixels");
+
             host_to_device(mem, n_view, raw_intervals, "intervals");
 
             host_to_device(mem, n_samp * n_det * 4, raw_quats, "quats");
@@ -345,7 +533,8 @@ PYBIND11_MODULE(pyomptarget, m) {
     m.def(
         "unstage_data", [](
             std::unordered_map <void *, void *> mem,
-            py::buffer quats
+            py::buffer quats,
+            py::buffer pixels
         ) {
             std::vector <int64_t> temp_shape(3);
             double * raw_quats = extract_buffer <double> (
@@ -353,7 +542,11 @@ PYBIND11_MODULE(pyomptarget, m) {
             );
             int64_t n_det = temp_shape[0];
             int64_t n_samp = temp_shape[1];
+            int64_t * raw_pixels = extract_buffer <int64_t> (
+                pixels, "pixels", 2, temp_shape, {n_det, n_samp}
+            );
             device_to_host(mem, n_samp * n_det * 4, raw_quats, "quats");
+            device_to_host(mem, n_samp * n_det, raw_pixels, "pixels");
 
             int target = omp_get_num_devices() - 1;
             for(auto & p : mem) {
@@ -447,6 +640,9 @@ PYBIND11_MODULE(pyomptarget, m) {
                             isamp <= raw_intervals[iview].last;
                             isamp++
                         ) {
+                            // Test stack allocation with a mapped variable
+                            //double dummy[n_det];
+
                             pointing_detector_inner(
                                 raw_quat_index,
                                 dev_flags,
@@ -465,4 +661,127 @@ PYBIND11_MODULE(pyomptarget, m) {
             return;
         });
 
+    m.def(
+        "pixels_healpix_nest", [](
+            std::unordered_map <void *, void *> mem,
+            py::buffer quat_index,
+            py::buffer quats,
+            py::buffer shared_flags,
+            py::buffer pixel_index,
+            py::buffer pixels,
+            py::buffer intervals,
+            py::buffer hit_submaps,
+            int64_t n_pix_submap,
+            int64_t nside
+        ) {
+            // This is used to return the actual shape of each buffer
+            std::vector <int64_t> temp_shape(3);
+
+            int32_t * raw_quat_index = extract_buffer <int32_t> (
+                quat_index, "quat_index", 1, temp_shape, {-1}
+            );
+            int64_t n_det = temp_shape[0];
+
+            int32_t * raw_pixel_index = extract_buffer <int32_t> (
+                pixel_index, "pixel_index", 1, temp_shape, {n_det}
+            );
+
+            int64_t * raw_pixels = extract_buffer <int64_t> (
+                pixels, "pixels", 2, temp_shape, {-1, -1}
+            );
+            int64_t n_samp = temp_shape[1];
+
+            double * raw_quats = extract_buffer <double> (
+                quats, "quats", 3, temp_shape, {-1, n_samp, 4}
+            );
+
+            Interval * raw_intervals = extract_buffer <Interval> (
+                intervals, "intervals", 1, temp_shape, {-1}
+            );
+            int64_t n_view = temp_shape[0];
+
+            uint8_t * raw_hsub = extract_buffer <uint8_t> (
+                hit_submaps, "hit_submaps", 1, temp_shape, {-1}
+            );
+            int64_t n_submap = temp_shape[0];
+
+            uint8_t * raw_flags = extract_buffer <uint8_t> (
+                shared_flags, "flags", 1, temp_shape, {n_samp}
+            );
+
+            uint64_t utab[256];
+            hpix_init_utab(utab);
+
+            int64_t factor = 0;
+            while (nside != (1ll << factor)) {
+                ++factor;
+            }
+
+            double * dev_quats = (double*)mem.at(raw_quats);
+            
+            int64_t * dev_pixels = (int64_t*)mem.at(raw_pixels);
+
+            Interval * dev_intervals = (Interval*)mem.at(raw_intervals);
+
+            uint8_t * dev_flags = (uint8_t*)mem.at(raw_flags);
+
+            // Make sure the lookup table exists on device
+            // if (mem.count(vutab) == 0) {
+            //     host_to_device(mem, 256, utab, "hpix_utab");
+            // }
+            // uint64_t * dev_utab = (uint64_t*)mem.at(utab);
+
+            # pragma omp target data  \
+            device(0)               \
+            map(to:                   \
+            utab[0:256],              \
+            raw_pixel_index[0:n_det], \
+            raw_quat_index[0:n_det],  \
+            n_pix_submap,             \
+            nside,                    \
+            factor,                   \
+            n_view,                   \
+            n_det,                    \
+            n_samp                   \
+            )                         \
+            map(tofrom: raw_hsub[0:n_submap])
+            {
+                # pragma omp target teams distribute collapse(2) \
+                is_device_ptr(                                   \
+                dev_pixels,                                      \
+                dev_quats,                                       \
+                dev_flags,                                       \
+                dev_intervals                                   \
+                )
+                for (int64_t idet = 0; idet < n_det; idet++) {
+                    for (int64_t iview = 0; iview < n_view; iview++) {
+                        # pragma omp parallel for default(shared)
+                        for (
+                            int64_t isamp = dev_intervals[iview].first;
+                            isamp <= dev_intervals[iview].last;
+                            isamp++
+                        ) {
+                            pixels_healpix_nest_inner(
+                                nside,
+                                factor,
+                                utab,
+                                raw_quat_index,
+                                raw_pixel_index,
+                                dev_quats,
+                                dev_flags,
+                                raw_hsub,
+                                dev_pixels,
+                                n_pix_submap,
+                                isamp,
+                                n_samp,
+                                idet
+                            );
+                        }
+                    }
+                }
+            }
+
+            return;
+        });
+
 }
diff --git a/etc/pyomptarget/run.py b/etc/pyomptarget/run.py
index 5e51908b9..3f4bc6660 100644
--- a/etc/pyomptarget/run.py
+++ b/etc/pyomptarget/run.py
@@ -52,6 +52,7 @@ def main():
     ).reshape((-1, 4))
 
     quat_index = np.arange(ndet, dtype=np.int32)
+    pixel_index = np.arange(ndet, dtype=np.int32)
 
     shared_flags = np.zeros(nsamp, dtype=np.uint8)
 
@@ -64,8 +65,9 @@ def main():
     intervals[0].stop = -1.0
 
     quats = np.zeros((ndet, nsamp, 4), dtype=np.float64)
+    pixels = np.zeros((ndet, nsamp), dtype=np.int64)
 
-    mem = pyomptarget.stage_data(boresight, quats, intervals, shared_flags)
+    mem = pyomptarget.stage_data(boresight, quats, intervals, shared_flags, pixels)
 
     # time.sleep(5)
 
@@ -79,9 +81,29 @@ def main():
         shared_flags,
     )
 
-    pyomptarget.unstage_data(mem, quats)
+    nside = 1024
+    nside_submap = 16
+    n_pix_submap = 12 * nside_submap**2
+    n_submap = (nside // nside_submap) ** 2
+    hit_submaps = np.zeros(n_submap, dtype=np.uint8)
 
-    # print(quats)
+    pyomptarget.pixels_healpix_nest(
+        mem,
+        quat_index,
+        quats,
+        shared_flags,
+        pixel_index,
+        pixels,
+        intervals,
+        hit_submaps,
+        n_pix_submap,
+        nside,
+    )
+
+    pyomptarget.unstage_data(mem, quats, pixels)
+
+    print(quats)
+    print(pixels)
 
 
 if __name__ == "__main__":
diff --git a/etc/pyomptarget/run_rocm.sh b/etc/pyomptarget/run_rocm.sh
new file mode 100755
index 000000000..7bb4a5259
--- /dev/null
+++ b/etc/pyomptarget/run_rocm.sh
@@ -0,0 +1,14 @@
+#!/bin/bash
+
+export OMP_TARGET_OFFLOAD=MANDATORY
+export OMP_DISPLAY_ENV=verbose
+
+#export LIBOMPTARGET_INFO=$((0x1 | 0x2 | 0x4 | 0x8 | 0x10 | 0x20))
+export LIBOMPTARGET_INFO=-1
+export LIBOMPTARGET_DEVICE_RTL_DEBUG=3
+
+#export HSA_IGNORE_SRAMECC_MISREPORT=1
+
+eval "$@"
+
+# OMP_NUM_THREADS=4 OMP_TARGET_OFFLOAD=MANDATORY OMP_DISPLAY_ENV=verbose LIBOMPTARGET_INFO=-1 LIBOMPTARGET_DEVICE_RTL_DEBUG=3 TOAST_LOGLEVEL=VERBOSE TOAST_GPU_OPENMP=1 toast_mini --node_mem_gb 7 2>&1 | tee log
\ No newline at end of file
diff --git a/platforms/cmbenv-rocm-dbg.sh b/platforms/cmbenv-rocm-dbg.sh
new file mode 100755
index 000000000..e2aefbb42
--- /dev/null
+++ b/platforms/cmbenv-rocm-dbg.sh
@@ -0,0 +1,31 @@
+#!/bin/bash
+
+# Pass extra cmake options to this script, including
+# things like -DCMAKE_INSTALL_PREFIX=/path/to/install, etc.
+
+opts="$@"
+
+cmake \
+    -DCMAKE_C_COMPILER="amdclang" \
+    -DCMAKE_CXX_COMPILER="amdclang++" \
+    -DCMAKE_BUILD_TYPE=Debug \
+    -DCMAKE_C_FLAGS="-O0 -g -fPIC -Wno-error=implicit-function-declaration" \
+    -DCMAKE_CXX_FLAGS="-O0 -g -fPIC" \
+    -DUSE_OPENMP_TARGET=TRUE \
+    -DOPENMP_TARGET_FLAGS="\
+    -gline-tables-only \
+    -fopenmp-target-debug=3 \
+    -fopenmp-targets=amdgcn-amd-amdhsa \
+    -Xopenmp-target=amdgcn-amd-amdhsa \
+    -march=gfx1030 \
+    " \
+    -DPYTHON_EXECUTABLE:FILEPATH=$(which python3) \
+    -DCMAKE_VERBOSE_MAKEFILE:BOOL=ON \
+    -DAATM_ROOT="${CMBENV_AUX_ROOT}" \
+    -DFFTW_ROOT="${CMBENV_AUX_ROOT}" \
+    -DBLAS_LIBRARIES=${CMBENV_AUX_ROOT}/lib/libopenblas.so \
+    -DLAPACK_LIBRARIES=${CMBENV_AUX_ROOT}/lib/libopenblas.so \
+    -DSUITESPARSE_INCLUDE_DIR_HINTS="${CMBENV_AUX_ROOT}/include" \
+    -DSUITESPARSE_LIBRARY_DIR_HINTS="${CMBENV_AUX_ROOT}/lib" \
+    ${opts} \
+    ..
diff --git a/platforms/cmbenv-rocm.sh b/platforms/cmbenv-rocm.sh
new file mode 100755
index 000000000..218af73d1
--- /dev/null
+++ b/platforms/cmbenv-rocm.sh
@@ -0,0 +1,28 @@
+#!/bin/bash
+
+# Pass extra cmake options to this script, including
+# things like -DCMAKE_INSTALL_PREFIX=/path/to/install, etc.
+
+opts="$@"
+
+cmake \
+    -DCMAKE_C_COMPILER="amdclang" \
+    -DCMAKE_CXX_COMPILER="amdclang++" \
+    -DCMAKE_C_FLAGS="-O3 -g -fPIC -Wno-error=implicit-function-declaration" \
+    -DCMAKE_CXX_FLAGS="-O3 -g -fPIC" \
+    -DUSE_OPENMP_TARGET=TRUE \
+    -DOPENMP_TARGET_FLAGS="\
+    -fopenmp-targets=amdgcn-amd-amdhsa \
+    -Xopenmp-target=amdgcn-amd-amdhsa \
+    -march=gfx1030 \
+    " \
+    -DPYTHON_EXECUTABLE:FILEPATH=$(which python3) \
+    -DCMAKE_VERBOSE_MAKEFILE:BOOL=ON \
+    -DAATM_ROOT="${CMBENV_AUX_ROOT}" \
+    -DFFTW_ROOT="${CMBENV_AUX_ROOT}" \
+    -DBLAS_LIBRARIES=${CMBENV_AUX_ROOT}/lib/libopenblas.so \
+    -DLAPACK_LIBRARIES=${CMBENV_AUX_ROOT}/lib/libopenblas.so \
+    -DSUITESPARSE_INCLUDE_DIR_HINTS="${CMBENV_AUX_ROOT}/include" \
+    -DSUITESPARSE_LIBRARY_DIR_HINTS="${CMBENV_AUX_ROOT}/lib" \
+    ${opts} \
+    ..
diff --git a/platforms/linux-llvm.sh b/platforms/linux-llvm.sh
new file mode 100755
index 000000000..ff22c5092
--- /dev/null
+++ b/platforms/linux-llvm.sh
@@ -0,0 +1,18 @@
+#!/bin/bash
+
+# Pass extra cmake options to this script, including
+# things like -DCMAKE_INSTALL_PREFIX=/path/to/install, etc.
+
+opts="$@"
+
+cmake \
+    -DCMAKE_C_COMPILER="clang" \
+    -DCMAKE_CXX_COMPILER="clang++" \
+    -DCMAKE_C_FLAGS="-O3 -g -fPIC -pthread" \
+    -DCMAKE_CXX_FLAGS="-O3 -g -fPIC -pthread -std=c++11" \
+    -DOPENMP_TARGET_FLAGS="-fopenmp -fopenmp-targets=nvptx64 -fopenmp-target-debug=3 --offload-arch=sm_80 " \
+    -DPYTHON_EXECUTABLE:FILEPATH=$(which python3) \
+    -DUSE_OPENMP_TARGET=TRUE \
+    -DCMAKE_VERBOSE_MAKEFILE:BOOL=ON \
+    ${opts} \
+    ..
diff --git a/platforms/linux-nvhpc-dbg.sh b/platforms/linux-nvhpc-dbg.sh
deleted file mode 100755
index 43f26837f..000000000
--- a/platforms/linux-nvhpc-dbg.sh
+++ /dev/null
@@ -1,19 +0,0 @@
-#!/bin/bash
-
-# Pass extra cmake options to this script, including
-# things like -DCMAKE_INSTALL_PREFIX=/path/to/install, etc.
-
-opts="$@"
-
-cmake \
-    -DCMAKE_BUILD_TYPE="Debug" \
-    -DCMAKE_C_COMPILER="nvc" \
-    -DCMAKE_CXX_COMPILER="nvc++" \
-    -DCMAKE_C_FLAGS="-O0 -g -fPIC -pthread" \
-    -DCMAKE_CXX_FLAGS="-O0 -g -fPIC -pthread -std=c++11" \
-    -DOPENMP_TARGET_FLAGS="-Minfo=mp" \
-    -DTOAST_STATIC_DEPS:BOOL=ON \
-    -DPYTHON_EXECUTABLE:FILEPATH=$(which python3) \
-    -DCMAKE_VERBOSE_MAKEFILE:BOOL=ON \
-    ${opts} \
-    ..
diff --git a/platforms/linux-nvhpc.sh b/platforms/linux-nvhpc.sh
index c2e9e3bcc..22dfbf900 100755
--- a/platforms/linux-nvhpc.sh
+++ b/platforms/linux-nvhpc.sh
@@ -5,14 +5,35 @@
 
 opts="$@"
 
+if [ "x${CMBENV_AUX_ROOT}" = "x" ]; then
+    echo "You must activate a cmbenv environment before using this script."
+    exit 1
+fi
+
+if [ "x${DEBUG}" = "x1" ]; then
+    CMAKE_BUILD_TYPE=Debug
+    oflags="-O0"
+else
+    CMAKE_BUILD_TYPE=Release
+    oflags="-O3"
+fi
+
 cmake \
+    -DCMAKE_BUILD_TYPE="${CMAKE_BUILD_TYPE}" \
     -DCMAKE_C_COMPILER="nvc" \
     -DCMAKE_CXX_COMPILER="nvc++" \
-    -DCMAKE_C_FLAGS="-O3 -g -fPIC -pthread" \
-    -DCMAKE_CXX_FLAGS="-O3 -g -fPIC -pthread -std=c++11" \
-    -DOPENMP_TARGET_FLAGS="-Minfo=mp" \
+    -DCMAKE_C_FLAGS="${oflags} -g -fPIC -pthread" \
+    -DCMAKE_CXX_FLAGS="${oflags} -g -fPIC -pthread -std=c++11" \
+    -DOPENMP_TARGET_FLAGS="-Minfo=mp -gpu=cc86" \
     -DPYTHON_EXECUTABLE:FILEPATH=$(which python3) \
     -DUSE_OPENMP_TARGET=TRUE \
-    -DCMAKE_VERBOSE_MAKEFILE:BOOL=ON \
-    ${opts} \
+    -DCMAKE_VERBOSE_MAKEFILE=1 \
+    -DCMAKE_INSTALL_PREFIX="${CMBENV_AUX_ROOT}" \
+    -DFFTW_ROOT="${CMBENV_AUX_ROOT}" \
+    -DDISABLE_FLAC=1 \
+    -DBLAS_LIBRARIES="-L${NVHPC_ROOT}/compilers/lib -lblas -lnvf -mp" \
+    -DLAPACK_LIBRARIES="-L${NVHPC_ROOT}/compilers/lib -llapack -lnvf -mp" \
+    -DSUITESPARSE_INCLUDE_DIR_HINTS="${CMBENV_AUX_ROOT}/include" \
+    -DSUITESPARSE_LIBRARY_DIR_HINTS="${CMBENV_AUX_ROOT}/lib" ${opts} \
     ..
+
diff --git a/src/libtoast/CMakeLists.txt b/src/libtoast/CMakeLists.txt
index 39d80bb29..3a5e5a1f9 100644
--- a/src/libtoast/CMakeLists.txt
+++ b/src/libtoast/CMakeLists.txt
@@ -164,26 +164,26 @@ endif()
 
 # Add the serial unit test executable
 
-add_executable(toast_test
-    toast_test.cpp
-)
+# add_executable(toast_test
+#     toast_test.cpp
+# )
 
-if(OpenMP_CXX_FOUND)
-    target_compile_options(toast_test PRIVATE "${OpenMP_CXX_FLAGS}")
-endif(OpenMP_CXX_FOUND)
+# if(OpenMP_CXX_FOUND)
+#     target_compile_options(toast_test PRIVATE "${OpenMP_CXX_FLAGS}")
+# endif(OpenMP_CXX_FOUND)
 
-target_include_directories(toast_test BEFORE PRIVATE
-    "${CMAKE_CURRENT_SOURCE_DIR}"
-    "${CMAKE_CURRENT_SOURCE_DIR}/include"
-    "${CMAKE_CURRENT_SOURCE_DIR}/tests"
-)
+# target_include_directories(toast_test BEFORE PRIVATE
+#     "${CMAKE_CURRENT_SOURCE_DIR}"
+#     "${CMAKE_CURRENT_SOURCE_DIR}/include"
+#     "${CMAKE_CURRENT_SOURCE_DIR}/tests"
+# )
 
-target_link_libraries(toast_test toast)
+# target_link_libraries(toast_test toast)
 
-target_link_options(toast_test PRIVATE
-    ${LAPACK_LINKER_FLAGS} ${BLAS_LINKER_FLAGS} ${OpenMP_CXX_FLAGS}
-)
+# target_link_options(toast_test PRIVATE
+#     ${LAPACK_LINKER_FLAGS} ${BLAS_LINKER_FLAGS} ${OpenMP_CXX_FLAGS}
+# )
 
-install(TARGETS toast_test DESTINATION ${CMAKE_INSTALL_BINDIR})
+# install(TARGETS toast_test DESTINATION ${CMAKE_INSTALL_BINDIR})
 
-add_test(NAME serial_tests COMMAND toast_test)
+# add_test(NAME serial_tests COMMAND toast_test)
diff --git a/src/toast/CMakeLists.txt b/src/toast/CMakeLists.txt
index e462806c3..8ed91cedd 100644
--- a/src/toast/CMakeLists.txt
+++ b/src/toast/CMakeLists.txt
@@ -55,6 +55,8 @@ pybind11_add_module(_libtoast MODULE
     _libtoast/ops_stokes_weights.cpp
     _libtoast/ops_pixels_healpix.cpp
     _libtoast/ops_mapmaker_utils.cpp
+    _libtoast/ops_noise_weight.cpp
+    _libtoast/ops_scan_map.cpp
 )
 
 # Handle recursion limit on NVC++ / PGI
diff --git a/src/toast/__init__.py b/src/toast/__init__.py
index dd59bfbf5..19ae022a0 100644
--- a/src/toast/__init__.py
+++ b/src/toast/__init__.py
@@ -26,6 +26,10 @@
       target offload.
     * Requires compile-time support for OpenMP 5.x features.
 
+TOAST_GPU_MEM_GB=<value>
+    * Value in GB of memory to use for OpenMP on each GPU.
+    * Conservative default is 2GB.
+
 TOAST_GPU_JAX=<value>
     * Values "1", "true", or "yes" will enable runtime support for jax.
     * Requires jax to be available / importable.
diff --git a/src/toast/_libtoast/accelerator.cpp b/src/toast/_libtoast/accelerator.cpp
index d26ff1e59..2888a7b02 100644
--- a/src/toast/_libtoast/accelerator.cpp
+++ b/src/toast/_libtoast/accelerator.cpp
@@ -10,12 +10,228 @@
 #include <utility>
 
 
+OmpPoolResource::OmpPoolResource(int target, size_t size, size_t align) {
+    if (target < 0) {
+        #ifdef HAVE_OPENMP_TARGET
+        target_ = omp_get_num_devices() - 1;
+        #endif
+    } else {
+        target_ = target;
+    }
+    pool_size_ = size;
+    alignment_ = align;
+    alloc();
+}
+
+OmpPoolResource::~OmpPoolResource() {
+    release();
+}
+
+bool OmpPoolResource::verbose() {
+    // Helper function to check the log level once and return a static value.
+    // This is useful to avoid expensive string operations inside deeply nested
+    // functions.
+    static bool called = false;
+    static bool verbose = false;
+    if (!called) {
+        // First time we were called
+        auto & env = toast::Environment::get();
+        std::string logval = env.log_level();
+        if (strncmp(logval.c_str(), "VERBOSE", 7) == 0) {
+            verbose = true;
+        }
+        called = true;
+    }
+    return verbose;
+}
+
+void OmpPoolResource::release() {
+    bool extra = verbose();
+
+    // Reset our block containers and free device memory.
+    blocks_.clear();
+    ptr_to_block_.clear();
+    if (raw_ != nullptr) {
+        if (extra) {
+            std::ostringstream o;
+            auto log = toast::Logger::get();
+            o << "OmpPoolResource:  Free " << pool_size_ << " bytes on device " <<
+                target_;
+            o << " at " << raw_;
+            log.verbose(o.str().c_str());
+        }
+        #ifdef HAVE_OPENMP_TARGET
+        omp_target_free(raw_, target_);
+        #endif
+    }
+    return;
+}
+
+void * OmpPoolResource::allocate(size_t bytes) {
+    bool extra = verbose();
+
+    // Compute the number of aligned bytes needed
+    size_t aligned_bytes = compute_aligned_bytes(bytes, alignment_);
+
+    // Check that we have available space
+    if (pool_used_ + aligned_bytes > pool_size_) {
+        std::ostringstream o;
+        auto log = toast::Logger::get();
+        o << "OmpPoolResource:  request of " << aligned_bytes << " aligned bytes ";
+        o << "would exceed pool size of " << pool_size_ << " bytes";
+        log.error(o.str().c_str());
+        throw std::runtime_error(o.str().c_str());
+    }
+
+    // Build a block starting at the end of the existing blocks
+    void * next_start = blocks_.back().end;
+    OmpBlock next_block;
+    next_block.start = next_start;
+    next_block.end = shift_void_pointer(next_start, aligned_bytes);
+    next_block.size = bytes;
+    next_block.aligned_size = aligned_bytes;
+    next_block.is_free = false;
+    if (extra) {
+        std::ostringstream o;
+        auto log = toast::Logger::get();
+        o << "OmpPoolResource:  Append block " << next_block.start << ":";
+        o << next_block.end << " with " << next_block.size << " bytes (";
+        o << aligned_bytes << " aligned)";
+        log.verbose(o.str().c_str());
+    }
+
+    // Store the block and its index
+    size_t bindex = blocks_.size();
+    ptr_to_block_[next_start] = bindex;
+    blocks_.push_back(next_block);
+
+    // Increment the used count
+    pool_used_ += aligned_bytes;
+
+    return next_start;
+}
+
+void OmpPoolResource::deallocate(void * p) {
+    bool extra = verbose();
+
+    // Get the memory block associated with this pointer
+    size_t bindex = block_index(p);
+    OmpBlock & block = blocks_[bindex];
+    if (extra) {
+        std::ostringstream o;
+        auto log = toast::Logger::get();
+        o << "OmpPoolResource:  Delete block " << block.start << ":";
+        o << block.end << " with " << block.size << " bytes (";
+        o << block.aligned_size << " aligned)";
+        log.verbose(o.str().c_str());
+    }
+
+    // Mark block as unused
+    block.is_free = true;
+
+    // Free trailing unused blocks
+    garbage_collection();
+}
+
+void OmpPoolResource::alloc() {
+    bool extra = verbose();
+
+    if (pool_size_ == 0) {
+        // Attempt to guess at the pool size
+        pool_size_ = 1024 * 1024 * 1024;
+        std::ostringstream o;
+        auto log = toast::Logger::get();
+        o << "OmpPoolResource:  Using " << pool_size_;
+        o << " bytes for default pool size";
+        log.warning(o.str().c_str());
+    }
+    #ifdef HAVE_OPENMP_TARGET
+    raw_ = omp_target_alloc(pool_size_, target_);
+    #endif
+    if (extra) {
+        std::ostringstream o;
+        auto log = toast::Logger::get();
+        o << "OmpPoolResource:  Allocated pool of " << pool_size_ << " bytes";
+        o << " on target " << target_;
+        log.verbose(o.str().c_str());
+    }
+
+    // Current bytes used
+    pool_used_ = 0;
+
+    // Large enough for any device?
+    alignment_ = 512;
+
+    // Initialize an empty starting block
+    OmpBlock first_block;
+    first_block.start = raw_;
+    first_block.end = raw_;
+    first_block.size = 0;
+    first_block.aligned_size = 0;
+    first_block.is_free = false;
+    blocks_.push_back(first_block);
+}
+
+size_t OmpPoolResource::block_index(void * ptr) {
+    auto entry = ptr_to_block_.find(ptr);
+    if (entry == ptr_to_block_.end()) {
+        std::ostringstream o;
+        auto log = toast::Logger::get();
+        o << "OmpPoolResource:  Pointer " << ptr << " is not mapped to a block";
+        log.error(o.str().c_str());
+        throw std::runtime_error(o.str().c_str());
+    }
+    return entry->second;
+}
+
+// Remove free blocks
+void OmpPoolResource::garbage_collection() {
+    bool extra = verbose();
+    OmpBlock & block = blocks_.back();
+    while (block.is_free) {
+        // Remove pointer from map
+        if (extra) {
+            std::ostringstream o;
+            auto log = toast::Logger::get();
+            o << "OmpPoolResource:  GC block " << block.start << ":";
+            o << block.end << " with " << block.size << " bytes (";
+            o << block.aligned_size << " aligned)";
+            log.verbose(o.str().c_str());
+        }
+        ptr_to_block_.erase(block.start);
+
+        // Shrink the currently used size
+        pool_used_ -= block.aligned_size;
+
+        // Delete block
+        blocks_.pop_back();
+
+        // Get next block
+        block = blocks_.back();
+    }
+}
+
+size_t OmpPoolResource::compute_aligned_bytes(size_t bytes, size_t alignment) {
+    size_t aligned_blocks = (size_t)(bytes / alignment);
+    if (aligned_blocks * alignment != bytes) {
+        aligned_blocks += 1;
+    }
+    return aligned_blocks * alignment;
+}
+
+void * OmpPoolResource::shift_void_pointer(void * ptr, size_t offset_bytes) {
+    uint8_t * bptr = static_cast <uint8_t *> (ptr);
+    void * new_ptr = static_cast <void *> (&(bptr[offset_bytes]));
+    return new_ptr;
+}
+
 OmpManager & OmpManager::get() {
     static OmpManager instance;
     return instance;
 }
 
-void OmpManager::assign_device(int node_procs, int node_rank, bool disabled) {
+void OmpManager::assign_device(int node_procs, int node_rank, float mem_gb,
+                               bool disabled) {
     std::ostringstream o;
     auto log = toast::Logger::get();
     if ((node_procs < 1) || (node_rank < 0)) {
@@ -32,7 +248,6 @@ void OmpManager::assign_device(int node_procs, int node_rank, bool disabled) {
 
     // Clear any existing memory buffers
     clear();
-    free_dummy();
 
     node_procs_ = node_procs;
     node_rank_ = node_rank;
@@ -73,7 +288,16 @@ void OmpManager::assign_device(int node_procs, int node_rank, bool disabled) {
     auto & env = toast::Environment::get();
     env.set_acc(n_target, proc_per_dev, target_dev_);
 
-    allocate_dummy(n_target);
+    // Create a memory pool on the target device
+
+    // if (n_target > 0) {
+    //     double one_gb = 1024.0 * 1024.0 * 1024.0;
+    //     size_t proc_bytes = (size_t)(mem_gb * one_gb / (double)proc_per_dev);
+    //     if (pool_ != nullptr) {
+    //         delete pool_;
+    //     }
+    //     pool_ = new OmpPoolResource(target_dev_, proc_bytes);
+    // }
 
     return;
 }
@@ -126,6 +350,8 @@ void * OmpManager::create(void * buffer, size_t nbytes) {
     log.verbose(o.str().c_str());
     mem_size_[buffer] = nbytes;
     mem_[buffer] = omp_target_alloc(nbytes, target_dev_);
+
+    // mem_[buffer] = pool_->allocate(nbytes);
     if (mem_.at(buffer) == NULL) {
         o.str("");
         o << "OmpManager:  on create, host ptr " << buffer
@@ -200,6 +426,8 @@ void OmpManager::remove(void * buffer, size_t nbytes) {
         throw std::runtime_error(o.str().c_str());
     }
     mem_size_.erase(buffer);
+
+    // pool_->deallocate(mem_.at(buffer));
     omp_target_free(mem_.at(buffer), target_dev_);
     mem_.erase(buffer);
 
@@ -364,59 +592,6 @@ void OmpManager::dump() {
     return;
 }
 
-void OmpManager::allocate_dummy(int n_target) {
-    // Allocate a small dummy device buffer that can be used to represent
-    // a NULL host pointer when passing to device functions with use_device_ptr().
-
-    auto log = toast::Logger::get();
-    std::ostringstream o;
-
-    null = malloc(sizeof(int64_t));
-    dev_null_ = NULL;
-
-    if (n_target > 0) {
-        #ifdef HAVE_OPENMP_TARGET
-        dev_null_ = omp_target_alloc(sizeof(int64_t), target_dev_);
-        if (dev_null_ == NULL) {
-            o.str("");
-            o << "OmpManager:  failed to allocate dummy dev pointer";
-            log.error(o.str().c_str());
-            throw std::runtime_error(o.str().c_str());
-        }
-        int failed = omp_target_associate_ptr(
-            null, dev_null_, sizeof(int64_t), 0, target_dev_
-        );
-        if (failed != 0) {
-            o.str("");
-            o << "OmpManager:  failed to associate dev null pointer";
-            log.error(o.str().c_str());
-            throw std::runtime_error(o.str().c_str());
-        }
-        #endif // ifdef HAVE_OPENMP_TARGET
-    }
-}
-
-void OmpManager::free_dummy() {
-    auto log = toast::Logger::get();
-    std::ostringstream o;
-
-    if (dev_null_ != NULL) {
-        #ifdef HAVE_OPENMP_TARGET
-        int failed = omp_target_disassociate_ptr(null, target_dev_);
-        if (failed != 0) {
-            o.str("");
-            o << "OmpManager:  destructor failed to disassociate dev null pointer";
-            log.error(o.str().c_str());
-            throw std::runtime_error(o.str().c_str());
-        }
-        omp_target_free(dev_null_, target_dev_);
-        #endif // ifdef HAVE_OPENMP_TARGET
-    }
-    if (null != NULL) {
-        free(null);
-    }
-}
-
 OmpManager::OmpManager() : mem_size_(), mem_() {
     auto log = toast::Logger::get();
     std::ostringstream o;
@@ -427,25 +602,28 @@ OmpManager::OmpManager() : mem_size_(), mem_() {
     target_dev_ = -2;
     node_procs_ = -1;
     node_rank_ = -1;
-    null = NULL;
-    dev_null_ = NULL;
+    pool_ = nullptr;
 }
 
 OmpManager::~OmpManager() {
     auto log = toast::Logger::get();
     std::ostringstream o;
     clear();
-    free_dummy();
+    if (pool_ != nullptr) {
+        delete pool_;
+    }
 }
 
 void OmpManager::clear() {
     #ifdef HAVE_OPENMP_TARGET
+    std::vector <void *> to_clear;
     for (auto & p : mem_) {
-        omp_target_free(p.second, target_dev_);
+        to_clear.push_back(p.first);
+    }
+    for (auto & p : to_clear) {
+        remove(p, mem_size_.at(p));
     }
     #endif // ifdef HAVE_OPENMP_TARGET
-    mem_size_.clear();
-    mem_.clear();
 }
 
 void extract_accel_buffer(py::buffer_info const & info, void ** host_ptr,
@@ -457,13 +635,13 @@ void extract_accel_buffer(py::buffer_info const & info, void ** host_ptr,
         // stride " << info.strides[d] << " (itemsize = " << info.itemsize << ") raw = "
         // << (*host_ptr) << std::endl;
         (*n_elem) *= info.shape[d];
-        if (info.strides[d] != info.itemsize) {
-            auto log = toast::Logger::get();
-            std::ostringstream o;
-            o << "Cannot use python buffers with stride != itemsize.";
-            log.error(o.str().c_str());
-            throw std::runtime_error(o.str().c_str());
-        }
+    }
+    if (info.strides[info.ndim - 1] != info.itemsize) {
+        auto log = toast::Logger::get();
+        std::ostringstream o;
+        o << "Cannot use python buffers with stride of final dimension != itemsize.";
+        log.error(o.str().c_str());
+        throw std::runtime_error(o.str().c_str());
     }
     (*n_bytes) = (*n_elem) * info.itemsize;
     return;
@@ -486,10 +664,11 @@ void init_accelerator(py::module & m) {
         )");
 
     m.def(
-        "accel_assign_device", [](int node_procs, int node_rank, bool disabled)
+        "accel_assign_device",
+        [](int node_procs, int node_rank, float mem_gb, bool disabled)
         {
             auto & omgr = OmpManager::get();
-            omgr.assign_device(node_procs, node_rank, disabled);
+            omgr.assign_device(node_procs, node_rank, mem_gb, disabled);
             return;
         },
         R"(
@@ -608,6 +787,68 @@ void init_accelerator(py::module & m) {
 
     )");
 
+    m.def(
+        "accel_reset", [](py::buffer data)
+        {
+            auto & log = toast::Logger::get();
+            py::buffer_info info = data.request();
+            void * p_host;
+            size_t n_elem;
+            size_t n_bytes;
+            extract_accel_buffer(info, &p_host, &n_elem, &n_bytes);
+
+            std::ostringstream o;
+            #ifndef HAVE_OPENMP_TARGET
+            o.str("");
+            o << "TOAST not built with OpenMP target support";
+            log.error(o.str().c_str());
+            throw std::runtime_error(o.str().c_str());
+            #endif // ifndef HAVE_OPENMP_TARGET
+
+            auto & omgr = OmpManager::get();
+            int present = omgr.present(p_host, n_bytes);
+            if (present == 0) {
+                o.str("");
+                o << "Data is not present on device, cannot reset.";
+                log.error(o.str().c_str());
+                throw std::runtime_error(o.str().c_str());
+            }
+
+            o.str(
+                "");
+            o << "reset device with host pointer " << p_host << " (" << n_bytes << " bytes)";
+            log.verbose(o.str().c_str());
+
+            uint8_t * p_dev = omgr.device_ptr <uint8_t> (
+                static_cast <uint8_t *> (p_host)
+            );
+            #ifndef HAVE_OPENMP_TARGET
+            #pragma omp target teams distribute parallel for default(shared) \
+            map(to: n_bytes)                    \
+            is_device_ptr(p_dev)
+            for (size_t i = 0; i < n_bytes; ++i) {
+                p_dev[i] = 0;
+            }
+            #endif // ifndef HAVE_OPENMP_TARGET
+            o.str("");
+            o << "reset device with host pointer " << p_host << " DONE";
+            log.verbose(o.str().c_str());
+            return;
+        }, py::arg(
+            "data"),
+        R"(
+        Reset device copy of the data to zero.
+
+        This is done directly on the device, without copying zeros.
+
+        Args:
+            data (array):  The host data.
+
+        Returns:
+            None
+
+    )");
+
     m.def(
         "accel_update_device", [](py::buffer data)
         {
@@ -775,9 +1016,10 @@ void init_accelerator(py::module & m) {
             void * dev_raw = omgr.device_ptr(raw);
 
             #ifdef HAVE_OPENMP_TARGET
-            # pragma omp target data device(dev) use_device_ptr(dev_raw)
+            # pragma omp target data device(dev)
             {
-                # pragma omp target teams distribute parallel for collapse(2)
+                # pragma omp target teams distribute parallel for collapse(2) \
+                is_device_ptr(dev_raw)
                 for (int64_t i = 0; i < n_det; i++) {
                     for (int64_t j = 0; j < n_samp; j++) {
                         raw[i * n_samp + j] *= 2.0;
@@ -803,9 +1045,10 @@ void init_accelerator(py::module & m) {
             void * dev_raw = omgr.device_ptr(raw);
 
             #ifdef HAVE_OPENMP_TARGET
-            # pragma omp target data device(dev) use_device_ptr(dev_raw)
+            # pragma omp target data device(dev)
             {
-                # pragma omp target teams distribute parallel for collapse(2)
+                # pragma omp target teams distribute parallel for collapse(2) \
+                is_device_ptr(dev_raw)
                 for (size_t i = 0; i < n_det; i++) {
                     for (size_t j = 0; j < n_samp; j++) {
                         raw[i * n_samp + j] *= 2.0;
diff --git a/src/toast/_libtoast/accelerator.hpp b/src/toast/_libtoast/accelerator.hpp
index befd7a450..19601ddf0 100644
--- a/src/toast/_libtoast/accelerator.hpp
+++ b/src/toast/_libtoast/accelerator.hpp
@@ -12,11 +12,62 @@
 #endif // ifdef _OPENMP
 
 
-// Declarations for OpenMP target offload helpers
+// Declarations for OpenMP target offload helpers.
 
-// This helper class stores a mapping between host and device pointers.
-// Once the OpenMP 5.1 standard is widely implemented across compilers,
-// We can use the omp_get_mapped_ptr() method.
+typedef struct {
+    // Pointers to the beginning and end of the block
+    void * start;
+    void * end;
+
+    // Size of the usable block in bytes, not including alignment padding.
+    size_t size;
+
+    // Size including alignment
+    size_t aligned_size;
+
+    // Whether that block has been freed
+    bool is_free;
+} OmpBlock;
+
+
+class OmpPoolResource {
+    public:
+
+        OmpPoolResource(int target, size_t size, size_t align = 256);
+        OmpPoolResource() : OmpPoolResource(-1, 0) {}
+
+        OmpPoolResource(int target) : OmpPoolResource(target, 0) {}
+
+        OmpPoolResource(OmpPoolResource const &) = delete;
+
+        ~OmpPoolResource();
+
+        void release();
+
+        void * allocate(std::size_t bytes);
+        void deallocate(void * p);
+
+    private:
+
+        int target_;
+        size_t pool_size_;
+        void * raw_;
+        size_t pool_used_;
+        size_t alignment_;
+
+        // Block handling
+        std::vector <OmpBlock> blocks_;
+        std::unordered_map <void *, size_t> ptr_to_block_;
+
+        void alloc();
+        bool verbose();
+        size_t block_index(void * ptr);
+        size_t compute_aligned_bytes(size_t bytes, size_t alignment);
+        void * shift_void_pointer(void * ptr, size_t offset_bytes);
+
+        // Remove free blocks
+        void garbage_collection();
+};
 
 
 class OmpManager {
@@ -24,7 +75,8 @@ class OmpManager {
 
         static OmpManager & get();
 
-        void assign_device(int node_procs, int node_rank, bool disabled = false);
+        void assign_device(int node_procs, int node_rank, float mem_gb,
+                           bool disabled = false);
         int get_device();
         bool device_is_host();
 
@@ -39,7 +91,17 @@ class OmpManager {
 
         ~OmpManager();
 
-        void * null;
+        template <typename T>
+        T * null_ptr() {
+            static T instance = T();
+            if (!present(static_cast <void *> (&instance), sizeof(T))) {
+                // Create device copy on demand
+                void * dummy = create(
+                    static_cast <void *> (&instance), sizeof(T)
+                );
+            }
+            return &instance;
+        }
 
         template <typename T>
         T * device_ptr(T * buffer) {
@@ -75,8 +137,6 @@ class OmpManager {
 
         OmpManager();
         void clear();
-        void allocate_dummy(int n_target);
-        void free_dummy();
 
         std::unordered_map <void *, size_t> mem_size_;
         std::unordered_map <void *, void *> mem_;
@@ -84,7 +144,8 @@ class OmpManager {
         int target_dev_;
         int node_procs_;
         int node_rank_;
-        void * dev_null_;
+
+        OmpPoolResource * pool_;
 };
 
 
diff --git a/src/toast/_libtoast/module.cpp b/src/toast/_libtoast/module.cpp
index d95306c02..712c70aa3 100644
--- a/src/toast/_libtoast/module.cpp
+++ b/src/toast/_libtoast/module.cpp
@@ -49,6 +49,8 @@ PYBIND11_MODULE(_libtoast, m) {
     init_ops_stokes_weights(m);
     init_ops_pixels_healpix(m);
     init_ops_mapmaker_utils(m);
+    init_ops_noise_weight(m);
+    init_ops_scan_map(m);
     init_io_compression_flac(m);
 
     // Internal unit test runner
diff --git a/src/toast/_libtoast/module.hpp b/src/toast/_libtoast/module.hpp
index a677122d4..33e3f1d4a 100644
--- a/src/toast/_libtoast/module.hpp
+++ b/src/toast/_libtoast/module.hpp
@@ -422,6 +422,8 @@ void init_ops_pointing_detector(py::module & m);
 void init_ops_stokes_weights(py::module & m);
 void init_ops_pixels_healpix(py::module & m);
 void init_ops_mapmaker_utils(py::module & m);
+void init_ops_noise_weight(py::module & m);
+void init_ops_scan_map(py::module & m);
 void init_io_compression_flac(py::module & m);
 
 #endif // ifndef LIBTOAST_HPP
diff --git a/src/toast/_libtoast/ops_mapmaker_utils.cpp b/src/toast/_libtoast/ops_mapmaker_utils.cpp
index 41c8961f5..5a4d8b36b 100644
--- a/src/toast/_libtoast/ops_mapmaker_utils.cpp
+++ b/src/toast/_libtoast/ops_mapmaker_utils.cpp
@@ -24,8 +24,7 @@ void build_noise_weighted_inner(
     int64_t const * pixels,
     double const * weights,
     double const * det_scale,
-    double * zmap_val,
-    int64_t * zoff,
+    double * zmap,
     int64_t isamp,
     int64_t n_samp,
     int64_t idet,
@@ -34,7 +33,8 @@ void build_noise_weighted_inner(
     uint8_t shared_mask,
     int64_t n_pix_submap,
     bool use_shared_flags,
-    bool use_det_flags) {
+    bool use_det_flags
+) {
     int32_t w_indx = weight_index[idet];
     int32_t p_indx = pixel_index[idet];
     int32_t f_indx = flag_index[idet];
@@ -49,6 +49,7 @@ void build_noise_weighted_inner(
     double scaled_data;
     int64_t local_submap;
     int64_t global_submap;
+    int64_t zoff;
 
     uint8_t det_check = 0;
     if (use_det_flags) {
@@ -62,26 +63,23 @@ void build_noise_weighted_inner(
     if (
         (pixels[off_p] >= 0) &&
         (det_check == 0) &&
-        (shared_check == 0)) {
+        (shared_check == 0)
+    ) {
         // Good data, accumulate
         global_submap = (int64_t)(pixels[off_p] / n_pix_submap);
 
         local_submap = global2local[global_submap];
 
         isubpix = pixels[off_p] - global_submap * n_pix_submap;
-        (*zoff) = nnz * (local_submap * n_pix_submap + isubpix);
+        zoff = nnz * (local_submap * n_pix_submap + isubpix);
 
         off_wt = nnz * off_w;
 
         scaled_data = data[off_d] * det_scale[idet];
 
         for (int64_t iweight = 0; iweight < nnz; iweight++) {
-            zmap_val[iweight] = scaled_data * weights[off_wt + iweight];
-        }
-    } else   {
-        (*zoff) = -1;
-        for (int64_t iweight = 0; iweight < nnz; iweight++) {
-            zmap_val[iweight] = 0.0;
+            # pragma omp atomic update
+            zmap[zoff + iweight] += scaled_data * weights[off_wt + iweight];
         }
     }
     return;
@@ -109,8 +107,8 @@ void init_ops_mapmaker_utils(py::module & m) {
             py::buffer intervals,
             py::buffer shared_flags,
             uint8_t shared_flag_mask,
-            bool use_accel)
-        {
+            bool use_accel
+        ) {
             auto & omgr = OmpManager::get();
             int dev = omgr.get_device();
             bool offload = (!omgr.device_is_host()) && use_accel;
@@ -132,10 +130,21 @@ void init_ops_mapmaker_utils(py::module & m) {
                 weight_index, "weight_index", 1, temp_shape, {n_det}
             );
 
-            double * raw_weights = extract_buffer <double> (
-                weights, "weights", 3, temp_shape, {-1, n_samp, -1}
-            );
-            int64_t nnz = temp_shape[2];
+            // Handle the case of either 2 or 3 dims
+            auto winfo = weights.request();
+            double * raw_weights;
+            int64_t nnz;
+            if (winfo.ndim == 2) {
+                nnz = 1;
+                raw_weights = extract_buffer <double> (
+                    weights, "weights", 2, temp_shape, {-1, n_samp}
+                );
+            } else {
+                raw_weights = extract_buffer <double> (
+                    weights, "weights", 3, temp_shape, {-1, n_samp, -1}
+                );
+                nnz = temp_shape[2];
+            }
 
             int32_t * raw_data_index = extract_buffer <int32_t> (
                 data_index, "data_index", 1, temp_shape, {n_det}
@@ -174,7 +183,7 @@ void init_ops_mapmaker_utils(py::module & m) {
                 shared_flags, "flags", 1, temp_shape, {-1}
             );
             if (temp_shape[0] != n_samp) {
-                raw_shared_flags = (uint8_t *)omgr.null;
+                raw_shared_flags = omgr.null_ptr <uint8_t> ();
                 use_shared_flags = false;
             }
 
@@ -183,7 +192,7 @@ void init_ops_mapmaker_utils(py::module & m) {
                 det_flags, "det_flags", 2, temp_shape, {-1, -1}
             );
             if (temp_shape[1] != n_samp) {
-                raw_det_flags = (uint8_t *)omgr.null;
+                raw_det_flags = omgr.null_ptr <uint8_t> ();
                 use_det_flags = false;
             }
 
@@ -198,32 +207,25 @@ void init_ops_mapmaker_utils(py::module & m) {
                 uint8_t * dev_shared_flags = omgr.device_ptr(raw_shared_flags);
                 uint8_t * dev_det_flags = omgr.device_ptr(raw_det_flags);
 
-                # pragma omp target data              \
-                device(dev)                           \
-                map(to                                \
-                :                                     \
-                raw_weight_index [0:n_det],           \
-                raw_pixel_index [0:n_det],            \
-                raw_flag_index [0:n_det],             \
-                raw_data_index [0:n_det],             \
-                raw_det_scale [0:n_det],              \
-                raw_global2local [0:n_global_submap], \
-                n_view,                               \
-                n_det,                                \
-                n_samp,                               \
-                nnz,                                  \
-                n_pix_submap,                         \
-                det_flag_mask,                        \
-                shared_flag_mask,                     \
-                use_shared_flags,                     \
-                use_det_flags)                        \
-                use_device_ptr(                       \
-                raw_weight_index,                     \
-                raw_pixel_index,                      \
-                raw_flag_index,                       \
-                raw_data_index,                       \
-                raw_det_scale,                        \
-                raw_global2local)
+                # pragma omp target data             \
+                device(dev)                          \
+                map(to:                              \
+                raw_weight_index[0:n_det],           \
+                raw_pixel_index[0:n_det],            \
+                raw_flag_index[0:n_det],             \
+                raw_data_index[0:n_det],             \
+                raw_det_scale[0:n_det],              \
+                raw_global2local[0:n_global_submap], \
+                n_view,                              \
+                n_det,                               \
+                n_samp,                              \
+                nnz,                                 \
+                n_pix_submap,                        \
+                det_flag_mask,                       \
+                shared_flag_mask,                    \
+                use_shared_flags,                    \
+                use_det_flags                        \
+                )
                 {
                     # pragma omp target teams distribute collapse(2) \
                     is_device_ptr(                                   \
@@ -233,44 +235,41 @@ void init_ops_mapmaker_utils(py::module & m) {
                     dev_det_flags,                                   \
                     dev_intervals,                                   \
                     dev_shared_flags,                                \
-                    dev_zmap)
+                    dev_zmap                                         \
+                    )
                     for (int64_t idet = 0; idet < n_det; idet++) {
                         for (int64_t iview = 0; iview < n_view; iview++) {
-                            # pragma omp parallel for default(shared)
-                            for (
-                                int64_t isamp = dev_intervals[iview].first;
-                                isamp <= dev_intervals[iview].last;
-                                isamp++
-                            ) {
-                                double zmap_val[nnz];
-                                int64_t zoff;
-                                build_noise_weighted_inner(
-                                    raw_pixel_index,
-                                    raw_weight_index,
-                                    raw_flag_index,
-                                    raw_data_index,
-                                    raw_global2local,
-                                    dev_det_data,
-                                    dev_det_flags,
-                                    dev_shared_flags,
-                                    dev_pixels,
-                                    dev_weights,
-                                    raw_det_scale,
-                                    zmap_val,
-                                    &zoff,
-                                    isamp,
-                                    n_samp,
-                                    idet,
-                                    nnz,
-                                    det_flag_mask,
-                                    shared_flag_mask,
-                                    n_pix_submap,
-                                    use_shared_flags,
-                                    use_det_flags
-                                );
-                                for (int64_t iw = 0; iw < nnz; iw++) {
-                                    # pragma omp atomic
-                                    dev_zmap[zoff + iw] += zmap_val[iw];
+                            # pragma omp parallel default(shared)
+                            {
+                                # pragma omp for
+                                for (
+                                    int64_t isamp = dev_intervals[iview].first;
+                                    isamp <= dev_intervals[iview].last;
+                                    isamp++
+                                ) {
+                                    build_noise_weighted_inner(
+                                        raw_pixel_index,
+                                        raw_weight_index,
+                                        raw_flag_index,
+                                        raw_data_index,
+                                        raw_global2local,
+                                        dev_det_data,
+                                        dev_det_flags,
+                                        dev_shared_flags,
+                                        dev_pixels,
+                                        dev_weights,
+                                        raw_det_scale,
+                                        dev_zmap,
+                                        isamp,
+                                        n_samp,
+                                        idet,
+                                        nnz,
+                                        det_flag_mask,
+                                        shared_flag_mask,
+                                        n_pix_submap,
+                                        use_shared_flags,
+                                        use_det_flags
+                                    );
                                 }
                             }
                         }
@@ -283,9 +282,6 @@ void init_ops_mapmaker_utils(py::module & m) {
                     for (int64_t iview = 0; iview < n_view; iview++) {
                         #pragma omp parallel default(shared)
                         {
-                            double zmap_val[nnz];
-                            int64_t zoff;
-                            
                             #pragma omp for
                             for (
                                 int64_t isamp = raw_intervals[iview].first;
@@ -304,8 +300,7 @@ void init_ops_mapmaker_utils(py::module & m) {
                                     raw_pixels,
                                     raw_weights,
                                     raw_det_scale,
-                                    zmap_val,
-                                    &zoff,
+                                    raw_zmap,
                                     isamp,
                                     n_samp,
                                     idet,
@@ -316,11 +311,6 @@ void init_ops_mapmaker_utils(py::module & m) {
                                     use_shared_flags,
                                     use_det_flags
                                 );
-
-                                for (int64_t iw = 0; iw < nnz; iw++) {
-                                    #pragma omp atomic
-                                    raw_zmap[zoff + iw] += zmap_val[iw];
-                                }
                             }
                         }
                     }
diff --git a/src/toast/_libtoast/ops_noise_weight.cpp b/src/toast/_libtoast/ops_noise_weight.cpp
new file mode 100644
index 000000000..9fb621d02
--- /dev/null
+++ b/src/toast/_libtoast/ops_noise_weight.cpp
@@ -0,0 +1,107 @@
+// Copyright (c) 2023-2023 by the parties listed in the AUTHORS file.
+// All rights reserved.  Use of this source code is governed by
+// a BSD-style license that can be found in the LICENSE file.
+
+#include <module.hpp>
+
+#include <intervals.hpp>
+
+#include <accelerator.hpp>
+
+
+void init_ops_noise_weight(py::module & m) {
+    m.def(
+        "noise_weight", [](
+            py::buffer det_data,
+            py::buffer data_index,
+            py::buffer intervals,
+            py::buffer detector_weights,
+            bool use_accel
+        ) {
+            auto & omgr = OmpManager::get();
+            int dev = omgr.get_device();
+            bool offload = (!omgr.device_is_host()) && use_accel;
+
+            // This is used to return the actual shape of each buffer
+            std::vector <int64_t> temp_shape(3);
+
+            int32_t * raw_data_index = extract_buffer <int32_t> (
+                data_index, "data_index", 1, temp_shape, {-1}
+            );
+            int64_t n_det = temp_shape[0];
+
+            double * raw_det_data = extract_buffer <double> (
+                det_data, "det_data", 2, temp_shape, {-1, -1}
+            );
+            int64_t n_samp = temp_shape[1];
+
+            Interval * raw_intervals = extract_buffer <Interval> (
+                intervals, "intervals", 1, temp_shape, {-1}
+            );
+            int64_t n_view = temp_shape[0];
+
+            double * raw_det_weights = extract_buffer <double> (
+                detector_weights, "detector_weights", 1, temp_shape, {n_det}
+            );
+
+            if (offload) {
+                #ifdef HAVE_OPENMP_TARGET
+
+                double * dev_det_data = omgr.device_ptr(raw_det_data);
+                Interval * dev_intervals = omgr.device_ptr(raw_intervals);
+
+                # pragma omp target data  \
+                device(dev)               \
+                map(to:                   \
+                raw_data_index[0:n_det],  \
+                raw_det_weights[0:n_det], \
+                n_view,                   \
+                n_det,                    \
+                n_samp                    \
+                )
+                {
+                    # pragma omp target teams distribute collapse(2) \
+                    is_device_ptr(                                   \
+                    dev_det_data,                                    \
+                    dev_intervals                                    \
+                    )
+                    for (int64_t idet = 0; idet < n_det; idet++) {
+                        for (int64_t iview = 0; iview < n_view; iview++) {
+                            # pragma omp parallel
+                            {
+                                # pragma omp for default(shared)
+                                for (
+                                    int64_t isamp = dev_intervals[iview].first;
+                                    isamp <= dev_intervals[iview].last;
+                                    isamp++
+                                ) {
+                                    int32_t d_indx = raw_data_index[idet];
+                                    int64_t off_d = d_indx * n_samp + isamp;
+                                    dev_det_data[off_d] *= raw_det_weights[idet];
+                                }
+                            }
+                        }
+                    }
+                }
+
+                #endif // ifdef HAVE_OPENMP_TARGET
+            } else {
+                for (int64_t idet = 0; idet < n_det; idet++) {
+                    for (int64_t iview = 0; iview < n_view; iview++) {
+                        #pragma omp parallel for default(shared)
+                        for (
+                            int64_t isamp = raw_intervals[iview].first;
+                            isamp <= raw_intervals[iview].last;
+                            isamp++
+                        ) {
+                            int32_t d_indx = raw_data_index[idet];
+                            int64_t off_d = d_indx * n_samp + isamp;
+                            raw_det_data[off_d] *= raw_det_weights[idet];
+                        }
+                    }
+                }
+            }
+
+            return;
+        });
+}
diff --git a/src/toast/_libtoast/ops_pixels_healpix.cpp b/src/toast/_libtoast/ops_pixels_healpix.cpp
index ca117d892..453f2fcc5 100644
--- a/src/toast/_libtoast/ops_pixels_healpix.cpp
+++ b/src/toast/_libtoast/ops_pixels_healpix.cpp
@@ -17,15 +17,20 @@
 // 2/3
 #define TWOTHIRDS 0.66666666666666666667
 
+// Helper table initialization
+void hpix_init_utab(uint64_t * utab) {
+    for (uint64_t m = 0; m < 0x100; ++m) {
+        utab[m] = (m & 0x1) | ((m & 0x2) << 1) | ((m & 0x4) << 2) |
+                  ((m & 0x8) << 3) | ((m & 0x10) << 4) | ((m & 0x20) << 5) |
+                  ((m & 0x40) << 6) | ((m & 0x80) << 7);
+    }
+    return;
+}
 
 #ifdef HAVE_OPENMP_TARGET
 # pragma omp declare target
 #endif // ifdef HAVE_OPENMP_TARGET
 
-// FIXME:  this ridiculous code duplication is due to nvc++
-// not supporting loadable device objects in shared libraries.
-// So we must duplicate this across compilation units.
-
 void pixels_healpix_qa_rotate(double const * q_in, double const * v_in,
                               double * v_out) {
     // The input quaternion has already been normalized on the host.
@@ -52,69 +57,16 @@ void pixels_healpix_qa_rotate(double const * q_in, double const * v_in,
     return;
 }
 
-typedef struct {
-    int64_t nside;
-    int64_t npix;
-    int64_t ncap;
-    double dnside;
-    int64_t twonside;
-    int64_t fournside;
-    int64_t nsideplusone;
-    int64_t nsideminusone;
-    double halfnside;
-    double tqnside;
-    int64_t factor;
-    int64_t jr[12];
-    int64_t jp[12];
-    uint64_t utab[0x100];
-    uint64_t ctab[0x100];
-} hpix;
-
-void hpix_init(hpix * hp, int64_t nside) {
-    hp->nside = nside;
-    hp->ncap = 2 * (nside * nside - nside);
-    hp->npix = 12 * nside * nside;
-    hp->dnside = static_cast <double> (nside);
-    hp->twonside = 2 * nside;
-    hp->fournside = 4 * nside;
-    hp->nsideplusone = nside + 1;
-    hp->halfnside = 0.5 * (hp->dnside);
-    hp->tqnside = 0.75 * (hp->dnside);
-    hp->factor = 0;
-    hp->nsideminusone = nside - 1;
-    while (nside != (1ll << hp->factor)) {
-        ++hp->factor;
-    }
-
-    int64_t init_jr[12] = {2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4};
-    memcpy(hp->jr, init_jr, sizeof(init_jr));
-
-    int64_t init_jp[12] = {1, 3, 5, 7, 0, 2, 4, 6, 1, 3, 5, 7};
-    memcpy(hp->jp, init_jp, sizeof(init_jp));
-
-    for (uint64_t m = 0; m < 0x100; ++m) {
-        hp->utab[m] = (m & 0x1) | ((m & 0x2) << 1) | ((m & 0x4) << 2) |
-                      ((m & 0x8) << 3) | ((m & 0x10) << 4) | ((m & 0x20) << 5) |
-                      ((m & 0x40) << 6) | ((m & 0x80) << 7);
-
-        hp->ctab[m] = (m & 0x1) | ((m & 0x2) << 7) | ((m & 0x4) >> 1) |
-                      ((m & 0x8) << 6) | ((m & 0x10) >> 2) | ((m & 0x20) << 5) |
-                      ((m & 0x40) >> 3) | ((m & 0x80) << 4);
-    }
-    return;
-}
-
-uint64_t hpix_xy2pix(hpix * hp, uint64_t x, uint64_t y) {
-    return hp->utab[x & 0xff] | (hp->utab[(x >> 8) & 0xff] << 16) |
-           (hp->utab[(x >> 16) & 0xff] << 32) |
-           (hp->utab[(x >> 24) & 0xff] << 48) |
-           (hp->utab[y & 0xff] << 1) | (hp->utab[(y >> 8) & 0xff] << 17) |
-           (hp->utab[(y >> 16) & 0xff] << 33) |
-           (hp->utab[(y >> 24) & 0xff] << 49);
+uint64_t hpix_xy2pix(uint64_t * utab, uint64_t x, uint64_t y) {
+    return utab[x & 0xff] | (utab[(x >> 8) & 0xff] << 16) |
+           (utab[(x >> 16) & 0xff] << 32) |
+           (utab[(x >> 24) & 0xff] << 48) |
+           (utab[y & 0xff] << 1) | (utab[(y >> 8) & 0xff] << 17) |
+           (utab[(y >> 16) & 0xff] << 33) |
+           (utab[(y >> 24) & 0xff] << 49);
 }
 
-void hpix_vec2zphi(hpix * hp, double const * vec,
-                   double * phi, int * region, double * z,
+void hpix_vec2zphi(double const * vec, double * phi, int * region, double * z,
                    double * rtz) {
     // region encodes BOTH the sign of Z and whether its
     // absolute value is greater than 2/3.
@@ -127,8 +79,8 @@ void hpix_vec2zphi(hpix * hp, double const * vec,
     return;
 }
 
-void hpix_zphi2nest(hpix * hp, double phi, int region, double z,
-                    double rtz, int64_t * pix) {
+void hpix_zphi2nest(int64_t nside, int64_t factor, uint64_t * utab, double phi,
+                    int region, double z, double rtz, int64_t * pix) {
     double tt = (phi >= 0.0) ? phi * TWOINVPI : phi * TWOINVPI + 4.0;
     int64_t x;
     int64_t y;
@@ -142,15 +94,21 @@ void hpix_zphi2nest(hpix * hp, double phi, int region, double z,
     int64_t ntt;
     double tp;
 
+    double dnside = static_cast <double> (nside);
+    int64_t twonside = 2 * nside;
+    double halfnside = 0.5 * dnside;
+    double tqnside = 0.75 * dnside;
+    int64_t nsideminusone = nside - 1;
+
     if ((region == 1) || (region == -1)) {
-        temp1 = hp->halfnside + hp->dnside * tt;
-        temp2 = hp->tqnside * z;
+        temp1 = halfnside + dnside * tt;
+        temp2 = tqnside * z;
 
         jp = (int64_t)(temp1 - temp2);
         jm = (int64_t)(temp1 + temp2);
 
-        ifp = jp >> hp->factor;
-        ifm = jm >> hp->factor;
+        ifp = jp >> factor;
+        ifm = jm >> factor;
 
         if (ifp == ifm) {
             face = (ifp == 4) ? (int64_t)4 : ifp + 4;
@@ -160,29 +118,29 @@ void hpix_zphi2nest(hpix * hp, double phi, int region, double z,
             face = ifm + 8;
         }
 
-        x = jm & hp->nsideminusone;
-        y = hp->nsideminusone - (jp & hp->nsideminusone);
+        x = jm & nsideminusone;
+        y = nsideminusone - (jp & nsideminusone);
     } else {
         ntt = (int64_t)tt;
 
         tp = tt - (double)ntt;
 
-        temp1 = hp->dnside * rtz;
+        temp1 = dnside * rtz;
 
         jp = (int64_t)(tp * temp1);
         jm = (int64_t)((1.0 - tp) * temp1);
 
-        if (jp >= hp->nside) {
-            jp = hp->nsideminusone;
+        if (jp >= nside) {
+            jp = nsideminusone;
         }
-        if (jm >= hp->nside) {
-            jm = hp->nsideminusone;
+        if (jm >= nside) {
+            jm = nsideminusone;
         }
 
         if (z >= 0) {
             face = ntt;
-            x = hp->nsideminusone - jm;
-            y = hp->nsideminusone - jp;
+            x = nsideminusone - jm;
+            y = nsideminusone - jp;
         } else {
             face = ntt + 8;
             x = jp;
@@ -190,14 +148,14 @@ void hpix_zphi2nest(hpix * hp, double phi, int region, double z,
         }
     }
 
-    uint64_t sipf = hpix_xy2pix(hp, (uint64_t)x, (uint64_t)y);
+    uint64_t sipf = hpix_xy2pix(utab, (uint64_t)x, (uint64_t)y);
 
-    (*pix) = (int64_t)sipf + (face << (2 * hp->factor));
+    (*pix) = (int64_t)sipf + (face << (2 * factor));
 
     return;
 }
 
-void hpix_zphi2ring(hpix * hp, double phi, int region, double z,
+void hpix_zphi2ring(int64_t nside, int64_t factor, double phi, int region, double z,
                     double rtz, int64_t * pix) {
     double tt = (phi >= 0.0) ? phi * TWOINVPI : phi * TWOINVPI + 4.0;
     double tp;
@@ -210,24 +168,32 @@ void hpix_zphi2ring(hpix * hp, double phi, int region, double z,
     int64_t ir;
     int64_t kshift;
 
+    double dnside = static_cast <double> (nside);
+    int64_t fournside = 4 * nside;
+    double halfnside = 0.5 * dnside;
+    double tqnside = 0.75 * dnside;
+    int64_t nsideplusone = nside + 1;
+    int64_t ncap = 2 * (nside * nside - nside);
+    int64_t npix = 12 * nside * nside;
+
     if ((region == 1) || (region == -1)) {
-        temp1 = hp->halfnside + hp->dnside * tt;
-        temp2 = hp->tqnside * z;
+        temp1 = halfnside + dnside * tt;
+        temp2 = tqnside * z;
 
         jp = (int64_t)(temp1 - temp2);
         jm = (int64_t)(temp1 + temp2);
 
-        ir = hp->nsideplusone + jp - jm;
+        ir = nsideplusone + jp - jm;
         kshift = 1 - (ir & 1);
 
-        ip = (jp + jm - hp->nside + kshift + 1) >> 1;
-        ip = ip % hp->fournside;
+        ip = (jp + jm - nside + kshift + 1) >> 1;
+        ip = ip % fournside;
 
-        (*pix) = hp->ncap + ((ir - 1) * hp->fournside + ip);
+        (*pix) = ncap + ((ir - 1) * fournside + ip);
     } else {
         tp = tt - floor(tt);
 
-        temp1 = hp->dnside * rtz;
+        temp1 = dnside * rtz;
 
         jp = (int64_t)(tp * temp1);
         jm = (int64_t)((1.0 - tp) * temp1);
@@ -237,14 +203,16 @@ void hpix_zphi2ring(hpix * hp, double phi, int region, double z,
         ip -= longpart;
 
         (*pix) = (region > 0) ? (2 * ir * (ir - 1) + ip)
-                 : (hp->npix - 2 * ir * (ir + 1) + ip);
+                 : (npix - 2 * ir * (ir + 1) + ip);
     }
 
     return;
 }
 
 void pixels_healpix_nest_inner(
-    hpix & hp,
+    int64_t nside,
+    int64_t factor,
+    uint64_t * utab,
     int32_t const * quat_index,
     int32_t const * pixel_index,
     double const * quats,
@@ -271,8 +239,8 @@ void pixels_healpix_nest_inner(
     int64_t sub_map;
 
     pixels_healpix_qa_rotate(&(quats[qoff]), zaxis, dir);
-    hpix_vec2zphi(&hp, dir, &phi, &region, &z, &rtz);
-    hpix_zphi2nest(&hp, phi, region, z, rtz, &(pixels[poff]));
+    hpix_vec2zphi(dir, &phi, &region, &z, &rtz);
+    hpix_zphi2nest(nside, factor, utab, phi, region, z, rtz, &(pixels[poff]));
     if (use_flags && ((flags[isamp] & mask) != 0)) {
         pixels[poff] = -1;
     } else {
@@ -284,7 +252,8 @@ void pixels_healpix_nest_inner(
 }
 
 void pixels_healpix_ring_inner(
-    hpix & hp,
+    int64_t nside,
+    int64_t factor,
     int32_t const * quat_index,
     int32_t const * pixel_index,
     double const * quats,
@@ -311,8 +280,8 @@ void pixels_healpix_ring_inner(
     int64_t sub_map;
 
     pixels_healpix_qa_rotate(&(quats[qoff]), zaxis, dir);
-    hpix_vec2zphi(&hp, dir, &phi, &region, &z, &rtz);
-    hpix_zphi2ring(&hp, phi, region, z, rtz, &(pixels[poff]));
+    hpix_vec2zphi(dir, &phi, &region, &z, &rtz);
+    hpix_zphi2ring(nside, factor, phi, region, z, rtz, &(pixels[poff]));
     if (use_flags && ((flags[isamp] & mask) != 0)) {
         pixels[poff] = -1;
     } else {
@@ -376,6 +345,7 @@ void init_ops_pixels_healpix(py::module & m) {
             uint8_t * raw_hsub = extract_buffer <uint8_t> (
                 hit_submaps, "hit_submaps", 1, temp_shape, {-1}
             );
+            int64_t n_submap = temp_shape[0];
 
             // Optionally use flags
             bool use_flags = true;
@@ -383,10 +353,18 @@ void init_ops_pixels_healpix(py::module & m) {
                 shared_flags, "flags", 1, temp_shape, {-1}
             );
             if (temp_shape[0] != n_samp) {
-                raw_flags = (uint8_t *)omgr.null;
+                raw_flags = omgr.null_ptr <uint8_t> ();
                 use_flags = false;
             }
 
+            static uint64_t utab[0x100];
+            hpix_init_utab(utab);
+
+            int64_t factor = 0;
+            while (nside != (1ll << factor)) {
+                ++factor;
+            }
+
             if (offload) {
                 #ifdef HAVE_OPENMP_TARGET
 
@@ -395,6 +373,16 @@ void init_ops_pixels_healpix(py::module & m) {
                 Interval * dev_intervals = omgr.device_ptr(raw_intervals);
                 uint8_t * dev_flags = omgr.device_ptr(raw_flags);
 
+                // Make sure the lookup table exists on device
+                size_t utab_bytes = 0x100 * sizeof(uint64_t);
+                int present_utab = omgr.present(utab, utab_bytes);
+                if (present_utab == 0) {
+                    void * vptr = omgr.create((void *)utab, utab_bytes);
+                    omgr.update_device((void *)utab, utab_bytes);
+                }
+
+                uint64_t * dev_utab = omgr.device_ptr(utab);
+
                 # pragma omp target data  \
                 device(dev)               \
                 map(to:                   \
@@ -402,6 +390,7 @@ void init_ops_pixels_healpix(py::module & m) {
                 raw_quat_index[0:n_det],  \
                 n_pix_submap,             \
                 nside,                    \
+                factor,                   \
                 nest,                     \
                 n_view,                   \
                 n_det,                    \
@@ -409,46 +398,45 @@ void init_ops_pixels_healpix(py::module & m) {
                 shared_flag_mask,         \
                 use_flags                 \
                 )                         \
-                map(tofrom: raw_hsub)     \
-                use_device_ptr(           \
-                raw_pixel_index,          \
-                raw_quat_index,           \
-                raw_hsub                  \
-                )
+                map(tofrom: raw_hsub[0:n_submap])
                 {
-                    hpix hp;
-                    hpix_init(&hp, nside);
                     if (nest) {
                         # pragma omp target teams distribute collapse(2) \
                         is_device_ptr(                                   \
                         dev_pixels,                                      \
                         dev_quats,                                       \
                         dev_flags,                                       \
-                        dev_intervals                                    \
+                        dev_intervals,                                   \
+                        dev_utab                                         \
                         )
                         for (int64_t idet = 0; idet < n_det; idet++) {
                             for (int64_t iview = 0; iview < n_view; iview++) {
-                                # pragma omp parallel for default(shared)
-                                for (
-                                    int64_t isamp = dev_intervals[iview].first;
-                                    isamp <= dev_intervals[iview].last;
-                                    isamp++
-                                ) {
-                                    pixels_healpix_nest_inner(
-                                        hp,
-                                        raw_quat_index,
-                                        raw_pixel_index,
-                                        dev_quats,
-                                        dev_flags,
-                                        raw_hsub,
-                                        dev_pixels,
-                                        n_pix_submap,
-                                        isamp,
-                                        n_samp,
-                                        idet,
-                                        shared_flag_mask,
-                                        use_flags
-                                    );
+                                # pragma omp parallel
+                                {
+                                    # pragma omp for default(shared)
+                                    for (
+                                        int64_t isamp = dev_intervals[iview].first;
+                                        isamp <= dev_intervals[iview].last;
+                                        isamp++
+                                    ) {
+                                        pixels_healpix_nest_inner(
+                                            nside,
+                                            factor,
+                                            dev_utab,
+                                            raw_quat_index,
+                                            raw_pixel_index,
+                                            dev_quats,
+                                            dev_flags,
+                                            raw_hsub,
+                                            dev_pixels,
+                                            n_pix_submap,
+                                            isamp,
+                                            n_samp,
+                                            idet,
+                                            shared_flag_mask,
+                                            use_flags
+                                        );
+                                    }
                                 }
                             }
                         }
@@ -458,31 +446,36 @@ void init_ops_pixels_healpix(py::module & m) {
                         dev_pixels,                                      \
                         dev_quats,                                       \
                         dev_flags,                                       \
-                        dev_intervals                                    \
+                        dev_intervals,                                   \
+                        dev_utab                                         \
                         )
                         for (int64_t idet = 0; idet < n_det; idet++) {
                             for (int64_t iview = 0; iview < n_view; iview++) {
-                                # pragma omp parallel for default(shared)
-                                for (
-                                    int64_t isamp = dev_intervals[iview].first;
-                                    isamp <= dev_intervals[iview].last;
-                                    isamp++
-                                ) {
-                                    pixels_healpix_ring_inner(
-                                        hp,
-                                        raw_quat_index,
-                                        raw_pixel_index,
-                                        dev_quats,
-                                        dev_flags,
-                                        raw_hsub,
-                                        dev_pixels,
-                                        n_pix_submap,
-                                        isamp,
-                                        n_samp,
-                                        idet,
-                                        shared_flag_mask,
-                                        use_flags
-                                    );
+                                # pragma omp parallel
+                                {
+                                    # pragma omp for default(shared)
+                                    for (
+                                        int64_t isamp = dev_intervals[iview].first;
+                                        isamp <= dev_intervals[iview].last;
+                                        isamp++
+                                    ) {
+                                        pixels_healpix_ring_inner(
+                                            nside,
+                                            factor,
+                                            raw_quat_index,
+                                            raw_pixel_index,
+                                            dev_quats,
+                                            dev_flags,
+                                            raw_hsub,
+                                            dev_pixels,
+                                            n_pix_submap,
+                                            isamp,
+                                            n_samp,
+                                            idet,
+                                            shared_flag_mask,
+                                            use_flags
+                                        );
+                                    }
                                 }
                             }
                         }
@@ -491,8 +484,6 @@ void init_ops_pixels_healpix(py::module & m) {
 
                 #endif // ifdef HAVE_OPENMP_TARGET
             } else {
-                hpix hp;
-                hpix_init(&hp, nside);
                 if (nest) {
                     for (int64_t idet = 0; idet < n_det; idet++) {
                         for (int64_t iview = 0; iview < n_view; iview++) {
@@ -503,7 +494,9 @@ void init_ops_pixels_healpix(py::module & m) {
                                 isamp++
                             ) {
                                 pixels_healpix_nest_inner(
-                                    hp,
+                                    nside,
+                                    factor,
+                                    utab,
                                     raw_quat_index,
                                     raw_pixel_index,
                                     raw_quats,
@@ -530,7 +523,8 @@ void init_ops_pixels_healpix(py::module & m) {
                                 isamp++
                             ) {
                                 pixels_healpix_ring_inner(
-                                    hp,
+                                    nside,
+                                    factor,
                                     raw_quat_index,
                                     raw_pixel_index,
                                     raw_quats,
diff --git a/src/toast/_libtoast/ops_pointing_detector.cpp b/src/toast/_libtoast/ops_pointing_detector.cpp
index 9f5948688..c46948146 100644
--- a/src/toast/_libtoast/ops_pointing_detector.cpp
+++ b/src/toast/_libtoast/ops_pointing_detector.cpp
@@ -124,7 +124,7 @@ void init_ops_pointing_detector(py::module & m) {
                 shared_flags, "flags", 1, temp_shape, {-1}
             );
             if (temp_shape[0] != n_samp) {
-                raw_flags = (uint8_t *)omgr.null;
+                raw_flags = omgr.null_ptr <uint8_t> ();
                 use_flags = false;
             }
 
@@ -146,10 +146,6 @@ void init_ops_pointing_detector(py::module & m) {
                 n_det,                     \
                 n_samp,                    \
                 use_flags                  \
-                )                          \
-                use_device_ptr(            \
-                raw_focalplane,            \
-                raw_quat_index             \
                 )
                 {
                     # pragma omp target teams distribute collapse(2) \
@@ -161,24 +157,27 @@ void init_ops_pointing_detector(py::module & m) {
                     )
                     for (int64_t idet = 0; idet < n_det; idet++) {
                         for (int64_t iview = 0; iview < n_view; iview++) {
-                            # pragma omp parallel for default(shared)
-                            for (
-                                int64_t isamp = dev_intervals[iview].first;
-                                isamp <= dev_intervals[iview].last;
-                                isamp++
-                            ) {
-                                pointing_detector_inner(
-                                    raw_quat_index,
-                                    dev_flags,
-                                    dev_boresight,
-                                    raw_focalplane,
-                                    dev_quats,
-                                    isamp,
-                                    n_samp,
-                                    idet,
-                                    shared_flag_mask,
-                                    use_flags
-                                );
+                            # pragma omp parallel
+                            {
+                                # pragma omp for default(shared)
+                                for (
+                                    int64_t isamp = dev_intervals[iview].first;
+                                    isamp <= dev_intervals[iview].last;
+                                    isamp++
+                                ) {
+                                    pointing_detector_inner(
+                                        raw_quat_index,
+                                        dev_flags,
+                                        dev_boresight,
+                                        raw_focalplane,
+                                        dev_quats,
+                                        isamp,
+                                        n_samp,
+                                        idet,
+                                        shared_flag_mask,
+                                        use_flags
+                                    );
+                                }
                             }
                         }
                     }
diff --git a/src/toast/_libtoast/ops_scan_map.cpp b/src/toast/_libtoast/ops_scan_map.cpp
new file mode 100644
index 000000000..0c9bc962b
--- /dev/null
+++ b/src/toast/_libtoast/ops_scan_map.cpp
@@ -0,0 +1,278 @@
+// Copyright (c) 2023-2023 by the parties listed in the AUTHORS file.
+// All rights reserved.  Use of this source code is governed by
+// a BSD-style license that can be found in the LICENSE file.
+
+#include <module.hpp>
+
+#include <accelerator.hpp>
+
+#include <intervals.hpp>
+
+#ifdef HAVE_OPENMP_TARGET
+# pragma omp declare target
+#endif // ifdef HAVE_OPENMP_TARGET
+
+template <typename T>
+void scan_map_inner(
+    int32_t const * pixel_index,
+    int32_t const * weight_index,
+    int32_t const * data_index,
+    int64_t const * global2local,
+    double * data,
+    int64_t const * pixels,
+    double const * weights,
+    T const * mapdata,
+    double data_scale,
+    int64_t isamp,
+    int64_t n_samp,
+    int64_t idet,
+    int64_t nnz,
+    int64_t n_pix_submap,
+    bool should_zero,
+    bool should_subtract,
+    bool should_scale
+) {
+    int32_t w_indx = weight_index[idet];
+    int32_t p_indx = pixel_index[idet];
+    int32_t d_indx = data_index[idet];
+
+    int64_t off_p = p_indx * n_samp + isamp;
+    int64_t off_w = w_indx * n_samp + isamp;
+    int64_t off_d = d_indx * n_samp + isamp;
+    int64_t isubpix;
+    int64_t off_wt;
+    double tod_val;
+    int64_t local_submap;
+    int64_t global_submap;
+    int64_t map_off;
+
+    if (should_zero) {
+        data[off_d] = 0.0;
+    }
+
+    if (pixels[off_p] >= 0) {
+        // Good data, accumulate
+        global_submap = (int64_t)(pixels[off_p] / n_pix_submap);
+
+        local_submap = global2local[global_submap];
+
+        isubpix = pixels[off_p] - global_submap * n_pix_submap;
+        map_off = nnz * (local_submap * n_pix_submap + isubpix);
+
+        off_wt = nnz * off_w;
+
+        tod_val = 0.0;
+        for (int64_t iweight = 0; iweight < nnz; iweight++) {
+            tod_val += weights[off_wt + iweight] * mapdata[map_off + iweight];
+        }
+        tod_val *= data_scale;
+
+        if (should_subtract) {
+            data[off_d] -= tod_val;
+        } else if (should_scale) {
+            data[off_d] *= tod_val;
+        } else {
+            data[off_d] += tod_val;
+        }
+    }
+    return;
+}
+
+#ifdef HAVE_OPENMP_TARGET
+# pragma omp end declare target
+#endif // ifdef HAVE_OPENMP_TARGET
+
+template <typename T>
+void register_ops_scan_map(py::module & m, char const * name) {
+    m.def(name,
+          [](
+              py::buffer global2local,
+              int64_t n_pix_submap,
+              py::buffer mapdata,
+              py::buffer det_data,
+              py::buffer data_index,
+              py::buffer pixels,
+              py::buffer pixel_index,
+              py::buffer weights,
+              py::buffer weight_index,
+              py::buffer intervals,
+              double data_scale,
+              bool should_zero,
+              bool should_subtract,
+              bool should_scale,
+              bool use_accel
+          ) {
+              auto & omgr = OmpManager::get();
+              int dev = omgr.get_device();
+              bool offload = (!omgr.device_is_host()) && use_accel;
+
+              // This is used to return the actual shape of each buffer
+              std::vector <int64_t> temp_shape(3);
+
+              int32_t * raw_pixel_index = extract_buffer <int32_t> (
+                  pixel_index, "pixel_index", 1, temp_shape, {-1}
+              );
+              int64_t n_det = temp_shape[0];
+
+              int64_t * raw_pixels = extract_buffer <int64_t> (
+                  pixels, "pixels", 2, temp_shape, {-1, -1}
+              );
+              int64_t n_samp = temp_shape[1];
+
+              int32_t * raw_weight_index = extract_buffer <int32_t> (
+                  weight_index, "weight_index", 1, temp_shape, {n_det}
+              );
+
+              // Handle the case of either 2 or 3 dims
+              auto winfo = weights.request();
+              double * raw_weights;
+              int64_t nnz;
+              if (winfo.ndim == 2) {
+                  nnz = 1;
+                  raw_weights = extract_buffer <double> (
+                      weights, "weights", 2, temp_shape, {-1, n_samp}
+                  );
+              } else {
+                  raw_weights = extract_buffer <double> (
+                      weights, "weights", 3, temp_shape, {-1, n_samp, -1}
+                  );
+                  nnz = temp_shape[2];
+              }
+
+              int32_t * raw_data_index = extract_buffer <int32_t> (
+                  data_index, "data_index", 1, temp_shape, {n_det}
+              );
+              double * raw_det_data = extract_buffer <double> (
+                  det_data, "det_data", 2, temp_shape, {-1, n_samp}
+              );
+
+              Interval * raw_intervals = extract_buffer <Interval> (
+                  intervals, "intervals", 1, temp_shape, {-1}
+              );
+              int64_t n_view = temp_shape[0];
+
+              int64_t * raw_global2local = extract_buffer <int64_t> (
+                  global2local, "global2local", 1, temp_shape, {-1}
+              );
+              int64_t n_global_submap = temp_shape[0];
+
+              T * raw_mapdata = extract_buffer <T> (
+                  mapdata, "mapdata", 3, temp_shape, {-1, n_pix_submap, nnz}
+              );
+              int64_t n_local_submap = temp_shape[0];
+
+              if (offload) {
+                  #ifdef HAVE_OPENMP_TARGET
+
+                  int64_t * dev_pixels = omgr.device_ptr(raw_pixels);
+                  double * dev_weights = omgr.device_ptr(raw_weights);
+                  double * dev_det_data = omgr.device_ptr(raw_det_data);
+                  Interval * dev_intervals = omgr.device_ptr(raw_intervals);
+                  T * dev_mapdata = omgr.device_ptr(raw_mapdata);
+
+                  # pragma omp target data             \
+                  device(dev)                          \
+                  map(to:                              \
+                  raw_weight_index[0:n_det],           \
+                  raw_pixel_index[0:n_det],            \
+                  raw_data_index[0:n_det],             \
+                  raw_global2local[0:n_global_submap], \
+                  n_view,                              \
+                  n_det,                               \
+                  n_samp,                              \
+                  nnz,                                 \
+                  n_pix_submap,                        \
+                  data_scale,                          \
+                  should_scale,                        \
+                  should_subtract,                     \
+                  should_zero                          \
+                  )
+                  {
+                      # pragma omp target teams distribute collapse(2) \
+                      is_device_ptr(                                   \
+                      dev_pixels,                                      \
+                      dev_weights,                                     \
+                      dev_det_data,                                    \
+                      dev_intervals,                                   \
+                      dev_mapdata                                      \
+                      )
+                      for (int64_t idet = 0; idet < n_det; idet++) {
+                          for (int64_t iview = 0; iview < n_view; iview++) {
+                              # pragma omp parallel
+                              {
+                                  # pragma omp for default(shared)
+                                  for (
+                                      int64_t isamp = dev_intervals[iview].first;
+                                      isamp <= dev_intervals[iview].last;
+                                      isamp++
+                                  ) {
+                                      scan_map_inner <T> (
+                                          raw_pixel_index,
+                                          raw_weight_index,
+                                          raw_data_index,
+                                          raw_global2local,
+                                          dev_det_data,
+                                          dev_pixels,
+                                          dev_weights,
+                                          dev_mapdata,
+                                          data_scale,
+                                          isamp,
+                                          n_samp,
+                                          idet,
+                                          nnz,
+                                          n_pix_submap,
+                                          should_zero,
+                                          should_subtract,
+                                          should_scale
+                                      );
+                                  }
+                              }
+                          }
+                      }
+                  }
+
+                  #endif // ifdef HAVE_OPENMP_TARGET
+              } else {
+                  for (int64_t idet = 0; idet < n_det; idet++) {
+                      for (int64_t iview = 0; iview < n_view; iview++) {
+                          #pragma omp parallel for default(shared)
+                          for (
+                              int64_t isamp = raw_intervals[iview].first;
+                              isamp <= raw_intervals[iview].last;
+                              isamp++
+                          ) {
+                              scan_map_inner <T> (
+                                  raw_pixel_index,
+                                  raw_weight_index,
+                                  raw_data_index,
+                                  raw_global2local,
+                                  raw_det_data,
+                                  raw_pixels,
+                                  raw_weights,
+                                  raw_mapdata,
+                                  data_scale,
+                                  isamp,
+                                  n_samp,
+                                  idet,
+                                  nnz,
+                                  n_pix_submap,
+                                  should_zero,
+                                  should_subtract,
+                                  should_scale
+                              );
+                          }
+                      }
+                  }
+              }
+
+              return;
+          });
+    return;
+}
+
+void init_ops_scan_map(py::module & m) {
+    register_ops_scan_map <double> (m, "ops_scan_map_float64");
+    register_ops_scan_map <float> (m, "ops_scan_map_float32");
+    register_ops_scan_map <int64_t> (m, "ops_scan_map_int64");
+    register_ops_scan_map <int32_t> (m, "ops_scan_map_int32");
+}
diff --git a/src/toast/_libtoast/ops_stokes_weights.cpp b/src/toast/_libtoast/ops_stokes_weights.cpp
index c3c0c45ef..e54dfea98 100644
--- a/src/toast/_libtoast/ops_stokes_weights.cpp
+++ b/src/toast/_libtoast/ops_stokes_weights.cpp
@@ -160,7 +160,10 @@ void init_ops_stokes_weights(py::module & m) {
                 double * dev_quats = omgr.device_ptr(raw_quats);
                 double * dev_weights = omgr.device_ptr(raw_weights);
                 Interval * dev_intervals = omgr.device_ptr(raw_intervals);
-                double * dev_hwp = omgr.device_ptr(raw_hwp);
+                double * dev_hwp = NULL;
+                if (raw_hwp != NULL) {
+                    dev_hwp = omgr.device_ptr(raw_hwp);
+                }
 
                 # pragma omp target data   \
                 device(dev)                \
@@ -172,11 +175,6 @@ void init_ops_stokes_weights(py::module & m) {
                 n_view,                    \
                 n_det,                     \
                 n_samp                     \
-                )                          \
-                use_device_ptr(            \
-                raw_weight_index,          \
-                raw_quat_index,            \
-                raw_epsilon                \
                 )
                 {
                     # pragma omp target teams distribute collapse(2) \
@@ -188,24 +186,27 @@ void init_ops_stokes_weights(py::module & m) {
                     )
                     for (int64_t idet = 0; idet < n_det; idet++) {
                         for (int64_t iview = 0; iview < n_view; iview++) {
-                            # pragma omp parallel for default(shared)
-                            for (
-                                int64_t isamp = dev_intervals[iview].first;
-                                isamp <= dev_intervals[iview].last;
-                                isamp++
-                            ) {
-                                stokes_weights_IQU_inner(
-                                    cal,
-                                    raw_quat_index,
-                                    raw_weight_index,
-                                    dev_quats,
-                                    dev_hwp,
-                                    raw_epsilon,
-                                    dev_weights,
-                                    isamp,
-                                    n_samp,
-                                    idet
-                                );
+                            # pragma omp parallel
+                            {
+                                # pragma omp for default(shared)
+                                for (
+                                    int64_t isamp = dev_intervals[iview].first;
+                                    isamp <= dev_intervals[iview].last;
+                                    isamp++
+                                ) {
+                                    stokes_weights_IQU_inner(
+                                        cal,
+                                        raw_quat_index,
+                                        raw_weight_index,
+                                        dev_quats,
+                                        dev_hwp,
+                                        raw_epsilon,
+                                        dev_weights,
+                                        isamp,
+                                        n_samp,
+                                        idet
+                                    );
+                                }
                             }
                         }
                     }
@@ -285,10 +286,8 @@ void init_ops_stokes_weights(py::module & m) {
                 raw_weight_index[0:n_det], \
                 n_view,                    \
                 n_det,                     \
-                n_samp                     \
-                )                          \
-                use_device_ptr(            \
-                raw_weight_index           \
+                n_samp,                    \
+                cal                        \
                 )
                 {
                     # pragma omp target teams distribute collapse(2) \
@@ -298,15 +297,18 @@ void init_ops_stokes_weights(py::module & m) {
                     )
                     for (int64_t idet = 0; idet < n_det; idet++) {
                         for (int64_t iview = 0; iview < n_view; iview++) {
-                            # pragma omp parallel for default(shared)
-                            for (
-                                int64_t isamp = dev_intervals[iview].first;
-                                isamp <= dev_intervals[iview].last;
-                                isamp++
-                            ) {
-                                int32_t w_indx = raw_weight_index[idet];
-                                int64_t off = (w_indx * n_samp) + isamp;
-                                dev_weights[off] = cal;
+                            # pragma omp parallel
+                            {
+                                # pragma omp for default(shared)
+                                for (
+                                    int64_t isamp = dev_intervals[iview].first;
+                                    isamp <= dev_intervals[iview].last;
+                                    isamp++
+                                ) {
+                                    int32_t w_indx = raw_weight_index[idet];
+                                    int64_t off = (w_indx * n_samp) + isamp;
+                                    dev_weights[off] = cal;
+                                }
                             }
                         }
                     }
diff --git a/src/toast/_libtoast/template_offset.cpp b/src/toast/_libtoast/template_offset.cpp
index 96f734fc7..855725418 100644
--- a/src/toast/_libtoast/template_offset.cpp
+++ b/src/toast/_libtoast/template_offset.cpp
@@ -57,15 +57,16 @@ void init_template_offset(py::module & m) {
                 Interval * dev_intervals = omgr.device_ptr(raw_intervals);
                 double * dev_amplitudes = omgr.device_ptr(raw_amplitudes);
 
-                # pragma omp target data \
-                device(dev)              \
-                map(to                   \
-                : n_view,                \
-                n_samp,                  \
-                data_index,              \
-                step_length,             \
-                amp_offset,              \
-                raw_n_amp_views)
+                # pragma omp target data  \
+                device(dev)               \
+                map(to:                   \
+                n_view,                   \
+                n_samp,                   \
+                data_index,               \
+                step_length,              \
+                amp_offset,               \
+                raw_n_amp_views[0:n_view] \
+                )
                 {
                     int64_t offset = amp_offset;
                     # pragma omp target teams distribute firstprivate(offset) \
@@ -74,17 +75,20 @@ void init_template_offset(py::module & m) {
                     dev_det_data,                                             \
                     dev_intervals)
                     for (int64_t iview = 0; iview < n_view; iview++) {
-                        int64_t first_samp = dev_intervals[iview].first;
-                        # pragma omp parallel for default(shared)
-                        for (
-                            int64_t isamp = dev_intervals[iview].first;
-                            isamp <= dev_intervals[iview].last;
-                            isamp++
-                        ) {
-                            int64_t d = data_index * n_samp + isamp;
-                            int64_t amp = offset +
-                                          (int64_t)((isamp - first_samp) / step_length);
-                            dev_det_data[d] += dev_amplitudes[amp];
+                        # pragma omp parallel
+                        {
+                            # pragma omp for default(shared)
+                            for (
+                                int64_t isamp = dev_intervals[iview].first;
+                                isamp <= dev_intervals[iview].last;
+                                isamp++
+                            ) {
+                                int64_t d = data_index * n_samp + isamp;
+                                int64_t amp = offset + (int64_t)(
+                                    (isamp - dev_intervals[iview].first) / step_length
+                                );
+                                dev_det_data[d] += dev_amplitudes[amp];
+                            }
                         }
                         offset += raw_n_amp_views[iview];
                     }
@@ -94,7 +98,6 @@ void init_template_offset(py::module & m) {
             } else {
                 int64_t offset = amp_offset;
                 for (int64_t iview = 0; iview < n_view; iview++) {
-                    int64_t first_samp = raw_intervals[iview].first;
                     #pragma omp parallel for default(shared)
                     for (
                         int64_t isamp = raw_intervals[iview].first;
@@ -102,8 +105,9 @@ void init_template_offset(py::module & m) {
                         isamp++
                     ) {
                         int64_t d = data_index * n_samp + isamp;
-                        int64_t amp = offset +
-                                      (int64_t)((isamp - first_samp) / step_length);
+                        int64_t amp = offset + (int64_t)(
+                            (isamp - raw_intervals[iview].first) / step_length
+                        );
                         raw_det_data[d] += raw_amplitudes[amp];
                     }
                     offset += raw_n_amp_views[iview];
@@ -155,7 +159,7 @@ void init_template_offset(py::module & m) {
 
             // Optionally use flags
             bool use_flags = false;
-            uint8_t * raw_det_flags = (uint8_t *)omgr.null;
+            uint8_t * raw_det_flags = omgr.null_ptr <uint8_t> ();
             if (flag_index >= 0) {
                 raw_det_flags = extract_buffer <uint8_t> (
                     flag_data, "flag_data", 2, temp_shape, {-1, n_samp}
@@ -171,17 +175,18 @@ void init_template_offset(py::module & m) {
                 Interval * dev_intervals = omgr.device_ptr(raw_intervals);
                 double * dev_amplitudes = omgr.device_ptr(raw_amplitudes);
 
-                # pragma omp target data \
-                device(dev)              \
-                map(to                   \
-                : n_view,                \
-                n_samp,                  \
-                data_index,              \
-                flag_index,              \
-                step_length,             \
-                amp_offset,              \
-                raw_n_amp_views,         \
-                use_flags)
+                # pragma omp target data   \
+                device(dev)                \
+                map(to:                    \
+                n_view,                    \
+                n_samp,                    \
+                data_index,                \
+                flag_index,                \
+                step_length,               \
+                amp_offset,                \
+                raw_n_amp_views[0:n_view], \
+                use_flags                  \
+                )
                 {
                     int64_t offset = amp_offset;
                     # pragma omp target teams distribute firstprivate(offset) \
@@ -191,28 +196,31 @@ void init_template_offset(py::module & m) {
                     dev_det_flags,                                            \
                     dev_intervals)
                     for (int64_t iview = 0; iview < n_view; iview++) {
-                        int64_t first_samp = dev_intervals[iview].first;
-                        # pragma omp parallel for default(shared)
-                        for (
-                            int64_t isamp = dev_intervals[iview].first;
-                            isamp <= dev_intervals[iview].last;
-                            isamp++
-                        ) {
-                            int64_t d = data_index * n_samp + isamp;
-                            int64_t amp = offset +
-                                          (int64_t)((isamp - first_samp) / step_length);
-                            double contrib = 0.0;
-                            if (use_flags) {
-                                int64_t f = flag_index * n_samp + isamp;
-                                uint8_t check = dev_det_flags[f] & flag_mask;
-                                if (check == 0) {
-                                    contrib = dev_det_data[d];
+                        # pragma omp parallel
+                        {
+                            # pragma omp for default(shared)
+                            for (
+                                int64_t isamp = dev_intervals[iview].first;
+                                isamp <= dev_intervals[iview].last;
+                                isamp++
+                            ) {
+                                int64_t d = data_index * n_samp + isamp;
+                                int64_t amp = offset + (int64_t)(
+                                    (isamp - dev_intervals[iview].first) / step_length
+                                );
+                                double contrib = 0.0;
+                                if (use_flags) {
+                                    int64_t f = flag_index * n_samp + isamp;
+                                    uint8_t check = dev_det_flags[f] & flag_mask;
+                                    if (check == 0) {
+                                        contrib = dev_det_data[d];
+                                    }
+                                } else {
+                                    contrib = raw_det_data[d];
                                 }
-                            } else {
-                                contrib = raw_det_data[d];
+                                # pragma omp atomic update
+                                dev_amplitudes[amp] += contrib;
                             }
-                            # pragma omp atomic
-                            dev_amplitudes[amp] += contrib;
                         }
                         offset += raw_n_amp_views[iview];
                     }
@@ -222,7 +230,6 @@ void init_template_offset(py::module & m) {
             } else {
                 int64_t offset = amp_offset;
                 for (int64_t iview = 0; iview < n_view; iview++) {
-                    int64_t first_samp = raw_intervals[iview].first;
                     #pragma omp parallel for default(shared)
                     for (
                         int64_t isamp = raw_intervals[iview].first;
@@ -230,8 +237,9 @@ void init_template_offset(py::module & m) {
                         isamp++
                     ) {
                         int64_t d = data_index * n_samp + isamp;
-                        int64_t amp = offset +
-                                      (int64_t)((isamp - first_samp) / step_length);
+                        int64_t amp = offset + (int64_t)(
+                            (isamp - raw_intervals[iview].first) / step_length
+                        );
                         double contrib = 0.0;
                         if (use_flags) {
                             int64_t f = flag_index * n_samp + isamp;
@@ -242,7 +250,7 @@ void init_template_offset(py::module & m) {
                         } else {
                             contrib = raw_det_data[d];
                         }
-                        #pragma omp atomic
+                        #pragma omp atomic update
                         raw_amplitudes[amp] += contrib;
                     }
                     offset += raw_n_amp_views[iview];
@@ -296,10 +304,13 @@ void init_template_offset(py::module & m) {
                     dev_amp_out,        \
                     dev_offset_var)
                     {
-                        # pragma omp parallel for default(shared)
-                        for (int64_t iamp = 0; iamp < n_amp; iamp++) {
-                            dev_amp_out[iamp] = dev_amp_in[iamp];
-                            dev_amp_out[iamp] *= dev_offset_var[iamp];
+                        # pragma omp parallel
+                        {
+                            # pragma omp for default(shared)
+                            for (int64_t iamp = 0; iamp < n_amp; iamp++) {
+                                dev_amp_out[iamp] = dev_amp_in[iamp];
+                                dev_amp_out[iamp] *= dev_offset_var[iamp];
+                            }
                         }
                     }
                 }
@@ -315,72 +326,5 @@ void init_template_offset(py::module & m) {
             return;
         });
 
-    // m.def(
-    //     "template_offset_add_to_signal", [](int64_t step_length, py::buffer
-    // amplitudes,
-    //                                         py::buffer data) {
-    //         pybuffer_check_1D <double> (amplitudes);
-    //         pybuffer_check_1D <double> (data);
-    //         py::buffer_info info_amplitudes = amplitudes.request();
-    //         py::buffer_info info_data = data.request();
-    //         int64_t n_amp = info_amplitudes.size;
-    //         int64_t n_data = info_data.size;
-    //         double * raw_amplitudes = reinterpret_cast <double *>
-    // (info_amplitudes.ptr);
-    //         double * raw_data = reinterpret_cast <double *> (info_data.ptr);
-    //         toast::template_offset_add_to_signal(step_length, n_amp, raw_amplitudes,
-    //                                              n_data, raw_data);
-    //         return;
-    //     }, py::arg("step_length"), py::arg("amplitudes"), py::arg(
-    //         "data"), R"(
-    //     Accumulate offset amplitudes to timestream data.
-
-    //     Each amplitude value is accumulated to `step_length` number of samples.  The
-    //     final offset will be at least this many samples, but may be more if the step
-    //     size does not evenly divide into the number of samples.
-
-    //     Args:
-    //         step_length (int64):  The minimum number of samples for each offset.
-    //         amplitudes (array):  The float64 amplitude values.
-    //         data (array):  The float64 timestream values to accumulate.
-
-    //     Returns:
-    //         None.
-
-    // )");
-
-    // m.def(
-    //     "template_offset_project_signal", [](int64_t step_length, py::buffer data,
-    //                                          py::buffer amplitudes) {
-    //         pybuffer_check_1D <double> (amplitudes);
-    //         pybuffer_check_1D <double> (data);
-    //         py::buffer_info info_amplitudes = amplitudes.request();
-    //         py::buffer_info info_data = data.request();
-    //         int64_t n_amp = info_amplitudes.size;
-    //         int64_t n_data = info_data.size;
-    //         double * raw_amplitudes = reinterpret_cast <double *>
-    // (info_amplitudes.ptr);
-    //         double * raw_data = reinterpret_cast <double *> (info_data.ptr);
-    //         toast::template_offset_project_signal(step_length, n_data, raw_data,
-    //                                               n_amp, raw_amplitudes);
-    //         return;
-    //     }, py::arg("step_length"), py::arg("data"), py::arg(
-    //         "amplitudes"), R"(
-    //     Accumulate timestream data into offset amplitudes.
-
-    //     Chunks of `step_length` number of samples are accumulated into the offset
-    //     amplitudes.  If step_length does not evenly divide into the total number of
-    //     samples, the final amplitude will be extended to include the remainder.
-
-    //     Args:
-    //         step_length (int64):  The minimum number of samples for each offset.
-    //         data (array):  The float64 timestream values.
-    //         amplitudes (array):  The float64 amplitude values.
-
-    //     Returns:
-    //         None.
-
-    // )");
-
     return;
 }
diff --git a/src/toast/_libtoast/todmap_scanning.cpp b/src/toast/_libtoast/todmap_scanning.cpp
index b68dd3d50..4cb4a5b93 100644
--- a/src/toast/_libtoast/todmap_scanning.cpp
+++ b/src/toast/_libtoast/todmap_scanning.cpp
@@ -5,6 +5,8 @@
 
 #include <module.hpp>
 
+// FIXME:  This whole file should be removed and functionality placed
+// in ops_scan_map.cpp
 
 template <typename T>
 void register_scan_map(py::module & m, char const * name) {
diff --git a/src/toast/accelerator/__init__.py b/src/toast/accelerator/__init__.py
index 02759bdc6..ce069b5b7 100644
--- a/src/toast/accelerator/__init__.py
+++ b/src/toast/accelerator/__init__.py
@@ -8,6 +8,7 @@
     AcceleratorObject,
     accel_assign_device,
     accel_data_create,
+    accel_data_reset,
     accel_data_delete,
     accel_data_present,
     accel_data_update_device,
diff --git a/src/toast/accelerator/accel.py b/src/toast/accelerator/accel.py
index 293b2e9e9..5088c031f 100644
--- a/src/toast/accelerator/accel.py
+++ b/src/toast/accelerator/accel.py
@@ -7,6 +7,7 @@
 from .._libtoast import Logger
 from .._libtoast import accel_assign_device as omp_accel_assign_device
 from .._libtoast import accel_create as omp_accel_create
+from .._libtoast import accel_reset as omp_accel_reset
 from .._libtoast import accel_delete as omp_accel_delete
 from .._libtoast import accel_enabled as omp_accel_enabled
 from .._libtoast import accel_get_device as omp_accel_get_device
@@ -33,6 +34,7 @@
 if ("TOAST_GPU_JAX" in os.environ) and (os.environ["TOAST_GPU_JAX"] in enable_vals):
     try:
         import jax
+        import jax.numpy as jnp
 
         from ..jax.device import jax_accel_assign_device, jax_accel_get_device
         from ..jax.intervals import INTERVALS_JAX
@@ -80,25 +82,29 @@ def accel_get_device():
         return -1
 
 
-def accel_assign_device(node_procs, node_rank, disabled):
+def accel_assign_device(node_procs, node_rank, mem_gb, disabled):
     """
     Assign processes to target devices.
 
-    NOTE:
-    One can pick devices visible to processes using Slurm and teh following commands
-    `--gpus-per-task=1 --gpu-bind=single:1`
+    NOTE for NERSC:
+    One can pick devices visible to processes using Slurm and the following
+    commands:
+
+        `--gpus-per-task=1 --gpu-bind=single:1`
 
     Args:
         node_procs (int): number of processes per node
         node_rank (int): rank of the current process, within the node
+        mem_gb (float):  The GPU device memory in GB.
         disabled (bool): gpu computing is disabled
 
     Returns:
         None: the device is stored in a backend specific global variable
+
     """
-    # FIXME some functions (such as poiting_detector) require the omp device to have been assigned
-    # so it should be called even when using JAX or running on CPU
-    omp_accel_assign_device(node_procs, node_rank, disabled)
+    # This function should always be called to avoid extra logic
+    # in the compiled kernels.
+    omp_accel_assign_device(node_procs, node_rank, mem_gb, disabled)
     if use_accel_jax:
         jax_accel_assign_device(node_procs, node_rank, disabled)
 
@@ -133,6 +139,7 @@ def accel_data_present(data):
         log.warning("Accelerator support not enabled, data not present")
         return False
 
+
 @function_timer
 def accel_data_create(data):
     """Create device buffers.
@@ -157,6 +164,31 @@ def accel_data_create(data):
         log = Logger.get()
         log.warning("Accelerator support not enabled, cannot create")
 
+
+@function_timer
+def accel_data_reset(data):
+    """Reset device buffers.
+
+    Using the input data array, reset to zero the corresponding device array.
+    For OpenMP target offload, this runs a small device kernel that sets the
+    memory to zero.  For jax arrays, this is a no-op, since those
+    arrays are mapped and managed elsewhere.
+
+    Args:
+        data (array):  The host array.
+    Returns:
+        None
+
+    """
+    if use_accel_omp:
+        omp_accel_reset(data)
+    elif use_accel_jax:
+        data.data = jnp.zeroes_like(data.data)
+    else:
+        log = Logger.get()
+        log.warning("Accelerator support not enabled, cannot reset")
+
+
 @function_timer
 def accel_data_update_device(data):
     """Update device buffers.
@@ -182,6 +214,7 @@ def accel_data_update_device(data):
         log.warning("Accelerator support not enabled, not updating device")
         return None
 
+
 @function_timer
 def accel_data_update_host(data):
     """Update host buffers.
@@ -207,6 +240,7 @@ def accel_data_update_host(data):
         log = Logger.get()
         log.warning("Accelerator support not enabled, not updating host")
 
+
 @function_timer
 def accel_data_delete(data):
     """Delete device copy of the data.
@@ -390,3 +424,23 @@ def accel_delete(self):
             raise RuntimeError(msg)
         self._accel_delete()
         self._accel_used = False
+
+    def _accel_reset(self):
+        msg = f"The _accel_reset function was not defined for this class."
+        raise RuntimeError(msg)
+
+    def accel_reset(self):
+        """Reset to zero the data on the accelerator.
+
+        Returns:
+            None
+
+        """
+        if not accel_enabled():
+            return
+        if not self.accel_in_use():
+            log = Logger.get()
+            msg = f"Device data not in use, cannot reset"
+            log.error(msg)
+            raise RuntimeError(msg)
+        self._accel_reset()
diff --git a/src/toast/accelerator/data_localization.py b/src/toast/accelerator/data_localization.py
index 0a4161422..8688e70b1 100644
--- a/src/toast/accelerator/data_localization.py
+++ b/src/toast/accelerator/data_localization.py
@@ -138,4 +138,5 @@ def f_result(*args, **kwargs):
 
 def display_datamovement():
     """Displays all recorded data movement so far"""
-    print(dataMovementTracker)
+    log = Logger.get()
+    log.info(str(dataMovementTracker))
diff --git a/src/toast/accelerator/kernel_registry.py b/src/toast/accelerator/kernel_registry.py
index c4a4e7ea0..ff2751c0a 100644
--- a/src/toast/accelerator/kernel_registry.py
+++ b/src/toast/accelerator/kernel_registry.py
@@ -80,7 +80,10 @@ def dispatch(
         if (
             (impl == ImplementationType.JAX) or (impl == ImplementationType.COMPILED)
         ) and accel_enabled():
-            registry[name][impl] = function_timer(function_datamovementtracker(f))
+            # FIXME:  disabled for now because it does not work for arbitrary
+            # data objects (e.g. scalars) when running with omp target offload.
+            # registry[name][impl] = function_timer(function_datamovementtracker(f))
+            registry[name][impl] = function_timer(f)
         else:
             registry[name][impl] = function_timer(f)
 
diff --git a/src/toast/data.py b/src/toast/data.py
index c75283cfb..d3ec9a2d8 100644
--- a/src/toast/data.py
+++ b/src/toast/data.py
@@ -690,8 +690,9 @@ def accel_delete(self, names):
         for key in names["global"]:
             val = self._internal.get(key, None)
             if isinstance(val, AcceleratorObject):
-                log.verbose(f"Calling Data accel_delete for '{key}'")
-                val.accel_delete()
+                if val.accel_exists():
+                    log.verbose(f"Calling Data accel_delete for '{key}'")
+                    val.accel_delete()
             else:
                 msg = f"Data accel_delete: '{key}' ({type(val)}) "
                 msg += "is not an AcceleratorObject"
@@ -708,7 +709,8 @@ def accel_clear(self):
 
         for key, val in self._internal.items():
             if isinstance(val, AcceleratorObject):
-                val.accel_delete()
+                if val.accel_exists():
+                    val.accel_delete()
             else:
                 msg = f"Data accel_clear: '{key}' ({type(val)}) "
                 msg += "is not an AcceleratorObject"
diff --git a/src/toast/mpi.py b/src/toast/mpi.py
index bdb336b4a..b1aa189d9 100644
--- a/src/toast/mpi.py
+++ b/src/toast/mpi.py
@@ -49,19 +49,22 @@ class MPI_Comm:
 if use_accel_omp or use_accel_jax:
     node_procs = 1
     node_rank = 0
+    accel_gb = 2.0
+    if "TOAST_GPU_MEM_GB" in os.environ:
+        accel_gb = float(os.environ["TOAST_GPU_MEM_GB"])
     if use_mpi:
         # We need to compute which process goes to which device
         nodecomm = MPI.COMM_WORLD.Split_type(MPI.COMM_TYPE_SHARED, 0)
         node_procs = nodecomm.size
         node_rank = nodecomm.rank
-        accel_assign_device(node_procs, node_rank, False)
+        accel_assign_device(node_procs, node_rank, accel_gb, False)
         nodecomm.Free()
         del nodecomm
     else:
-        accel_assign_device(node_procs, node_rank, False)
+        accel_assign_device(node_procs, node_rank, accel_gb, False)
 else:
     # Disabled == True
-    accel_assign_device(1, 0, True)
+    accel_assign_device(1, 0, 0.0, True)
 
 # We put other imports *after* the MPI check, since usually the MPI initialization
 # is time sensitive and may timeout the job if it does not happen quickly enough.
diff --git a/src/toast/observation_data.py b/src/toast/observation_data.py
index 9e0992e7b..2ab742cd2 100644
--- a/src/toast/observation_data.py
+++ b/src/toast/observation_data.py
@@ -12,6 +12,7 @@
 from .accelerator import (
     AcceleratorObject,
     accel_data_create,
+    accel_data_reset,
     accel_data_delete,
     accel_data_present,
     accel_data_update_device,
@@ -263,14 +264,7 @@ def change_detectors(self, detectors):
             self._data = self._flatdata.reshape(self._shape)
             # move things to GPU if needed
             if should_update_GPU:
-                if use_accel_jax:
-                    # NOTE: we create an array of zeroes directly on GPU to avoid useless data movement
-                    self._data = MutableJaxArray(
-                        self._data, gpu_data=jax.numpy.zeros_like(self._data)
-                    )
-                else:
-                    # FIXME: Should we zero the device memory?
-                    pass
+                self.accel_reset()
             realloced = False
         return realloced
 
@@ -535,6 +529,18 @@ def _accel_delete(self):
         elif use_accel_jax:
             self._data = accel_data_delete(self._data)
 
+    def _accel_reset(self):
+        log = Logger.get()
+        log.verbose(f"DetectorData _accel_reset")
+        if use_accel_omp:
+            accel_data_reset(self._raw)
+        elif use_accel_jax:
+            # NOTE: we create an array of zeroes directly on GPU to avoid
+            # useless data movement
+            self._data = MutableJaxArray(
+                self._data, gpu_data=jax.numpy.zeros_like(self._data)
+            )
+
 
 class DetDataManager(MutableMapping):
     """Class used to manage DetectorData objects in an Observation.
diff --git a/src/toast/ops/mapmaker.py b/src/toast/ops/mapmaker.py
index ff585cac9..d09945090 100644
--- a/src/toast/ops/mapmaker.py
+++ b/src/toast/ops/mapmaker.py
@@ -375,7 +375,7 @@ def _exec(self, data, detectors=None, use_accel=None, **kwargs):
             out_cleaned = self.clean_name
             if self.save_cleaned and self.overwrite_cleaned:
                 # Modify data in place
-                out_cleaned = self.det_data
+                out_cleaned = None
 
             amplitudes_apply = ApplyAmplitudes(
                 op="subtract",
diff --git a/src/toast/ops/mapmaker_templates.py b/src/toast/ops/mapmaker_templates.py
index 26a6edb2f..0ba4e04a4 100644
--- a/src/toast/ops/mapmaker_templates.py
+++ b/src/toast/ops/mapmaker_templates.py
@@ -168,9 +168,6 @@ def reset_templates(self):
     def _exec(self, data, detectors=None, use_accel=None, **kwargs):
         log = Logger.get()
 
-        # Kernel selection
-        implementation = self.select_kernels(use_accel=use_accel)
-
         # Check that the detector data is set
         if self.det_data is None:
             raise RuntimeError("You must set the det_data trait before calling exec()")
@@ -193,7 +190,9 @@ def _exec(self, data, detectors=None, use_accel=None, **kwargs):
         if not self._initialized:
             if use_accel:
                 # fail when a user tries to run the initialization pipeline on GPU
-                raise RuntimeError("You cannot currently initialize templates on device (please disable accel for this operator/pipeline).") 
+                raise RuntimeError(
+                    "You cannot currently initialize templates on device (please disable accel for this operator/pipeline)."
+                )
             for tmpl in self.templates:
                 tmpl.view = self.view
                 tmpl.det_data_units = self.det_data_units
@@ -265,8 +264,13 @@ def _exec(self, data, detectors=None, use_accel=None, **kwargs):
                 )
                 ob.detdata[self.det_data].update_units(self.det_data_units)
 
-                for d in dets:
-                    ob.detdata[self.det_data][d, :] = 0
+                if use_accel:
+                    # We are running on the accelerator, so our output data must exist
+                    # on the device and will be used there.
+                    ob.detdata[self.det_data].accel_reset()
+                else:
+                    for d in dets:
+                        ob.detdata[self.det_data][d, :] = 0
 
             for d in all_dets:
                 for tmpl in self.templates:
@@ -922,7 +926,7 @@ def _exec(self, data, detectors=None, **kwargs):
         self.template_matrix.det_flags = save_tmpl_flags
         self.template_matrix.det_flag_mask = save_tmpl_mask
         # FIXME: this reset does not seem needed
-        #if not self.mc_mode:
+        # if not self.mc_mode:
         #    self.template_matrix.reset_templates()
 
         memreport.prefix = "End of amplitude solve"
@@ -1037,7 +1041,9 @@ def _exec(self, data, detectors=None, use_accel=None, **kwargs):
 
         if self.output is not None:
             # We just copy the input here, since it will be overwritten
-            Copy(detdata=[(self.det_data, self.output)]).apply(data, use_accel=use_accel)
+            Copy(detdata=[(self.det_data, self.output)]).apply(
+                data, use_accel=use_accel
+            )
 
         # Projecting amplitudes to timestreams
         self.template_matrix.transpose = False
@@ -1052,8 +1058,6 @@ def _exec(self, data, detectors=None, use_accel=None, **kwargs):
         else:
             combine.result = self.output
 
-        apply_input = np.array(data.obs[0].detdata[self.det_data][:])
-
         # Project and operate, one detector at a time
         pipe = Pipeline(
             detector_sets=["SINGLE"],
diff --git a/src/toast/ops/mapmaker_utils/mapmaker_utils.py b/src/toast/ops/mapmaker_utils/mapmaker_utils.py
index 316484377..d9e4b6833 100644
--- a/src/toast/ops/mapmaker_utils/mapmaker_utils.py
+++ b/src/toast/ops/mapmaker_utils/mapmaker_utils.py
@@ -101,7 +101,7 @@ def _exec(self, data, detectors=None, use_accel=None, **kwargs):
         log = Logger.get()
 
         # Kernel selection
-        implementation = self.select_kernels(use_accel=use_accel)
+        implementation, use_accel = self.select_kernels(use_accel=use_accel)
 
         if self.pixel_dist is None:
             raise RuntimeError(
@@ -321,7 +321,7 @@ def _exec(self, data, detectors=None, use_accel=None, **kwargs):
         log = Logger.get()
 
         # Kernel selection
-        implementation = self.select_kernels(use_accel=use_accel)
+        implementation, use_accel = self.select_kernels(use_accel=use_accel)
 
         if self.pixel_dist is None:
             raise RuntimeError(
@@ -613,7 +613,7 @@ def _exec(self, data, detectors=None, use_accel=None, **kwargs):
         log = Logger.get()
 
         # Kernel selection
-        implementation = self.select_kernels(use_accel=use_accel)
+        implementation, use_accel = self.select_kernels(use_accel=use_accel)
 
         if self.pixel_dist is None:
             raise RuntimeError(
@@ -756,9 +756,7 @@ def _exec(self, data, detectors=None, use_accel=None, **kwargs):
                 pix_indx,
                 ob.detdata[self.pixels].data,
                 weight_indx,
-                ob.detdata[self.weights].data.reshape(
-                    (n_weight_dets, ob.n_local_samples, -1)
-                ),
+                ob.detdata[self.weights].data,
                 data_indx,
                 ob.detdata[self.det_data].data,
                 flag_indx,
@@ -788,12 +786,6 @@ def _finalize(self, data, use_accel=None, **kwargs):
                 data[self.zmap].sync_alltoallv()
             else:
                 data[self.zmap].sync_allreduce()
-            if use_accel:
-                log.verbose_rank(
-                    f"Operator {self.name} finalize calling zmap update device",
-                    comm=data.comm.comm_group,
-                )
-                data[self.zmap].accel_update_device()
 
             zmap_good = data[self.zmap].data[:, :, 0] != 0.0
             zmap_min = np.zeros((data[self.zmap].n_value), dtype=np.float64)
@@ -811,6 +803,13 @@ def _finalize(self, data, use_accel=None, **kwargs):
                 for m in range(data[self.zmap].n_value):
                     msg += f"    map {m} {zmap_min[m]:1.3e} ... {zmap_max[m]:1.3e}"
                 log.debug(msg)
+
+            if use_accel:
+                log.verbose_rank(
+                    f"Operator {self.name} finalize calling zmap update device",
+                    comm=data.comm.comm_group,
+                )
+                data[self.zmap].accel_update_device()
         return
 
     def _requires(self):
diff --git a/src/toast/ops/noise_weight/kernels.py b/src/toast/ops/noise_weight/kernels.py
index deb2d8bc6..682e88377 100644
--- a/src/toast/ops/noise_weight/kernels.py
+++ b/src/toast/ops/noise_weight/kernels.py
@@ -4,7 +4,7 @@
 
 import numpy as np
 
-# from ..._libtoast import noise_weight as libtoast_noise_weight
+from ..._libtoast import noise_weight as libtoast_noise_weight
 from ...accelerator import ImplementationType, kernel, use_accel_jax
 from .kernels_numpy import noise_weight_numpy
 
@@ -12,6 +12,11 @@
     from .kernels_jax import noise_weight_jax
 
 
+@kernel(impl=ImplementationType.COMPILED, name="noise_weight")
+def noise_weight_compiled(*args, use_accel=False):
+    return libtoast_noise_weight(*args, use_accel)
+
+
 @kernel(impl=ImplementationType.DEFAULT)
 def noise_weight(
     det_data,
@@ -33,19 +38,10 @@ def noise_weight(
         None
 
     """
-    # FIXME to be replaced by a *direct* call to libtoast_noise_weight once it is ported
-    return noise_weight(
+    return libtoast_noise_weight(
         det_data,
         det_data_index,
         intervals,
         detector_weights,
-        impl=ImplementationType.NUMPY,
-        use_accel=use_accel,
+        use_accel,
     )
-
-
-# FIXME:  Add this to _libtoast
-
-# @kernel(impl=ImplementationType.COMPILED, name="noise_weight")
-# def noise_weight_compiled(*args, use_accel=False):
-#     return libtoast_noise_weight(*args, use_accel)
diff --git a/src/toast/ops/noise_weight/noise_weight.py b/src/toast/ops/noise_weight/noise_weight.py
index 121bec1e3..24b10a06b 100644
--- a/src/toast/ops/noise_weight/noise_weight.py
+++ b/src/toast/ops/noise_weight/noise_weight.py
@@ -47,7 +47,7 @@ def _exec(self, data, detectors=None, use_accel=None, **kwargs):
         log = Logger.get()
 
         # Kernel selection
-        implementation = self.select_kernels(use_accel=use_accel)
+        implementation, use_accel = self.select_kernels(use_accel=use_accel)
 
         for ob in data.obs:
             # Get the detectors we are using for this observation
@@ -69,10 +69,13 @@ def _exec(self, data, detectors=None, use_accel=None, **kwargs):
 
             # Compute the noise for each detector (using the correct units)
             noise = ob[self.noise_model]
-            detector_weights = [
-                noise.detector_weight(detector).to(data_invcov_units).value
-                for detector in dets
-            ]
+            detector_weights = np.array(
+                [
+                    noise.detector_weight(detector).to(data_invcov_units).value
+                    for detector in dets
+                ],
+                dtype=np.float64,
+            )
 
             # Multiply detectors by their respective noise weight
             intervals = ob.intervals[self.view].data
diff --git a/src/toast/ops/pixels_healpix/pixels_healpix.py b/src/toast/ops/pixels_healpix/pixels_healpix.py
index cd014133b..7a8d8bc37 100644
--- a/src/toast/ops/pixels_healpix/pixels_healpix.py
+++ b/src/toast/ops/pixels_healpix/pixels_healpix.py
@@ -148,7 +148,7 @@ def _exec(self, data, detectors=None, use_accel=None, **kwargs):
         log = Logger.get()
 
         # Kernel selection
-        implementation = self.select_kernels(use_accel=use_accel)
+        implementation, use_accel = self.select_kernels(use_accel=use_accel)
 
         if self.detector_pointing is None:
             raise RuntimeError("The detector_pointing trait must be set")
@@ -157,13 +157,7 @@ def _exec(self, data, detectors=None, use_accel=None, **kwargs):
             self._local_submaps = np.zeros(self._n_submap, dtype=np.uint8)
 
         # Expand detector pointing
-        if self.quats is not None:
-            quats_name = self.quats
-        else:
-            if self.detector_pointing.quats is not None:
-                quats_name = self.detector_pointing.quats
-            else:
-                quats_name = defaults.quats
+        quats_name = self.detector_pointing.quats
 
         view = self.view
         if view is None:
@@ -171,7 +165,6 @@ def _exec(self, data, detectors=None, use_accel=None, **kwargs):
             view = self.detector_pointing.view
 
         # Expand detector pointing
-        self.detector_pointing.quats = quats_name
         self.detector_pointing.apply(data, detectors=detectors, use_accel=use_accel)
 
         for ob in data.obs:
@@ -251,9 +244,9 @@ def _exec(self, data, detectors=None, use_accel=None, **kwargs):
                     ob.detdata.accel_create(self.pixels)
 
             # Get the flags if needed.  Use the same flags as
-            # detector pointing.  Use empty array in absence of flags
+            # detector pointing.
             if self.detector_pointing.shared_flags is None:
-                flags = np.array([], dtype=np.uint8)
+                flags = np.zeros(1, dtype=np.uint8)
             else:
                 flags = ob.shared[self.detector_pointing.shared_flags].data
 
@@ -322,11 +315,8 @@ def _requires(self):
         return req
 
     def _provides(self):
-        # FIXME quats generated are keyed in self.detector_pointing.quats (and not self.quats)
-        prov = {
-            "detdata": [self.pixels, self.detector_pointing.quats],
-            "global": list(),
-        }
+        prov = self.detector_pointing.provides()
+        prov["detdata"].append(self.pixels)
         if self.create_dist is not None:
             prov["global"].append(self.create_dist)
         return prov
@@ -340,4 +330,9 @@ def _implementations(self):
         ]
 
     def _supports_accel(self):
-        return True
+        if (self.detector_pointing is not None) and (
+            self.detector_pointing.supports_accel()
+        ):
+            return True
+        else:
+            return False
diff --git a/src/toast/ops/pointing_detector/pointing_detector.py b/src/toast/ops/pointing_detector/pointing_detector.py
index 98a5139c4..74df6b60a 100644
--- a/src/toast/ops/pointing_detector/pointing_detector.py
+++ b/src/toast/ops/pointing_detector/pointing_detector.py
@@ -90,7 +90,7 @@ def _exec(self, data, detectors=None, use_accel=None, **kwargs):
         log = Logger.get()
 
         # Kernel selection
-        implementation = self.select_kernels(use_accel=use_accel)
+        implementation, use_accel = self.select_kernels(use_accel=use_accel)
 
         coord_rot = None
         if self.coord_in is None:
diff --git a/src/toast/ops/polyfilter/polyfilter.py b/src/toast/ops/polyfilter/polyfilter.py
index 898fce659..6b5619e60 100644
--- a/src/toast/ops/polyfilter/polyfilter.py
+++ b/src/toast/ops/polyfilter/polyfilter.py
@@ -112,7 +112,7 @@ def _exec(self, data, detectors=None, use_accel=None, **kwargs):
         gt = GlobalTimers.get()
 
         # Kernel selection
-        implementation = self.select_kernels(use_accel=use_accel)
+        implementation, use_accel = self.select_kernels(use_accel=use_accel)
 
         if detectors is not None:
             raise RuntimeError("PolyFilter2D cannot be run on subsets of detectors")
@@ -493,7 +493,7 @@ def _exec(self, data, detectors=None, use_accel=None, **kwargs):
         log = Logger.get()
 
         # Kernel selection
-        implementation = self.select_kernels(use_accel=use_accel)
+        implementation, use_accel = self.select_kernels(use_accel=use_accel)
 
         if self.pattern is None:
             pat = None
@@ -758,7 +758,7 @@ def _exec(self, data, detectors=None, use_accel=None, **kwargs):
         pat = re.compile(self.pattern)
 
         # Kernel selection
-        implementation = self.select_kernels(use_accel=use_accel)
+        implementation, use_accel = self.select_kernels(use_accel=use_accel)
 
         for obs in data.obs:
             comm, temp_ob = self._redistribute(data, obs, timer, log)
diff --git a/src/toast/ops/scan_map/kernels.py b/src/toast/ops/scan_map/kernels.py
index 49a02fe29..58212f02f 100644
--- a/src/toast/ops/scan_map/kernels.py
+++ b/src/toast/ops/scan_map/kernels.py
@@ -6,11 +6,10 @@
 import numpy as np
 
 from ... import qarray as qa
-from ..._libtoast import AlignedF64, AlignedI64
-from ..._libtoast import scan_map_float32 as libtoast_scan_map_float32
-from ..._libtoast import scan_map_float64 as libtoast_scan_map_float64
-from ..._libtoast import scan_map_int32 as libtoast_scan_map_int32
-from ..._libtoast import scan_map_int64 as libtoast_scan_map_int64
+from ..._libtoast import ops_scan_map_float32 as libtoast_scan_map_float32
+from ..._libtoast import ops_scan_map_float64 as libtoast_scan_map_float64
+from ..._libtoast import ops_scan_map_int32 as libtoast_scan_map_int32
+from ..._libtoast import ops_scan_map_int64 as libtoast_scan_map_int64
 from ...accelerator import ImplementationType, kernel, use_accel_jax
 from .kernels_numpy import scan_map_numpy
 
@@ -30,10 +29,10 @@ def scan_map(
     weights,
     weight_index,
     intervals,
-    data_scale=1.0,
-    should_zero=False,
-    should_subtract=False,
-    should_scale=False,
+    data_scale,
+    should_zero,
+    should_subtract,
+    should_scale,
     use_accel=False,
 ):
     """Kernel for scanning a map into timestreams.
@@ -64,7 +63,7 @@ def scan_map(
         None.
 
     """
-    return scan_map_pseudocompiled(
+    return scan_map_compiled(
         global2local,
         n_pix_submap,
         mapdata,
@@ -79,14 +78,16 @@ def scan_map(
         should_zero,
         should_subtract,
         should_scale,
-        use_accel)
+        use_accel=use_accel,
+    )
+
 
 @kernel(impl=ImplementationType.COMPILED, name="scan_map")
 def scan_map_compiled(*args, use_accel=False):
-    return scan_map_pseudocompiled(*args, use_accel)
+    return scan_map_compiled(*args, use_accel)
+
 
-# FIXME:  This "compiled" kernel will migrate fully to the _libtoast extension.
-def scan_map_pseudocompiled(
+def scan_map_compiled(
     global2local,
     n_pix_submap,
     mapdata,
@@ -97,10 +98,10 @@ def scan_map_pseudocompiled(
     weights,
     weight_index,
     intervals,
-    data_scale=1.0,
-    should_zero=False,
-    should_subtract=False,
-    should_scale=False,
+    data_scale,
+    should_zero,
+    should_subtract,
+    should_scale,
     use_accel=False,
 ):
     if mapdata.dtype.char == "d":
@@ -116,55 +117,28 @@ def scan_map_pseudocompiled(
         msg += f"type '{mapdata.dtype.char}'"
         raise NotImplementedError(msg)
 
-    nmap = 1
-    if len(weights.shape) > 2:
-        nmap = weights.shape[2]
-    for idet in range(len(det_data_index)):
-        pix = pixels[pixels_index[idet]]
-        wts = weights[weight_index[idet]]
-        tod = det_data[det_data_index[idet]]
-        for view in intervals:
-            view_samples = view.last + 1 - view.first
-            vslice = slice(view.first, view.last + 1)
-
-            # Temp buffers
-            local_tod_raw = AlignedF64.zeros(view_samples)
-            local_tod = local_tod_raw.array()
-            local_submap_pix_raw = AlignedI64.zeros(view_samples)
-            local_submap_pix = local_submap_pix_raw.array()
-            local_submap_pix[:] = -1
-            local_submap_raw = AlignedI64.zeros(view_samples)
-            local_submap = local_submap_raw.array()
-            local_submap[:] = -1
-
-            # Get local submaps and pixel indices within each submap
-            good = pix[vslice] >= 0
-            local_submap_pix[good] = pix[vslice][good] % n_pix_submap
-            local_submap[good] = global2local[pix[vslice][good] // n_pix_submap]
-
-            fcomp(
-                n_pix_submap,
-                nmap,
-                local_submap,
-                local_submap_pix,
-                mapdata.reshape((-1,)),
-                wts[vslice].reshape((-1,)),
-                local_tod,
-            )
+    wt_view = weights
+    if len(weights.shape) == 2:
+        # One value per sample, but our kernel always assumes a
+        # 3D array.
+        wdets = weights.shape[0]
+        wsamps = weights.shape[1]
+        wt_view = weights.reshape((wdets, wsamps, -1))
 
-            # Accumulate
-            if should_zero:
-                tod[vslice] = 0
-            if should_subtract:
-                tod[vslice][good] -= local_tod[good]
-            elif should_scale:
-                tod[vslice][good] *= local_tod[good]
-            else:
-                tod[vslice][good] += local_tod[good]
-
-            del local_submap
-            del local_submap_raw
-            del local_submap_pix
-            del local_submap_pix_raw
-            del local_tod
-            del local_tod_raw
+    fcomp(
+        global2local,
+        n_pix_submap,
+        mapdata,
+        det_data,
+        det_data_index,
+        pixels,
+        pixels_index,
+        wt_view,
+        weight_index,
+        intervals,
+        data_scale,
+        should_zero,
+        should_subtract,
+        should_scale,
+        use_accel,
+    )
diff --git a/src/toast/ops/scan_map/scan_map.py b/src/toast/ops/scan_map/scan_map.py
index 3d5760414..ec842c51a 100644
--- a/src/toast/ops/scan_map/scan_map.py
+++ b/src/toast/ops/scan_map/scan_map.py
@@ -70,7 +70,7 @@ def _exec(self, data, detectors=None, use_accel=None, **kwargs):
         log = Logger.get()
 
         # Kernel selection
-        implementation = self.select_kernels(use_accel=use_accel)
+        implementation, use_accel = self.select_kernels(use_accel=use_accel)
 
         # Check that the detector data is set
         if self.det_data is None:
@@ -144,10 +144,10 @@ def _exec(self, data, detectors=None, use_accel=None, **kwargs):
                 weights,
                 weight_indx,
                 intervals,
-                data_scale=data_scale,
-                should_zero=self.zero,
-                should_subtract=self.subtract,
-                should_scale=False,
+                data_scale,
+                self.zero,
+                self.subtract,
+                False,
                 impl=implementation,
                 use_accel=use_accel,
             )
@@ -235,7 +235,7 @@ def _exec(self, data, detectors=None, use_accel=None, **kwargs):
         log = Logger.get()
 
         # Kernel selection
-        implementation = self.select_kernels(use_accel=use_accel)
+        implementation, use_accel = self.select_kernels(use_accel=use_accel)
 
         # Check that the detector data is set
         if self.det_flags is None:
@@ -346,7 +346,7 @@ def _exec(self, data, detectors=None, use_accel=None, **kwargs):
         log = Logger.get()
 
         # Kernel selection
-        implementation = self.select_kernels(use_accel=use_accel)
+        implementation, use_accel = self.select_kernels(use_accel=use_accel)
 
         # Check that the detector data is set
         if self.det_data is None:
@@ -407,10 +407,10 @@ def _exec(self, data, detectors=None, use_accel=None, **kwargs):
                 weights,
                 weight_indx,
                 intervals,
-                data_scale=1.0,
-                should_zero=False,
-                should_subtract=False,
-                should_scale=True,
+                1.0,
+                False,
+                False,
+                True,
                 impl=implementation,
                 use_accel=use_accel,
             )
diff --git a/src/toast/ops/stokes_weights/kernels.py b/src/toast/ops/stokes_weights/kernels.py
index 60e767b06..57af09a30 100644
--- a/src/toast/ops/stokes_weights/kernels.py
+++ b/src/toast/ops/stokes_weights/kernels.py
@@ -36,7 +36,7 @@ def stokes_weights_I(
 
     Args:
         weight_index (array):  The index into the weights array for each detector.
-        weights (array):  The array of I, Q, and U weights at each sample for each
+        weights (array):  The array of intensity weights at each sample for each
             detector.
         intervals (array):  The array of sample intervals.
         use_accel (bool):  Whether to use the accelerator for this call (if supported).
diff --git a/src/toast/ops/stokes_weights/stokes_weights.py b/src/toast/ops/stokes_weights/stokes_weights.py
index e2671d61c..c67f49618 100644
--- a/src/toast/ops/stokes_weights/stokes_weights.py
+++ b/src/toast/ops/stokes_weights/stokes_weights.py
@@ -144,7 +144,7 @@ def _exec(self, data, detectors=None, use_accel=None, **kwargs):
         self._nnz = len(self.mode)
 
         # Kernel selection
-        implementation = self.select_kernels(use_accel=use_accel)
+        implementation, use_accel = self.select_kernels(use_accel=use_accel)
 
         if self.detector_pointing is None:
             raise RuntimeError("The detector_pointing trait must be set")
@@ -243,7 +243,7 @@ def _exec(self, data, detectors=None, use_accel=None, **kwargs):
 
             if self.mode == "IQU":
                 if self.hwp_angle is None:
-                    hwp_data = np.zeros((0,), dtype=np.float64)
+                    hwp_data = np.zeros(1, dtype=np.float64)
                 else:
                     hwp_data = ob.shared[self.hwp_angle].data
                 stokes_weights_IQU(
@@ -295,4 +295,9 @@ def _implementations(self):
         ]
 
     def _supports_accel(self):
-        return True
+        if (self.detector_pointing is not None) and (
+            self.detector_pointing.supports_accel()
+        ):
+            return True
+        else:
+            return False
diff --git a/src/toast/pixels_io_healpix.py b/src/toast/pixels_io_healpix.py
index 54b4d1627..878711dbd 100644
--- a/src/toast/pixels_io_healpix.py
+++ b/src/toast/pixels_io_healpix.py
@@ -268,12 +268,12 @@ def collect_healpix_submaps(pix, comm_bytes=10000000):
 
 @function_timer
 def write_healpix_fits(
-        pix,
-        path,
-        nest=True,
-        comm_bytes=10000000,
-        report_memory=False,
-        single_precision=False,
+    pix,
+    path,
+    nest=True,
+    comm_bytes=10000000,
+    report_memory=False,
+    single_precision=False,
 ):
     """Write pixel data to a HEALPix format FITS table.
 
diff --git a/src/toast/scripts/toast_benchmark_ground.py b/src/toast/scripts/toast_benchmark_ground.py
index 3d09e7611..0e22b3b28 100644
--- a/src/toast/scripts/toast_benchmark_ground.py
+++ b/src/toast/scripts/toast_benchmark_ground.py
@@ -175,6 +175,8 @@ def parse_arguments():
             write_invcov=False,
             write_rcond=False,
             keep_final_products=True,
+            save_cleaned=True,
+            overwrite_cleaned=True,
         ),
         toast.ops.PixelsHealpix(name="pixels_final", view="scanning", enabled=False),
         toast.ops.BinMap(
diff --git a/src/toast/scripts/toast_benchmark_satellite.py b/src/toast/scripts/toast_benchmark_satellite.py
index c06380962..4da4297d1 100644
--- a/src/toast/scripts/toast_benchmark_satellite.py
+++ b/src/toast/scripts/toast_benchmark_satellite.py
@@ -128,6 +128,8 @@ def parse_arguments():
             write_invcov=False,
             write_rcond=False,
             keep_final_products=True,
+            save_cleaned=True,
+            overwrite_cleaned=True,
         ),
         toast.ops.PixelsHealpix(name="pixels_final", enabled=False),
         toast.ops.BinMap(
@@ -353,6 +355,7 @@ def main():
     # display information on GPU data movement when running in debug/verbose mode
     display_datamovement()
 
+
 if __name__ == "__main__":
     world, procs, rank = toast.mpi.get_world()
     python_startup_time(rank)
diff --git a/src/toast/scripts/toast_merge_timings.py b/src/toast/scripts/toast_merge_timings.py
index 0175b7300..71c02509c 100755
--- a/src/toast/scripts/toast_merge_timings.py
+++ b/src/toast/scripts/toast_merge_timings.py
@@ -16,16 +16,18 @@
 import glob
 import pandas as pd
 
+
 def find_csv_files(folder, file_pattern):
     """
     Find all CSV files matching the specified pattern in the given folder and its subfolders.
-    
+
     :param folder: The root folder to start searching for CSV files.
     :param file_pattern: The pattern of the CSV files to search for.
     :return: A list of file paths matching the specified pattern.
     """
     return glob.glob(os.path.join(folder, file_pattern), recursive=True)
 
+
 def process_timer_path(s):
     """
     Process the input string and return three strings: simplified path, kernel type,
@@ -35,40 +37,41 @@ def process_timer_path(s):
     :return: A tuple containing the simplified path, kernel type, and operation name.
     """
     # Remove occurrences of '(function) ' and '._exec'
-    simplified_path = s.replace('(function) ', '').replace('._exec', '')
+    simplified_path = s.replace("(function) ", "").replace("._exec", "")
 
     # Determine kernel type
-    if s.endswith('_jax'):
-        kernel_type = 'JAX'
-        simplified_path = simplified_path.replace('_jax', '')
-    elif s.endswith('_compiled'):
-        kernel_type = 'COMPILED'
-        simplified_path = simplified_path.replace('_compiled', '')
-    elif s.endswith('_numpy'):
-        kernel_type = 'NUMPY'
-        simplified_path = simplified_path.replace('_numpy', '')
-    elif ('accel_data' in simplified_path) or ('INTERVALS_JAX' in simplified_path):
-        kernel_type = 'DATA_MOVEMENT'
-    elif '|dispatch|' in simplified_path:
-        kernel_type = 'DEFAULT'
+    if s.endswith("_jax"):
+        kernel_type = "JAX"
+        simplified_path = simplified_path.replace("_jax", "")
+    elif s.endswith("_compiled"):
+        kernel_type = "COMPILED"
+        simplified_path = simplified_path.replace("_compiled", "")
+    elif s.endswith("_numpy"):
+        kernel_type = "NUMPY"
+        simplified_path = simplified_path.replace("_numpy", "")
+    elif ("accel_data" in simplified_path) or ("INTERVALS_JAX" in simplified_path):
+        kernel_type = "DATA_MOVEMENT"
+    elif "|dispatch|" in simplified_path:
+        kernel_type = "DEFAULT"
     else:
         kernel_type = None
 
     # Extract operation name
-    operation_name = simplified_path.split('|')[-1]
-    if kernel_type == 'DATA_MOVEMENT':
+    operation_name = simplified_path.split("|")[-1]
+    if kernel_type == "DATA_MOVEMENT":
         # Name clean-up specific to INTERVALS_JAX data movement operations
-        if operation_name == 'INTERVALS_JAX.__init__':
-            operation_name = 'accel_data_update_device'
-        elif operation_name == 'INTERVALS_JAX.to_host':
-            operation_name = 'accel_data_update_host'
+        if operation_name == "INTERVALS_JAX.__init__":
+            operation_name = "accel_data_update_device"
+        elif operation_name == "INTERVALS_JAX.to_host":
+            operation_name = "accel_data_update_host"
 
     # Add '_kernel' at the end of simplified_path if it contains '|dispatch|'
-    if '|dispatch|' in simplified_path:
-        simplified_path += '_kernel'
+    if "|dispatch|" in simplified_path:
+        simplified_path += "_kernel"
 
     return simplified_path, kernel_type, operation_name
 
+
 def load_csv_files(file_paths):
     """
     Load CSV files from the given file paths, keeping only the specified columns
@@ -84,14 +87,16 @@ def load_csv_files(file_paths):
 
         # Process the index and extract simplified_path, kernel_type, and operation_name
         processed_indices = list(map(process_timer_path, df.index))
-        simplified_paths, kernel_types, operation_names = zip(*processed_indices) # unzip
+        simplified_paths, kernel_types, operation_names = zip(
+            *processed_indices
+        )  # unzip
 
         # Replace index with simplified_path
         df.index = simplified_paths
 
         # Add kernel_type and operation_name columns
-        df['kernel_type'] = kernel_types
-        df['operation_name'] = operation_names
+        df["kernel_type"] = kernel_types
+        df["operation_name"] = operation_names
 
         # Rename the 'Mean Time' column to the folder name
         df = df.rename(columns={"Mean Time": folder_name})
@@ -100,6 +105,7 @@ def load_csv_files(file_paths):
 
     return dataframes
 
+
 def combine_kernel_types(k1, k2):
     """
     Combine two kernel types into a single kernel type,
@@ -109,18 +115,19 @@ def combine_kernel_types(k1, k2):
     :param k2: The second kernel type (string or None).
     :return: The combined kernel type (string).
     """
-    if (k1 is None):
-        return k2 
+    if k1 is None:
+        return k2
     if (k2 is None) or (k1 == k2):
         return k1
-    return 'MULTIPLE'
+    return "MULTIPLE"
+
 
 def merge_dataframes(dataframes):
     """
     Merge the given DataFrames using the union of all row names and keeping only
     the specified columns. The columns are intermeshed as best as possible based
     on their original order in the input CSV files.
-    
+
     :param dataframes: A dictionary of DataFrames to merge.
     :return: A merged DataFrame with the union of all row names.
     """
@@ -129,17 +136,37 @@ def merge_dataframes(dataframes):
 
     # Merge the remaining DataFrames one-by-one
     for df in dataframes[1:]:
-        merged_df = pd.merge(merged_df, df, left_index=True, right_index=True, how='outer', suffixes=('_1', '_2'))
+        merged_df = pd.merge(
+            merged_df,
+            df,
+            left_index=True,
+            right_index=True,
+            how="outer",
+            suffixes=("_1", "_2"),
+        )
 
         # Combine kernel_type and operation_name columns
-        merged_df['kernel_type'] = merged_df['kernel_type_1'].combine(merged_df['kernel_type_2'], combine_kernel_types)
-        merged_df['operation_name'] = merged_df['operation_name_1'].combine_first(merged_df['operation_name_2'])
+        merged_df["kernel_type"] = merged_df["kernel_type_1"].combine(
+            merged_df["kernel_type_2"], combine_kernel_types
+        )
+        merged_df["operation_name"] = merged_df["operation_name_1"].combine_first(
+            merged_df["operation_name_2"]
+        )
 
         # Drop extra kernel_type and operation_name columns
-        merged_df.drop(columns=['kernel_type_1', 'kernel_type_2', 'operation_name_1', 'operation_name_2'], inplace=True)
+        merged_df.drop(
+            columns=[
+                "kernel_type_1",
+                "kernel_type_2",
+                "operation_name_1",
+                "operation_name_2",
+            ],
+            inplace=True,
+        )
 
     return merged_df
 
+
 def merge_kernel_rows(df):
     """
     Filters the input DataFrame to keep only rows where the `kernel_type` column is 'SEVERAL',
@@ -150,19 +177,22 @@ def merge_kernel_rows(df):
     :return: A filtered and grouped DataFrame.
     """
     # Filter DataFrame to keep only rows where kernel_type is 'MULTIPLE' or 'DATA_MOVEMENT'
-    df_filtered = df[df['kernel_type'].isin(['MULTIPLE', 'DATA_MOVEMENT'])]
+    df_filtered = df[df["kernel_type"].isin(["MULTIPLE", "DATA_MOVEMENT"])]
 
     # Drop the kernel_type column
-    df_filtered = df_filtered.drop(columns=['kernel_type'])
+    df_filtered = df_filtered.drop(columns=["kernel_type"])
 
     # Group by operation_name and sum the values of rows that get collapsed together
     # Some columns contains None/nan, hence the need for numeric_only=False
-    df_grouped = df_filtered.groupby('operation_name').apply(lambda x: x.sum(numeric_only=False))
+    df_grouped = df_filtered.groupby("operation_name").apply(
+        lambda x: x.sum(numeric_only=False)
+    )
     # Remove additional operation_name column added by the non-numeric-only summation
-    df_grouped.drop(columns=['operation_name'], inplace=True)
+    df_grouped.drop(columns=["operation_name"], inplace=True)
 
     return df_grouped
 
+
 if __name__ == "__main__":
     folder = "."
     file_pattern = "**/timing.csv"
@@ -177,7 +207,7 @@ def merge_kernel_rows(df):
     print(f"Merged kernel data saved to '{output_kernels_file}'")
 
     # Save the merged DataFrame to a CSV file
-    merged_df = merged_df.drop(columns=['kernel_type', 'operation_name'])
+    merged_df = merged_df.drop(columns=["kernel_type", "operation_name"])
     output_file = "merged_timings.csv"
     merged_df.to_csv(output_file)
     print(f"Merged data saved to '{output_file}'")
diff --git a/src/toast/templates/offset/kernels.py b/src/toast/templates/offset/kernels.py
index bdb6dbf14..ee545c586 100644
--- a/src/toast/templates/offset/kernels.py
+++ b/src/toast/templates/offset/kernels.py
@@ -62,8 +62,9 @@ def offset_add_to_signal(
         amplitudes,
         data_index,
         det_data,
-        intervals, 
-        use_accel)
+        intervals,
+        use_accel,
+    )
 
 
 @kernel(impl=ImplementationType.DEFAULT)
diff --git a/src/toast/templates/offset/kernels_jax.py b/src/toast/templates/offset/kernels_jax.py
index f967a4b29..f0ecaa65a 100644
--- a/src/toast/templates/offset/kernels_jax.py
+++ b/src/toast/templates/offset/kernels_jax.py
@@ -315,7 +315,8 @@ def offset_project_signal_jax(
         n_amp_views,
         intervals.first,
         intervals.last,
-        intervals_max_length)
+        intervals_max_length,
+    )
 
 
 @kernel(impl=ImplementationType.JAX, name="offset_apply_diag_precond")
diff --git a/src/toast/templates/offset/offset.py b/src/toast/templates/offset/offset.py
index 69fcb6c98..10458922c 100644
--- a/src/toast/templates/offset/offset.py
+++ b/src/toast/templates/offset/offset.py
@@ -494,7 +494,7 @@ def _add_to_signal(self, detector, amplitudes, use_accel=None, **kwargs):
         log = Logger.get()
 
         # Kernel selection
-        implementation = self.select_kernels(use_accel=use_accel)
+        implementation, use_accel = self.select_kernels(use_accel=use_accel)
 
         amp_offset = self._det_start[detector]
         for iob, ob in enumerate(self.data.obs):
@@ -528,7 +528,7 @@ def _project_signal(self, detector, amplitudes, use_accel=None, **kwargs):
         log = Logger.get()
 
         # Kernel selection
-        implementation = self.select_kernels(use_accel=use_accel)
+        implementation, use_accel = self.select_kernels(use_accel=use_accel)
 
         amp_offset = self._det_start[detector]
         for iob, ob in enumerate(self.data.obs):
@@ -642,7 +642,7 @@ def _apply_precond(self, amplitudes_in, amplitudes_out, use_accel=None, **kwargs
             # preconditioner is just the application of offset variance.
 
             # Kernel selection
-            implementation = self.select_kernels(use_accel=use_accel)
+            implementation, use_accel = self.select_kernels(use_accel=use_accel)
 
             offset_apply_diag_precond(
                 self._offsetvar,
diff --git a/src/toast/traits.py b/src/toast/traits.py
index 4c8b99433..526268447 100644
--- a/src/toast/traits.py
+++ b/src/toast/traits.py
@@ -514,25 +514,27 @@ def select_kernels(self, use_accel=None):
         This returns the kernel implementation that should be used
 
         Returns:
-            (ImplementationType):  The implementation type.
+            (tuple):  The (ImplementationType, bool use_accel) switches.
 
         """
         impls = self.implementations()
-        if use_accel:
+        if use_accel is None:
+            return ImplementationType.DEFAULT, False
+        elif use_accel:
             if use_accel_jax:
                 if ImplementationType.JAX not in impls:
                     msg = f"JAX accelerator use is enabled, "
                     msg += f"but not supported by {self.name}"
                     raise RuntimeError(msg)
-                return ImplementationType.JAX
+                return ImplementationType.JAX, True
             else:
                 if ImplementationType.COMPILED not in impls:
                     msg = f"OpenMP accelerator use is enabled, "
                     msg += f"but not supported by {self.name}"
                     raise RuntimeError(msg)
-                return ImplementationType.COMPILED
+                return ImplementationType.COMPILED, True
         else:
-            return ImplementationType.DEFAULT
+            return ImplementationType.DEFAULT, False
 
     def __eq__(self, other):
         if len(self.traits()) != len(other.traits()):
diff --git a/src/toast/utils.py b/src/toast/utils.py
index bd57153ae..d40152ed3 100644
--- a/src/toast/utils.py
+++ b/src/toast/utils.py
@@ -797,7 +797,7 @@ def __isub__(self, other):
         -= operation performing set difference on all keys
         `other` can be a normal dict
         """
-        for (key, value) in other.items():
+        for key, value in other.items():
             self[key] -= set(value)
         return self
 
@@ -806,7 +806,7 @@ def __ior__(self, other):
         |= operation performing set union on all keys
         `other` can be a normal dict
         """
-        for (key, value) in other.items():
+        for key, value in other.items():
             self[key] |= set(value)
         return self
 
@@ -815,14 +815,14 @@ def __iand__(self, other):
         &= operation performing set intersection on all keys
         `other` can be a normal dict
         """
-        for (key, value) in other.items():
+        for key, value in other.items():
             self[key] &= set(value)
         return self
 
     def __str__(self):
         """prints only the non-empty/None sets for brevity sake"""
         result = "{ "
-        for (k, v) in self.items():
+        for k, v in self.items():
             if (len(v) > 0) and not all(x is None for x in v):
                 result += f"{k}:{list(v)} "
         result += "}"
@@ -830,7 +830,7 @@ def __str__(self):
 
     def is_empty(self):
         """returns True if the container is empty or contains only None"""
-        for (k, v) in self.items():
+        for k, v in self.items():
             if (len(v) > 0) and not all(x is None for x in v):
                 return False
         return True


Test patch:
diff --git a/etc/pyomptarget/test.cpp b/etc/pyomptarget/test.cpp
new file mode 100644
index 000000000..7f214f72e
--- /dev/null
+++ b/etc/pyomptarget/test.cpp
@@ -0,0 +1,11 @@
+
+#include <iostream>
+#include <omp.h>
+
+int main() {
+  double s;
+#pragma omp target teams distribute parallel for reduction(+:s) map(tofrom:s)
+  for (int idx = 0; idx < 1000; ++idx) s+= idx;
+  std::cout << s << std::endl;
+}
+
diff --git a/src/toast/tests/_helpers.py b/src/toast/tests/_helpers.py
index bf1db88a9..210a69300 100644
--- a/src/toast/tests/_helpers.py
+++ b/src/toast/tests/_helpers.py
@@ -148,7 +148,7 @@ def create_boresight_telescope(group_size, sample_rate=10.0 * u.Hz):
         det_table[idet]["psd_fmin"] = 1.0e-5 * u.Hz
         det_table[idet]["psd_fknee"] = 0.05 * u.Hz
         det_table[idet]["psd_alpha"] = 1.0
-        det_table[idet]["psd_net"] = 100 * (u.K * np.sqrt(1.0 * u.second))
+        det_table[idet]["psd_net"] = 0.001 * (u.K * np.sqrt(1.0 * u.second))
 
     fp = Focalplane(
         detector_data=det_table,
@@ -409,11 +409,17 @@ def create_healpix_ring_satellite(mpicomm, obs_per_group=1, nside=64):
     for obs in range(obs_per_group):
         oname = "test-{}-{}".format(toastcomm.group, obs)
         oid = obs_per_group * toastcomm.group + obs
-        tele = create_space_telescope(toastcomm.group_size)
+        tele = create_space_telescope(
+            toastcomm.group_size, 
+            sample_rate=rate * u.Hz,
+            pixel_per_process=4,
+        )
 
         # Move all detectors to the boresight
         for row in tele.focalplane.detector_data:
             row["quat"] = np.array([0, 0, 0, 1], dtype=np.float64)
+            row["psd_net"] = 1.0 * (u.K * np.sqrt(1.0 * u.second))
+            row["psd_fknee"] = 0.0 * u.Hz
 
         # FIXME: for full testing we should set detranks as approximately the sqrt
         # of the grid size so that we test the row / col communicators.
@@ -446,7 +452,9 @@ def create_healpix_ring_satellite(mpicomm, obs_per_group=1, nside=64):
             dtype=np.float64,
         )
         ob.detdata.create(defaults.det_data, dtype=np.float64, units=u.K)
+        ob.detdata[defaults.det_data][:] = 0
         ob.detdata.create(defaults.det_flags, dtype=np.uint8)
+        ob.detdata[defaults.det_flags][:] = 0
         # Rank zero of each grid column creates the data
         stamps = None
         position = None
diff --git a/src/toast/tests/ops_mapmaker.py b/src/toast/tests/ops_mapmaker.py
index 56bba9579..192a94622 100644
--- a/src/toast/tests/ops_mapmaker.py
+++ b/src/toast/tests/ops_mapmaker.py
@@ -125,11 +125,11 @@ def test_offset(self):
             solve_rcond_threshold=1.0e-1,
             map_rcond_threshold=1.0e-1,
             write_hits=False,
-            write_map=False,
+            write_map=True,
             write_cov=False,
             write_rcond=False,
-            keep_solver_products=True,
-            keep_final_products=True,
+            keep_solver_products=False,
+            keep_final_products=False,
             output_dir=testdir,
         )
 
@@ -137,7 +137,12 @@ def test_offset(self):
         mapper.apply(data)
 
         # Check that we can also run in full-memory mode
-        use_accel = False
+        tmatrix.reset()
+        mapper.reset_pix_dist = True
+        pixels.apply(data)
+        weights.apply(data)
+
+        use_accel = None
         if accel_enabled() and (
             pixels.supports_accel()
             and weights.supports_accel()
@@ -151,8 +156,6 @@ def test_offset(self):
             data.accel_update_device(weights.requires())
             data.accel_update_device(mapper.requires())
 
-        pixels.apply(data, use_accel=use_accel)
-        weights.apply(data, use_accel=use_accel)
         binner.full_pointing = True
         mapper.name = "test2"
         mapper.apply(data, use_accel=use_accel)
diff --git a/src/toast/tests/ops_mapmaker_utils.py b/src/toast/tests/ops_mapmaker_utils.py
index 448054cba..cf36e39a1 100644
--- a/src/toast/tests/ops_mapmaker_utils.py
+++ b/src/toast/tests/ops_mapmaker_utils.py
@@ -9,6 +9,7 @@
 from astropy import units as u
 
 from .. import ops as ops
+from ..accelerator import accel_enabled
 from ..mpi import MPI
 from ..noise import Noise
 from ..observation import default_values as defaults
@@ -272,6 +273,32 @@ def test_zmap(self):
             build_zmap_corr.apply(data)
             zmap_corr[stype] = data[build_zmap_corr.zmap]
 
+        # Test one case on the accelerator if supported
+
+        build_zmap = ops.BuildNoiseWeighted(
+            pixel_dist="pixel_dist",
+            noise_model="noise_model",
+            det_data="noise",
+            zmap="zmap_gpu",
+            sync_type="allreduce",
+        )
+
+        use_accel = None
+        if accel_enabled() and (
+            pixels.supports_accel()
+            and weights.supports_accel()
+            and build_zmap.supports_accel()
+        ):
+            use_accel = True
+            data.accel_create(pixels.requires())
+            data.accel_create(weights.requires())
+            data.accel_create(build_zmap.requires())
+            data.accel_update_device(pixels.requires())
+            data.accel_update_device(weights.requires())
+            data.accel_update_device(build_zmap.requires())
+            build_zmap.apply(data)
+            zmap["gpu"] = data[build_zmap.zmap]
+
         # Manual check
 
         invnpp_units = 1.0 / (u.K**2)
@@ -342,4 +369,24 @@ def test_zmap(self):
                 failed = comm.allreduce(failed, op=MPI.LOR)
             self.assertFalse(failed)
 
+        if "gpu" in zmap:
+            comm = zmap["gpu"].distribution.comm
+
+            failed = False
+            for sm in range(zmap["gpu"].distribution.n_local_submap):
+                for px in range(zmap["gpu"].distribution.n_pix_submap):
+                    if zmap["gpu"].data[sm, px, 0] != 0:
+                        if not np.allclose(
+                            zmap["gpu"].data[sm, px],
+                            check_zmap.data[sm, px],
+                            atol=1.0e-6,
+                        ):
+                            print(
+                                f"zmap(gpu)[{sm}, {px}] = {zmap['gpu'].data[sm, px]}, check = {check_zmap.data[sm, px]}"
+                            )
+                            failed = True
+            if comm is not None:
+                failed = comm.allreduce(failed, op=MPI.LOR)
+            self.assertFalse(failed)
+
         close_data(data)
diff --git a/src/toast/tests/ops_pointing_healpix.py b/src/toast/tests/ops_pointing_healpix.py
index 2ebbc2f0c..e0eefc9af 100644
--- a/src/toast/tests/ops_pointing_healpix.py
+++ b/src/toast/tests/ops_pointing_healpix.py
@@ -10,7 +10,7 @@
 from .. import ops as ops
 from .. import qarray as qa
 from .._libtoast import healpix_pixels, stokes_weights
-from ..accelerator import ImplementationType
+from ..accelerator import ImplementationType, accel_enabled
 from ..healpix import HealpixPixels
 from ..intervals import IntervalList, interval_dtype
 from ..observation import default_values as defaults
@@ -255,7 +255,7 @@ def test_pixweight_pipe(self):
             detector_pointing=detpointing,
         )
         weights = ops.StokesWeights(
-            mode="IQU",
+            mode="I",
             hwp_angle=defaults.hwp_angle,
             detector_pointing=detpointing,
         )
@@ -293,6 +293,59 @@ def test_pixweight_pipe(self):
 
         close_data(data)
 
+    def test_pixweight_accel(self):
+        if not accel_enabled():
+            print("Accelerator not enabled, skipping test")
+            return
+
+        # Create a fake satellite data set for testing
+        data = create_satellite_data(self.comm)
+
+        detpointing = ops.PointingDetectorSimple()
+        pixels = ops.PixelsHealpix(
+            nside=64,
+            detector_pointing=detpointing,
+        )
+        weights = ops.StokesWeights(
+            mode="I",
+            hwp_angle=defaults.hwp_angle,
+            detector_pointing=detpointing,
+        )
+
+        pipe = ops.Pipeline(operators=[pixels, weights])
+
+        # Move data
+        # data.accel_create(pixels.requires())
+        # data.accel_create(weights.requires())
+        # data.accel_update_device(pixels.requires())
+        # data.accel_update_device(weights.requires())
+
+        # Compute
+        pipe.apply(data, use_accel=None)
+
+        # Move inputs back
+        data.accel_update_host(pixels.requires())
+        data.accel_update_host(weights.requires())
+
+        # Also make a copy using a python codepath
+        detpointing.kernel_implementation = ImplementationType.NUMPY
+        detpointing.quats = "pyquat"
+        pixels.kernel_implementation = ImplementationType.NUMPY
+        pixels.quats = "pyquat"
+        pixels.pixels = "pypixels"
+        pixels.apply(data)
+        weights.kernel_implementation = ImplementationType.NUMPY
+        weights.quats = "pyquat"
+        weights.weights = "pyweight"
+        weights.apply(data)
+
+        for ob in data.obs:
+            np.testing.assert_allclose(
+                ob.detdata[defaults.weights], ob.detdata["pyweight"]
+            )
+            np.testing.assert_equal(ob.detdata[defaults.pixels], ob.detdata["pypixels"])
+        close_data(data)
+
     def test_hpix_interval(self):
         data = create_satellite_data(self.comm)
 
diff --git a/src/toast/tests/ops_sim_tod_atm.py b/src/toast/tests/ops_sim_tod_atm.py
index 80e64beed..2f71a7164 100644
--- a/src/toast/tests/ops_sim_tod_atm.py
+++ b/src/toast/tests/ops_sim_tod_atm.py
@@ -123,6 +123,9 @@ def test_sim(self):
             binning=binner,
             output_dir=self.outdir,
             keep_final_products=True,
+            solve_rcond_threshold=1.0e-2,
+            map_rcond_threshold=1.0e-2,
+            write_noiseweighted_map=True,
         )
         mapmaker.apply(data)
 
@@ -171,6 +174,39 @@ def test_sim(self):
             plt.savefig(outfile)
             plt.close()
 
+            mapfile = os.path.join(self.outdir, f"{mapmaker.name}_noiseweighted_map.fits")
+            mdata = hp.read_map(mapfile, None, nest=False)
+            mdata[mdata == 0] = hp.UNSEEN
+
+            outfile = "{}.png".format(mapfile)
+            fig = plt.figure(figsize=[8 * 3, 12])
+            hp.gnomview(
+                mdata[0],
+                xsize=1600,
+                sub=[1, 3, 1],
+                rot=(42.0, -42.0),
+                reso=0.5,
+                nest=False,
+            )
+            hp.gnomview(
+                mdata[1],
+                xsize=1600,
+                sub=[1, 3, 2],
+                rot=(42.0, -42.0),
+                reso=0.5,
+                nest=False,
+            )
+            hp.gnomview(
+                mdata[2],
+                xsize=1600,
+                sub=[1, 3, 3],
+                rot=(42.0, -42.0),
+                reso=0.5,
+                nest=False,
+            )
+            plt.savefig(outfile)
+            plt.close()
+
             good = mdata[0] != hp.UNSEEN
             p = np.sqrt(mdata[1] ** 2 + mdata[2] ** 2)
             pfrac = np.median(p[good] / mdata[0][good])
diff --git a/src/toast/tests/ops_sim_tod_dipole.py b/src/toast/tests/ops_sim_tod_dipole.py
index 09e91067e..5fc89206c 100644
--- a/src/toast/tests/ops_sim_tod_dipole.py
+++ b/src/toast/tests/ops_sim_tod_dipole.py
@@ -109,7 +109,7 @@ def test_sim(self):
             pixel_pointing=pixels,
             stokes_weights=weights,
             noise_model=default_model.noise_model,
-            rcond_threshold=1.0e-6,
+            rcond_threshold=1.0e-2,
             sync_type="alltoallv",
         )
         cov_and_hits.apply(data)
@@ -153,7 +153,7 @@ def test_sim(self):
             plt.close()
 
             mapfile = os.path.join(self.outdir, "toast_bin.fits")
-            mdata = hp.read_map(mapfile, nest=False)
+            mdata = hp.read_map(mapfile, field=0, nest=False)
 
             outfile = "{}.png".format(mapfile)
             hp.mollview(mdata, xsize=1600, nest=False)
@@ -165,11 +165,11 @@ def test_sim(self):
 
             minmap = np.min(mdata)
             maxmap = np.max(mdata)
-            np.testing.assert_almost_equal(maxmap, self.dip_check, decimal=5)
-            np.testing.assert_almost_equal(minmap, -self.dip_check, decimal=5)
-
             minloc = np.argmin(mdata)
             maxloc = np.argmax(mdata)
+
+            np.testing.assert_almost_equal(maxmap, self.dip_check, decimal=5)
+            np.testing.assert_almost_equal(minmap, -self.dip_check, decimal=5)
             np.testing.assert_equal(minloc, self.dip_min_pix)
             np.testing.assert_equal(maxloc, self.dip_max_pix)
 
diff --git a/src/toast/tests/template_offset.py b/src/toast/tests/template_offset.py
index 82d1a9740..db1cd1f42 100644
--- a/src/toast/tests/template_offset.py
+++ b/src/toast/tests/template_offset.py
@@ -9,7 +9,7 @@
 from astropy import units as u
 
 from .. import ops
-from ..accelerator import ImplementationType
+from ..accelerator import ImplementationType, accel_enabled
 from ..observation import default_values as defaults
 from ..templates import Offset
 from ..utils import rate_from_times
@@ -41,6 +41,7 @@ def test_projection(self):
 
         tmpl = Offset(
             det_data=defaults.det_data,
+            det_flags=None,
             times=defaults.times,
             noise_model=noise_model.noise_model,
             step_time=step_seconds * u.second,
@@ -91,6 +92,90 @@ def test_projection(self):
 
         close_data(data)
 
+    def test_accel(self):
+        if not accel_enabled():
+            print("Accelerator use not enabled, skipping test")
+            return
+
+        # Create a fake satellite data set for testing
+        data = create_satellite_data(self.comm)
+
+        # Create a default noise model
+        noise_model = ops.DefaultNoiseModel()
+        noise_model.apply(data)
+
+        # Use 1/10 of an observation as the baseline length.  Make it not evenly
+        # divisible in order to test handling of the final amplitude.
+        ob_time = (
+            data.obs[0].shared[defaults.times][-1]
+            - data.obs[0].shared[defaults.times][0]
+        )
+        step_seconds = float(int(ob_time / 10.0))
+
+        tmpl = Offset(
+            det_data=defaults.det_data,
+            det_flags=None,
+            times=defaults.times,
+            noise_model=noise_model.noise_model,
+            step_time=step_seconds * u.second,
+        )
+        # Set the data
+        tmpl.data = data
+
+        # Get some amplitudes and set to one
+        amps = tmpl.zeros()
+        amps.local[:] = 1.0
+
+        data_names = {
+            "detdata": [defaults.det_data],
+            "shared": [],
+            "global": [],
+            "meta": [],
+            "intervals": [None],
+        }
+
+        data.accel_create(data_names)
+        data.accel_update_device(data_names)
+        amps.accel_create()
+        amps.accel_update_device()
+
+        # Project.
+        for det in tmpl.detectors():
+            for ob in data.obs:
+                tmpl.add_to_signal(det, amps, use_accel=True)
+
+        # Accumulate amplitudes
+        for det in tmpl.detectors():
+            for ob in data.obs:
+                tmpl.project_signal(det, amps, use_accel=True)
+
+        data.accel_update_host(data_names)
+        amps.accel_update_host()
+
+        # Verify
+        for ob in data.obs:
+            # Get the step boundaries
+            (rate, dt, dt_min, dt_max, dt_std) = rate_from_times(
+                ob.shared[defaults.times]
+            )
+            step_samples = int(step_seconds * rate)
+            n_step = ob.n_local_samples // step_samples
+            slices = [
+                slice(x * step_samples, (x + 1) * step_samples, 1)
+                for x in range(n_step - 1)
+            ]
+            sizes = [step_samples for x in range(n_step - 1)]
+            slices.append(slice((n_step - 1) * step_samples, ob.n_local_samples, 1))
+            sizes.append(ob.n_local_samples - (n_step - 1) * step_samples)
+
+            for det in ob.local_detectors:
+                for slc, sz in zip(slices, sizes):
+                    np.testing.assert_equal(
+                        np.sum(ob.detdata[defaults.det_data][det, slc]), 1.0 * sz
+                    )
+
+        close_data(data)
+
     def test_compare(self):
         # Create a fake satellite data set for testing
         data = create_satellite_data(self.comm)