From 0df6a6414ffb8388d47d2436bf0c48d8888af7c5 Mon Sep 17 00:00:00 2001
From: Dan Smith <dan.smith@maxar.com>
Date: Mon, 22 Jan 2024 12:04:02 -0500
Subject: [PATCH 01/13] There is a bug in VCL or MSVC that causes a compiler
 error w/C++20

---
 six/modules/c++/six/include/six/AmplitudeTable.h | 7 ++++++-
 1 file changed, 6 insertions(+), 1 deletion(-)

diff --git a/six/modules/c++/six/include/six/AmplitudeTable.h b/six/modules/c++/six/include/six/AmplitudeTable.h
index 15349047b..4003958ae 100644
--- a/six/modules/c++/six/include/six/AmplitudeTable.h
+++ b/six/modules/c++/six/include/six/AmplitudeTable.h
@@ -181,7 +181,12 @@ struct SIX_SIX_API AMP8I_PHS8I_t final
         // __has_include is part of C++17
         #if __has_include("../../../six.sicd/include/six/sicd/vectorclass/version2/vectorclass.h") || \
             __has_include("six/sicd/vectorclass/version2/vectorclass.h")
-        #define SIX_sicd_has_VCL 1
+            #if _MSC_VER
+            // Compiler error: bug in MSVC or VCL?
+            #define SIX_sicd_has_VCL !CODA_OSS_cpp20 // TODO: enable for C++20
+            #else
+            #define SIX_sicd_has_VCL 1
+            #endif
         #else
         #define SIX_sicd_has_VCL 0
         #endif // __has_include

From 847c063ff802162463ed936b8ff99c8e805b8eb1 Mon Sep 17 00:00:00 2001
From: Dan Smith <dan.smith@maxar.com>
Date: Mon, 22 Jan 2024 14:03:37 -0500
Subject: [PATCH 02/13] wasn't checking execution_policy::par_unseq

---
 six/modules/c++/six.sicd/source/NearestNeighbors.cpp | 7 +++++++
 1 file changed, 7 insertions(+)

diff --git a/six/modules/c++/six.sicd/source/NearestNeighbors.cpp b/six/modules/c++/six.sicd/source/NearestNeighbors.cpp
index c1af23a71..f1ed35ccb 100644
--- a/six/modules/c++/six.sicd/source/NearestNeighbors.cpp
+++ b/six/modules/c++/six.sicd/source/NearestNeighbors.cpp
@@ -127,6 +127,10 @@ void six::sicd::NearestNeighbors::nearest_neighbors(std::span<const zfloat> inpu
 void six::sicd::NearestNeighbors::nearest_neighbors(execution_policy policy, std::span<const zfloat> inputs, std::span<AMP8I_PHS8I_t> results) const
 {
     #if SIX_sicd_ComplexToAMP8IPHS8I_unseq    
+    if (policy == execution_policy::par_unseq)
+    {
+        return nearest_neighbors_par_unseq(inputs, results);
+    }
     if (policy == execution_policy::unseq)
     {
         return nearest_neighbors_unseq(inputs, results);
@@ -145,6 +149,9 @@ void six::sicd::NearestNeighbors::nearest_neighbors(execution_policy policy, std
     // https://en.cppreference.com/w/cpp/algorithm/execution_policy_tag_t
     // > If the implementation cannot parallelize or vectorize (e.g. due to lack of resources),
     // > all standard execution policies can fall back to sequential execution.
+    #if CODA_OSS_DEBUG
+    throw std::logic_error("Unhandled execution_policy value.");
+    #endif
     return nearest_neighbors(inputs, results); // no policy specified, "default policy"
 }
 

From 47951454e3003c288198d90fec37af3afab5be84 Mon Sep 17 00:00:00 2001
From: "J. Daniel Smith" <dan@jdanielsmith.me>
Date: Tue, 23 Jan 2024 09:16:54 -0500
Subject: [PATCH 03/13] build with C++17

---
 externals/coda-oss/modules/c++/coda-oss.vcxproj               | 1 -
 externals/coda-oss/modules/c++/coda-oss.vcxproj.filters       | 1 +
 six/modules/c++/samples/8AMPI_PHSI.dir/8AMPI_PHSI.dir.vcxproj | 2 ++
 3 files changed, 3 insertions(+), 1 deletion(-)

diff --git a/externals/coda-oss/modules/c++/coda-oss.vcxproj b/externals/coda-oss/modules/c++/coda-oss.vcxproj
index 4d2a9595d..00c993e28 100644
--- a/externals/coda-oss/modules/c++/coda-oss.vcxproj
+++ b/externals/coda-oss/modules/c++/coda-oss.vcxproj
@@ -560,7 +560,6 @@
   <PropertyGroup Label="UserMacros" />
   <PropertyGroup Condition="'$(Configuration)|$(Platform)'=='Debug|x64'">
     <LinkIncremental>true</LinkIncremental>
-    <RunCodeAnalysis>true</RunCodeAnalysis>
     <CodeAnalysisRuleSet>coda-oss.ruleset</CodeAnalysisRuleSet>
   </PropertyGroup>
   <PropertyGroup Condition="'$(Configuration)|$(Platform)'=='Release|x64'">
diff --git a/externals/coda-oss/modules/c++/coda-oss.vcxproj.filters b/externals/coda-oss/modules/c++/coda-oss.vcxproj.filters
index ce620b81e..2bba9d2ea 100644
--- a/externals/coda-oss/modules/c++/coda-oss.vcxproj.filters
+++ b/externals/coda-oss/modules/c++/coda-oss.vcxproj.filters
@@ -954,6 +954,7 @@
     </ClInclude>
     <ClInclude Include="coda_oss\include\coda_oss\mdspan.h" />
     <ClInclude Include="coda_oss\include\coda_oss\mdspan_.h" />
+    <ClInclude Include="sys\include\sys\Ximd.h" />
   </ItemGroup>
   <ItemGroup>
     <ClCompile Include="pch.cpp" />
diff --git a/six/modules/c++/samples/8AMPI_PHSI.dir/8AMPI_PHSI.dir.vcxproj b/six/modules/c++/samples/8AMPI_PHSI.dir/8AMPI_PHSI.dir.vcxproj
index b20b04052..4f968db4f 100644
--- a/six/modules/c++/samples/8AMPI_PHSI.dir/8AMPI_PHSI.dir.vcxproj
+++ b/six/modules/c++/samples/8AMPI_PHSI.dir/8AMPI_PHSI.dir.vcxproj
@@ -66,6 +66,7 @@
       <TreatWarningAsError>true</TreatWarningAsError>
       <RuntimeLibrary>MultiThreadedDebugDLL</RuntimeLibrary>
       <RuntimeTypeInfo>true</RuntimeTypeInfo>
+      <LanguageStandard>stdcpp17</LanguageStandard>
     </ClCompile>
     <Link>
       <SubSystem>Console</SubSystem>
@@ -91,6 +92,7 @@
       <FavorSizeOrSpeed>Speed</FavorSizeOrSpeed>
       <ControlFlowGuard>Guard</ControlFlowGuard>
       <RuntimeTypeInfo>true</RuntimeTypeInfo>
+      <LanguageStandard>stdcpp17</LanguageStandard>
     </ClCompile>
     <Link>
       <SubSystem>Console</SubSystem>

From 95ac5c38df9e02cb224cb1233df309a3eab30d87 Mon Sep 17 00:00:00 2001
From: "J. Daniel Smith" <dan@jdanielsmith.me>
Date: Tue, 23 Jan 2024 09:17:04 -0500
Subject: [PATCH 04/13] iterate over a mdspan

---
 .../include/six/sicd/NearestNeighbors.h       |   4 +-
 .../c++/six.sicd/source/NearestNeighbors.cpp  |   3 +-
 .../source/NearestNeighbors_unseq.cpp         | 168 +++++++++++++-----
 3 files changed, 130 insertions(+), 45 deletions(-)

diff --git a/six/modules/c++/six.sicd/include/six/sicd/NearestNeighbors.h b/six/modules/c++/six.sicd/include/six/sicd/NearestNeighbors.h
index ea5aaf159..6e53c3b81 100644
--- a/six/modules/c++/six.sicd/include/six/sicd/NearestNeighbors.h
+++ b/six/modules/c++/six.sicd/include/six/sicd/NearestNeighbors.h
@@ -69,8 +69,8 @@ struct NearestNeighbors final
     uint8_t getPhase(six::zfloat) const;
 
 #if SIX_sicd_ComplexToAMP8IPHS8I_unseq
-    template<typename ZFloatV, size_t N>
-    auto nearest_neighbors_unseq_T(const std::array<const zfloat, N>&) const; // TODO: std::span<T, N> ... ?
+    template<typename ZFloatV>
+    auto nearest_neighbors_unseq_T(std::span<const zfloat> p) const; // TODO: std::span<T, N> ... ?
     template<typename ZFloatV, int elements_per_iteration>
     void nearest_neighbors_unseq_(std::span<const zfloat> inputs, std::span<AMP8I_PHS8I_t> results) const;
 
diff --git a/six/modules/c++/six.sicd/source/NearestNeighbors.cpp b/six/modules/c++/six.sicd/source/NearestNeighbors.cpp
index f1ed35ccb..71f9fe62b 100644
--- a/six/modules/c++/six.sicd/source/NearestNeighbors.cpp
+++ b/six/modules/c++/six.sicd/source/NearestNeighbors.cpp
@@ -151,8 +151,9 @@ void six::sicd::NearestNeighbors::nearest_neighbors(execution_policy policy, std
     // > all standard execution policies can fall back to sequential execution.
     #if CODA_OSS_DEBUG
     throw std::logic_error("Unhandled execution_policy value.");
-    #endif
+    #else
     return nearest_neighbors(inputs, results); // no policy specified, "default policy"
+    #endif
 }
 
 AMP8I_PHS8I six::sicd::nearest_neighbor(const details::ComplexToAMP8IPHS8I& converter_, const zfloat& v)
diff --git a/six/modules/c++/six.sicd/source/NearestNeighbors_unseq.cpp b/six/modules/c++/six.sicd/source/NearestNeighbors_unseq.cpp
index 1c57eeb47..d34032c9e 100644
--- a/six/modules/c++/six.sicd/source/NearestNeighbors_unseq.cpp
+++ b/six/modules/c++/six.sicd/source/NearestNeighbors_unseq.cpp
@@ -29,6 +29,7 @@
 #include <algorithm>
 #include <functional>
 #include <std/numbers>
+#include <std/mdspan>
 
 #include <gsl/gsl.h>
 #include <sys/Span.h>
@@ -786,6 +787,8 @@ static auto lookup_and_find_nearest(const six::sicd::details::ComplexToAMP8IPHS8
 
 #if SIX_sicd_ComplexToAMP8IPHS8I_unseq
 
+/******************************************************************************************************/
+
 template<typename IntV>
 struct AMP8I_PHS8I_unseq final
 {
@@ -793,56 +796,123 @@ struct AMP8I_PHS8I_unseq final
     IntV phase;
 };
 
-// Cast a blob of T* into "chunks" of T[N]; "undecay" an array.
-// This isn't quite kosher, see: https://stackoverflow.com/questions/47924103/pointer-interconvertibility-vs-having-the-same-address
-template<size_t N, typename T>
-static inline auto array_cast(std::span<const T> data)
+// This isn't someplace more widely available because:
+// 1. there is no standard iteration for `mdspan` (not even in C++23), and
+// 2. the particular needs are quite varied and require a lot of parameterization.
+//
+// Keep things simple for now and do this "in place."
+template<typename T, typename TExtents, typename TValueType = std::span<T>>
+struct mdspan_iterator final
+{
+    using difference_type = std::ptrdiff_t;
+    using size_type = size_t;
+    using value_type = TValueType;
+    using pointer = TValueType*;
+    using const_pointer = const TValueType*;
+    using reference = TValueType&;
+    using const_reference = const TValueType&;
+    using iterator_category = std::random_access_iterator_tag;
+
+    mdspan_iterator() = default;
+    mdspan_iterator(coda_oss::mdspan<T, TExtents> md, difference_type row = 0) : md_(md), row_(row) {}
+
+    value_type operator*() const
+    {
+        return make_span();
+    }
+
+    mdspan_iterator& operator++() noexcept
+    {
+        ++row_;
+        return *this;
+    }
+
+    mdspan_iterator& operator+=(const difference_type n) noexcept
+    {
+        row_ += n;
+        return *this;
+    }
+    mdspan_iterator operator+(const difference_type n) const noexcept
+    {
+        auto ret = *this;
+        ret += n;
+        return ret;
+    }
+
+    difference_type operator-(const mdspan_iterator& rhs) const noexcept
+    {
+        return row_ - rhs.row_;
+    }
+
+    bool operator!=(const mdspan_iterator& rhs) const noexcept
+    {
+        return row_ != rhs.row_;
+    }
+
+private:
+    difference_type row_ = 0;
+    coda_oss::mdspan<T, TExtents> md_;
+
+    auto make_span() const
+    {
+        // We know our mdspan is contiguous, so this is OK
+        auto&& v = md_(row_, 0); // beginning of the the current row
+        return std::span<T>(&v, md_.extent(1)); // span for the whole row
+    }
+};
+template<typename TValueType, typename T, typename TExtents>
+auto begin(coda_oss::mdspan<T, TExtents> md)
+{
+    return mdspan_iterator<T, TExtents, TValueType>(md, 0);
+}
+template<typename T, typename TExtents>
+auto end(coda_oss::mdspan<T, TExtents> md)
+{
+    return mdspan_iterator<T, TExtents>(md, md.extent(0));
+}
+
+template<typename T, typename TExtents>
+using const_mdspan_iterator = mdspan_iterator<std::add_const_t<T>, TExtents>;
+template<typename T, typename TExtents>
+auto cbegin(coda_oss::mdspan<T, TExtents> md)
+{
+    return const_mdspan_iterator<T, TExtents>(md, 0);
+}
+template<typename T, typename TExtents>
+auto cend(coda_oss::mdspan<T, TExtents> md)
 {
-    using chunk_t = std::array<const T, N>;
-    const auto number_of_chunks = data.size() / N;
-    const void* const pData = data.data();
-    return sys::make_span(static_cast<const chunk_t*>(pData), number_of_chunks);
+    return const_mdspan_iterator<T, TExtents>(md, md.extent(0));
 }
 
-// Inside of `std::transform()`, there is a copy; something like
+// Inside of `std::transform()`, there is an assignment; something like
 // ```
 //    *dest = func(*first);
 // ```
-// By returning our own class from `func()`, we can take control of the assignment operator
-// and use that to copy from `AMP8I_PHS8I_unseq` to `std::array<AMP8I_PHS8I, N>&`.
+// By returning our own class from `func()`, we can take control of the assignment operator.
 // (Unlike most other operators, `operator=()` *must* be a member-function.)
-//
-// Using inheritance to avoid padding at the end of the `struct`. ... needed?
-template<typename IntV, size_t N>
-struct AMP8I_PHS8I_array final : public std::array<AMP8I_PHS8I, N>
+template<typename IntV>
+struct mdspan_iterator_value final
 {
-    AMP8I_PHS8I_array& operator=(const AMP8I_PHS8I_unseq<IntV>&) = delete; // should only be using move-assignment
-    AMP8I_PHS8I_array& operator=(AMP8I_PHS8I_unseq<IntV>&& other)
-    {
-        for (int i = 0; i < gsl::narrow<int>(N); i++)
+    std::span<AMP8I_PHS8I> p_;
+
+    mdspan_iterator_value(std::span<AMP8I_PHS8I> s) : p_(s) {}
+    mdspan_iterator_value& operator=(const AMP8I_PHS8I_unseq<IntV>& other) {
+        //assert(p_.size() <= size(other.amplitude));
+        for (int i = 0; i < p_.size(); i++)
         {
-            (*this)[i].phase = gsl::narrow<uint8_t>(other.phase[i]);
-            (*this)[i].amplitude = gsl::narrow<uint8_t>(other.amplitude[i]);
+            p_[i].amplitude = gsl::narrow<uint8_t>(other.amplitude[i]);
+            p_[i].phase = gsl::narrow<uint8_t>(other.phase[i]);
         }
         return *this;
     }
 };
-template<typename IntV, size_t N>
-static inline auto AMP8I_PHS8I_array_cast(std::span<AMP8I_PHS8I> data)
-{
-    // This isn't quite kosher, see `array_cast()` above.
-    using chunk_t = AMP8I_PHS8I_array<IntV, N>;
-    const auto number_of_chunks = data.size() / N;
-    void* const pDest = data.data();
-    return std::span<chunk_t>(static_cast<chunk_t*>(pDest), number_of_chunks);
-}
 
 template<typename ZFloatV>
 using IntV = decltype(::getPhase(ZFloatV{}, 0.0f));
 
 // The compiler can sometimes do better optimization with fixed-size structures.
-template<typename ZFloatV, size_t N>
-auto six::sicd::NearestNeighbors::nearest_neighbors_unseq_T(const std::array<const zfloat, N>& p) const // TODO: std::span<T, N> ... ?
+template<typename ZFloatV>
+auto six::sicd::NearestNeighbors::nearest_neighbors_unseq_T(std::span<const zfloat> p) const // TODO: std::span<T, N> ... ?
 {
     ZFloatV v;
     assert(p.size() == size(v));
@@ -906,16 +976,22 @@ void six::sicd::NearestNeighbors::nearest_neighbors_unseq_(std::span<const zfloa
     // View the data as chunks of *elements_per_iteration*.  This allows iterating
     // to go *elements_per_iteration* at a time; and each chunk can be processed
     // using `nearest_neighbors_unseq_T()`, above.
-    //
-    // This isn't quite kosher, see: https://stackoverflow.com/questions/47924103/pointer-interconvertibility-vs-having-the-same-address
-    auto const first = array_cast<elements_per_iteration>(inputs);
-    auto const dest = AMP8I_PHS8I_array_cast<intv_t, elements_per_iteration>(results);
+    using extents_t = coda_oss::dextents<size_t, 2>; // two dimensions: M×N
+    const extents_t extents{ inputs.size() / elements_per_iteration, elements_per_iteration };
+    const coda_oss::mdspan<const zfloat, extents_t> md_inputs(inputs.data(), extents);
+    assert(md_inputs.size() <= inputs.size());
+    auto const b = cbegin(md_inputs);
+    auto const e = cend(md_inputs);
+
+    const coda_oss::mdspan<AMP8I_PHS8I, extents_t> md_results(results.data(), extents);
+    assert(md_results.size() <= results.size());
+    auto const d = begin<mdspan_iterator_value<intv_t>>(md_results);
 
     const auto func = [&](const auto& v)
     {
         return nearest_neighbors_unseq_T<ZFloatV>(v);
     };
-    std::transform(first.begin(), first.end(), dest.begin(), func);
+    std::transform(/*std::execution::unseq,*/ b, e, d, func);
 
     // Then finish off anything left
     finish_nearest_neighbors_unseq<elements_per_iteration>(*this, inputs, results);
@@ -929,16 +1005,24 @@ void six::sicd::NearestNeighbors::nearest_neighbors_par_unseq_T(std::span<const
     // View the data as chunks of *elements_per_iteration*.  This allows iterating
     // to go *elements_per_iteration* at a time; and each chunk can be processed
     // using `nearest_neighbors_unseq_T()`, above.
-    //
-    // This isn't quite kosher, see: https://stackoverflow.com/questions/47924103/pointer-interconvertibility-vs-having-the-same-address
-    auto const first = array_cast<elements_per_iteration>(inputs);
-    auto const dest = AMP8I_PHS8I_array_cast<intv_t, elements_per_iteration>(results);
+    using extents_t = coda_oss::dextents<size_t, 2>; // two dimensions: M×N
+    const extents_t extents{ inputs.size() / elements_per_iteration, elements_per_iteration };
+    const coda_oss::mdspan<const zfloat, extents_t> md_inputs(inputs.data(), extents);
+    assert(md_inputs.size() <= inputs.size());
+    auto const b = cbegin(md_inputs);
+    auto const e = cend(md_inputs);
+
+    const coda_oss::mdspan<AMP8I_PHS8I, extents_t> md_results(results.data(), extents);
+    assert(md_results.size() <= results.size());
+    auto const d = begin<mdspan_iterator_value<intv_t>>(md_results);
 
     const auto func = [&](const auto& v)
     {
         return nearest_neighbors_unseq_T<ZFloatV>(v);
     };
-    mt::Transform_par(first.begin(), first.end(), dest.begin(), func);
+    //std::transform(std::execution::par_unseq, b, e, d, func);
+    mt::Transform_par(b, e, d, func);
+
 
     // Then finish off anything left
     finish_nearest_neighbors_unseq<elements_per_iteration>(*this, inputs, results);

From 4d1b599984bba36096e56fd7e7296794069959b0 Mon Sep 17 00:00:00 2001
From: "J. Daniel Smith" <dan@jdanielsmith.me>
Date: Tue, 23 Jan 2024 09:17:04 -0500
Subject: [PATCH 05/13] iterate over a mdspan

---
 .../include/six/sicd/NearestNeighbors.h       |   4 +-
 .../c++/six.sicd/source/NearestNeighbors.cpp  |   3 +-
 .../source/NearestNeighbors_unseq.cpp         | 168 +++++++++++++-----
 3 files changed, 130 insertions(+), 45 deletions(-)

diff --git a/six/modules/c++/six.sicd/include/six/sicd/NearestNeighbors.h b/six/modules/c++/six.sicd/include/six/sicd/NearestNeighbors.h
index ea5aaf159..6e53c3b81 100644
--- a/six/modules/c++/six.sicd/include/six/sicd/NearestNeighbors.h
+++ b/six/modules/c++/six.sicd/include/six/sicd/NearestNeighbors.h
@@ -69,8 +69,8 @@ struct NearestNeighbors final
     uint8_t getPhase(six::zfloat) const;
 
 #if SIX_sicd_ComplexToAMP8IPHS8I_unseq
-    template<typename ZFloatV, size_t N>
-    auto nearest_neighbors_unseq_T(const std::array<const zfloat, N>&) const; // TODO: std::span<T, N> ... ?
+    template<typename ZFloatV>
+    auto nearest_neighbors_unseq_T(std::span<const zfloat> p) const; // TODO: std::span<T, N> ... ?
     template<typename ZFloatV, int elements_per_iteration>
     void nearest_neighbors_unseq_(std::span<const zfloat> inputs, std::span<AMP8I_PHS8I_t> results) const;
 
diff --git a/six/modules/c++/six.sicd/source/NearestNeighbors.cpp b/six/modules/c++/six.sicd/source/NearestNeighbors.cpp
index f1ed35ccb..71f9fe62b 100644
--- a/six/modules/c++/six.sicd/source/NearestNeighbors.cpp
+++ b/six/modules/c++/six.sicd/source/NearestNeighbors.cpp
@@ -151,8 +151,9 @@ void six::sicd::NearestNeighbors::nearest_neighbors(execution_policy policy, std
     // > all standard execution policies can fall back to sequential execution.
     #if CODA_OSS_DEBUG
     throw std::logic_error("Unhandled execution_policy value.");
-    #endif
+    #else
     return nearest_neighbors(inputs, results); // no policy specified, "default policy"
+    #endif
 }
 
 AMP8I_PHS8I six::sicd::nearest_neighbor(const details::ComplexToAMP8IPHS8I& converter_, const zfloat& v)
diff --git a/six/modules/c++/six.sicd/source/NearestNeighbors_unseq.cpp b/six/modules/c++/six.sicd/source/NearestNeighbors_unseq.cpp
index 1c57eeb47..d34032c9e 100644
--- a/six/modules/c++/six.sicd/source/NearestNeighbors_unseq.cpp
+++ b/six/modules/c++/six.sicd/source/NearestNeighbors_unseq.cpp
@@ -29,6 +29,7 @@
 #include <algorithm>
 #include <functional>
 #include <std/numbers>
+#include <std/mdspan>
 
 #include <gsl/gsl.h>
 #include <sys/Span.h>
@@ -786,6 +787,8 @@ static auto lookup_and_find_nearest(const six::sicd::details::ComplexToAMP8IPHS8
 
 #if SIX_sicd_ComplexToAMP8IPHS8I_unseq
 
+/******************************************************************************************************/
+
 template<typename IntV>
 struct AMP8I_PHS8I_unseq final
 {
@@ -793,56 +796,123 @@ struct AMP8I_PHS8I_unseq final
     IntV phase;
 };
 
-// Cast a blob of T* into "chunks" of T[N]; "undecay" an array.
-// This isn't quite kosher, see: https://stackoverflow.com/questions/47924103/pointer-interconvertibility-vs-having-the-same-address
-template<size_t N, typename T>
-static inline auto array_cast(std::span<const T> data)
+// This isn't someplace more widely available because:
+// 1. there is no standard iteration for `mdspan` (not even in C++23), and
+// 2. the particular needs are quite varied and require a lot of parameterization.
+//
+// Keep things simple for now and do this "in place."
+template<typename T, typename TExtents, typename TValueType = std::span<T>>
+struct mdspan_iterator final
+{
+    using difference_type = std::ptrdiff_t;
+    using size_type = size_t;
+    using value_type = TValueType;
+    using pointer = TValueType*;
+    using const_pointer = const TValueType*;
+    using reference = TValueType&;
+    using const_reference = const TValueType&;
+    using iterator_category = std::random_access_iterator_tag;
+
+    mdspan_iterator() = default;
+    mdspan_iterator(coda_oss::mdspan<T, TExtents> md, difference_type row = 0) : md_(md), row_(row) {}
+
+    value_type operator*() const
+    {
+        return make_span();
+    }
+
+    mdspan_iterator& operator++() noexcept
+    {
+        ++row_;
+        return *this;
+    }
+
+    mdspan_iterator& operator+=(const difference_type n) noexcept
+    {
+        row_ += n;
+        return *this;
+    }
+    mdspan_iterator operator+(const difference_type n) const noexcept
+    {
+        auto ret = *this;
+        ret += n;
+        return ret;
+    }
+
+    difference_type operator-(const mdspan_iterator& rhs) const noexcept
+    {
+        return row_ - rhs.row_;
+    }
+
+    bool operator!=(const mdspan_iterator& rhs) const noexcept
+    {
+        return row_ != rhs.row_;
+    }
+
+private:
+    difference_type row_ = 0;
+    coda_oss::mdspan<T, TExtents> md_;
+
+    auto make_span() const
+    {
+        // We know our mdspan is contiguous, so this is OK
+        auto&& v = md_(row_, 0); // beginning of the the current row
+        return std::span<T>(&v, md_.extent(1)); // span for the whole row
+    }
+};
+template<typename TValueType, typename T, typename TExtents>
+auto begin(coda_oss::mdspan<T, TExtents> md)
+{
+    return mdspan_iterator<T, TExtents, TValueType>(md, 0);
+}
+template<typename T, typename TExtents>
+auto end(coda_oss::mdspan<T, TExtents> md)
+{
+    return mdspan_iterator<T, TExtents>(md, md.extent(0));
+}
+
+template<typename T, typename TExtents>
+using const_mdspan_iterator = mdspan_iterator<std::add_const_t<T>, TExtents>;
+template<typename T, typename TExtents>
+auto cbegin(coda_oss::mdspan<T, TExtents> md)
+{
+    return const_mdspan_iterator<T, TExtents>(md, 0);
+}
+template<typename T, typename TExtents>
+auto cend(coda_oss::mdspan<T, TExtents> md)
 {
-    using chunk_t = std::array<const T, N>;
-    const auto number_of_chunks = data.size() / N;
-    const void* const pData = data.data();
-    return sys::make_span(static_cast<const chunk_t*>(pData), number_of_chunks);
+    return const_mdspan_iterator<T, TExtents>(md, md.extent(0));
 }
 
-// Inside of `std::transform()`, there is a copy; something like
+// Inside of `std::transform()`, there is an assignment; something like
 // ```
 //    *dest = func(*first);
 // ```
-// By returning our own class from `func()`, we can take control of the assignment operator
-// and use that to copy from `AMP8I_PHS8I_unseq` to `std::array<AMP8I_PHS8I, N>&`.
+// By returning our own class from `func()`, we can take control of the assignment operator.
 // (Unlike most other operators, `operator=()` *must* be a member-function.)
-//
-// Using inheritance to avoid padding at the end of the `struct`. ... needed?
-template<typename IntV, size_t N>
-struct AMP8I_PHS8I_array final : public std::array<AMP8I_PHS8I, N>
+template<typename IntV>
+struct mdspan_iterator_value final
 {
-    AMP8I_PHS8I_array& operator=(const AMP8I_PHS8I_unseq<IntV>&) = delete; // should only be using move-assignment
-    AMP8I_PHS8I_array& operator=(AMP8I_PHS8I_unseq<IntV>&& other)
-    {
-        for (int i = 0; i < gsl::narrow<int>(N); i++)
+    std::span<AMP8I_PHS8I> p_;
+
+    mdspan_iterator_value(std::span<AMP8I_PHS8I> s) : p_(s) {}
+    mdspan_iterator_value& operator=(const AMP8I_PHS8I_unseq<IntV>& other) {
+        //assert(p_.size() <= size(other.amplitude));
+        for (int i = 0; i < p_.size(); i++)
         {
-            (*this)[i].phase = gsl::narrow<uint8_t>(other.phase[i]);
-            (*this)[i].amplitude = gsl::narrow<uint8_t>(other.amplitude[i]);
+            p_[i].amplitude = gsl::narrow<uint8_t>(other.amplitude[i]);
+            p_[i].phase = gsl::narrow<uint8_t>(other.phase[i]);
         }
         return *this;
     }
 };
-template<typename IntV, size_t N>
-static inline auto AMP8I_PHS8I_array_cast(std::span<AMP8I_PHS8I> data)
-{
-    // This isn't quite kosher, see `array_cast()` above.
-    using chunk_t = AMP8I_PHS8I_array<IntV, N>;
-    const auto number_of_chunks = data.size() / N;
-    void* const pDest = data.data();
-    return std::span<chunk_t>(static_cast<chunk_t*>(pDest), number_of_chunks);
-}
 
 template<typename ZFloatV>
 using IntV = decltype(::getPhase(ZFloatV{}, 0.0f));
 
 // The compiler can sometimes do better optimization with fixed-size structures.
-template<typename ZFloatV, size_t N>
-auto six::sicd::NearestNeighbors::nearest_neighbors_unseq_T(const std::array<const zfloat, N>& p) const // TODO: std::span<T, N> ... ?
+template<typename ZFloatV>
+auto six::sicd::NearestNeighbors::nearest_neighbors_unseq_T(std::span<const zfloat> p) const // TODO: std::span<T, N> ... ?
 {
     ZFloatV v;
     assert(p.size() == size(v));
@@ -906,16 +976,22 @@ void six::sicd::NearestNeighbors::nearest_neighbors_unseq_(std::span<const zfloa
     // View the data as chunks of *elements_per_iteration*.  This allows iterating
     // to go *elements_per_iteration* at a time; and each chunk can be processed
     // using `nearest_neighbors_unseq_T()`, above.
-    //
-    // This isn't quite kosher, see: https://stackoverflow.com/questions/47924103/pointer-interconvertibility-vs-having-the-same-address
-    auto const first = array_cast<elements_per_iteration>(inputs);
-    auto const dest = AMP8I_PHS8I_array_cast<intv_t, elements_per_iteration>(results);
+    using extents_t = coda_oss::dextents<size_t, 2>; // two dimensions: M×N
+    const extents_t extents{ inputs.size() / elements_per_iteration, elements_per_iteration };
+    const coda_oss::mdspan<const zfloat, extents_t> md_inputs(inputs.data(), extents);
+    assert(md_inputs.size() <= inputs.size());
+    auto const b = cbegin(md_inputs);
+    auto const e = cend(md_inputs);
+
+    const coda_oss::mdspan<AMP8I_PHS8I, extents_t> md_results(results.data(), extents);
+    assert(md_results.size() <= results.size());
+    auto const d = begin<mdspan_iterator_value<intv_t>>(md_results);
 
     const auto func = [&](const auto& v)
     {
         return nearest_neighbors_unseq_T<ZFloatV>(v);
     };
-    std::transform(first.begin(), first.end(), dest.begin(), func);
+    std::transform(/*std::execution::unseq,*/ b, e, d, func);
 
     // Then finish off anything left
     finish_nearest_neighbors_unseq<elements_per_iteration>(*this, inputs, results);
@@ -929,16 +1005,24 @@ void six::sicd::NearestNeighbors::nearest_neighbors_par_unseq_T(std::span<const
     // View the data as chunks of *elements_per_iteration*.  This allows iterating
     // to go *elements_per_iteration* at a time; and each chunk can be processed
     // using `nearest_neighbors_unseq_T()`, above.
-    //
-    // This isn't quite kosher, see: https://stackoverflow.com/questions/47924103/pointer-interconvertibility-vs-having-the-same-address
-    auto const first = array_cast<elements_per_iteration>(inputs);
-    auto const dest = AMP8I_PHS8I_array_cast<intv_t, elements_per_iteration>(results);
+    using extents_t = coda_oss::dextents<size_t, 2>; // two dimensions: M×N
+    const extents_t extents{ inputs.size() / elements_per_iteration, elements_per_iteration };
+    const coda_oss::mdspan<const zfloat, extents_t> md_inputs(inputs.data(), extents);
+    assert(md_inputs.size() <= inputs.size());
+    auto const b = cbegin(md_inputs);
+    auto const e = cend(md_inputs);
+
+    const coda_oss::mdspan<AMP8I_PHS8I, extents_t> md_results(results.data(), extents);
+    assert(md_results.size() <= results.size());
+    auto const d = begin<mdspan_iterator_value<intv_t>>(md_results);
 
     const auto func = [&](const auto& v)
     {
         return nearest_neighbors_unseq_T<ZFloatV>(v);
     };
-    mt::Transform_par(first.begin(), first.end(), dest.begin(), func);
+    //std::transform(std::execution::par_unseq, b, e, d, func);
+    mt::Transform_par(b, e, d, func);
+
 
     // Then finish off anything left
     finish_nearest_neighbors_unseq<elements_per_iteration>(*this, inputs, results);

From abd5cc7f5a96edca50adc7214753a58e3133fdaf Mon Sep 17 00:00:00 2001
From: Dan Smith <dan.smith@maxar.com>
Date: Tue, 23 Jan 2024 10:59:36 -0500
Subject: [PATCH 06/13] fix 'reorder' warning

---
 six/modules/c++/six.sicd/source/NearestNeighbors_unseq.cpp | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/six/modules/c++/six.sicd/source/NearestNeighbors_unseq.cpp b/six/modules/c++/six.sicd/source/NearestNeighbors_unseq.cpp
index d34032c9e..bd88056a9 100644
--- a/six/modules/c++/six.sicd/source/NearestNeighbors_unseq.cpp
+++ b/six/modules/c++/six.sicd/source/NearestNeighbors_unseq.cpp
@@ -850,8 +850,8 @@ struct mdspan_iterator final
     }
 
 private:
-    difference_type row_ = 0;
     coda_oss::mdspan<T, TExtents> md_;
+    difference_type row_ = 0;
 
     auto make_span() const
     {

From d87f0221bb7201f4adf1529b96f85d9b2d20a004 Mon Sep 17 00:00:00 2001
From: Dan Smith <dan.smith@maxar.com>
Date: Tue, 23 Jan 2024 10:59:36 -0500
Subject: [PATCH 07/13] fix 'reorder' warning

---
 six/modules/c++/six.sicd/source/NearestNeighbors_unseq.cpp | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/six/modules/c++/six.sicd/source/NearestNeighbors_unseq.cpp b/six/modules/c++/six.sicd/source/NearestNeighbors_unseq.cpp
index d34032c9e..bd88056a9 100644
--- a/six/modules/c++/six.sicd/source/NearestNeighbors_unseq.cpp
+++ b/six/modules/c++/six.sicd/source/NearestNeighbors_unseq.cpp
@@ -850,8 +850,8 @@ struct mdspan_iterator final
     }
 
 private:
-    difference_type row_ = 0;
     coda_oss::mdspan<T, TExtents> md_;
+    difference_type row_ = 0;
 
     auto make_span() const
     {

From 14c5fe7a29016a9618e423c7259ab86483140bf8 Mon Sep 17 00:00:00 2001
From: Dan Smith <dan.smith@maxar.com>
Date: Tue, 23 Jan 2024 11:24:26 -0500
Subject: [PATCH 08/13] fix 'reorder' warning

---
 six/modules/c++/six.sicd/source/NearestNeighbors_unseq.cpp | 7 ++++---
 1 file changed, 4 insertions(+), 3 deletions(-)

diff --git a/six/modules/c++/six.sicd/source/NearestNeighbors_unseq.cpp b/six/modules/c++/six.sicd/source/NearestNeighbors_unseq.cpp
index bd88056a9..dde1dbb02 100644
--- a/six/modules/c++/six.sicd/source/NearestNeighbors_unseq.cpp
+++ b/six/modules/c++/six.sicd/source/NearestNeighbors_unseq.cpp
@@ -898,10 +898,11 @@ struct mdspan_iterator_value final
     mdspan_iterator_value(std::span<AMP8I_PHS8I> s) : p_(s) {}
     mdspan_iterator_value& operator=(const AMP8I_PHS8I_unseq<IntV>& other) {
         //assert(p_.size() <= size(other.amplitude));
-        for (int i = 0; i < p_.size(); i++)
+        for (size_t i = 0; i < p_.size(); i++)
         {
-            p_[i].amplitude = gsl::narrow<uint8_t>(other.amplitude[i]);
-            p_[i].phase = gsl::narrow<uint8_t>(other.phase[i]);
+	    const auto i_ = gsl::narrow<ptrdiff_t>(i);
+            p_[i].amplitude = gsl::narrow<uint8_t>(other.amplitude[i_]);
+            p_[i].phase = gsl::narrow<uint8_t>(other.phase[i_]);
         }
         return *this;
     }

From a0265d595fbbb6dbf7c22da705e307f35fdbe5fe Mon Sep 17 00:00:00 2001
From: Dan Smith <dan.smith@maxar.com>
Date: Tue, 23 Jan 2024 11:24:26 -0500
Subject: [PATCH 09/13] fix 'reorder' warning

---
 six/modules/c++/six.sicd/source/NearestNeighbors_unseq.cpp | 7 ++++---
 1 file changed, 4 insertions(+), 3 deletions(-)

diff --git a/six/modules/c++/six.sicd/source/NearestNeighbors_unseq.cpp b/six/modules/c++/six.sicd/source/NearestNeighbors_unseq.cpp
index bd88056a9..023eddd76 100644
--- a/six/modules/c++/six.sicd/source/NearestNeighbors_unseq.cpp
+++ b/six/modules/c++/six.sicd/source/NearestNeighbors_unseq.cpp
@@ -898,10 +898,11 @@ struct mdspan_iterator_value final
     mdspan_iterator_value(std::span<AMP8I_PHS8I> s) : p_(s) {}
     mdspan_iterator_value& operator=(const AMP8I_PHS8I_unseq<IntV>& other) {
         //assert(p_.size() <= size(other.amplitude));
-        for (int i = 0; i < p_.size(); i++)
+        for (size_t i = 0; i < p_.size(); i++)
         {
-            p_[i].amplitude = gsl::narrow<uint8_t>(other.amplitude[i]);
-            p_[i].phase = gsl::narrow<uint8_t>(other.phase[i]);
+	    const auto i_ = gsl::narrow<int>(i);
+            p_[i].amplitude = gsl::narrow<uint8_t>(other.amplitude[i_]);
+            p_[i].phase = gsl::narrow<uint8_t>(other.phase[i_]);
         }
         return *this;
     }

From f9a6e1b271a3c0f2f0fce3e6f7cf232ca69cb6c7 Mon Sep 17 00:00:00 2001
From: Dan Smith <dan.smith@maxar.com>
Date: Tue, 23 Jan 2024 12:43:39 -0500
Subject: [PATCH 10/13] try making just operator=() a template

---
 .../c++/six.sicd/source/NearestNeighbors_unseq.cpp        | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/six/modules/c++/six.sicd/source/NearestNeighbors_unseq.cpp b/six/modules/c++/six.sicd/source/NearestNeighbors_unseq.cpp
index 023eddd76..08002db8a 100644
--- a/six/modules/c++/six.sicd/source/NearestNeighbors_unseq.cpp
+++ b/six/modules/c++/six.sicd/source/NearestNeighbors_unseq.cpp
@@ -890,17 +890,17 @@ auto cend(coda_oss::mdspan<T, TExtents> md)
 // ```
 // By returning our own class from `func()`, we can take control of the assignment operator.
 // (Unlike most other operators, `operator=()` *must* be a member-function.)
-template<typename IntV>
 struct mdspan_iterator_value final
 {
     std::span<AMP8I_PHS8I> p_;
 
     mdspan_iterator_value(std::span<AMP8I_PHS8I> s) : p_(s) {}
+    template<typename IntV>
     mdspan_iterator_value& operator=(const AMP8I_PHS8I_unseq<IntV>& other) {
         //assert(p_.size() <= size(other.amplitude));
         for (size_t i = 0; i < p_.size(); i++)
         {
-	    const auto i_ = gsl::narrow<int>(i);
+	        const auto i_ = gsl::narrow<int>(i);
             p_[i].amplitude = gsl::narrow<uint8_t>(other.amplitude[i_]);
             p_[i].phase = gsl::narrow<uint8_t>(other.phase[i_]);
         }
@@ -986,7 +986,7 @@ void six::sicd::NearestNeighbors::nearest_neighbors_unseq_(std::span<const zfloa
 
     const coda_oss::mdspan<AMP8I_PHS8I, extents_t> md_results(results.data(), extents);
     assert(md_results.size() <= results.size());
-    auto const d = begin<mdspan_iterator_value<intv_t>>(md_results);
+    auto const d = begin<mdspan_iterator_value>(md_results);
 
     const auto func = [&](const auto& v)
     {
@@ -1015,7 +1015,7 @@ void six::sicd::NearestNeighbors::nearest_neighbors_par_unseq_T(std::span<const
 
     const coda_oss::mdspan<AMP8I_PHS8I, extents_t> md_results(results.data(), extents);
     assert(md_results.size() <= results.size());
-    auto const d = begin<mdspan_iterator_value<intv_t>>(md_results);
+    auto const d = begin<mdspan_iterator_value>(md_results);
 
     const auto func = [&](const auto& v)
     {

From 4b0565a687dedc482032f03a9f3385d3c9a26e6d Mon Sep 17 00:00:00 2001
From: Dan Smith <dan.smith@maxar.com>
Date: Tue, 23 Jan 2024 12:53:36 -0500
Subject: [PATCH 11/13] trying to get iterator working in more situations

---
 .../c++/six.sicd/source/NearestNeighbors_unseq.cpp        | 8 ++------
 1 file changed, 2 insertions(+), 6 deletions(-)

diff --git a/six/modules/c++/six.sicd/source/NearestNeighbors_unseq.cpp b/six/modules/c++/six.sicd/source/NearestNeighbors_unseq.cpp
index 08002db8a..330ff0d8f 100644
--- a/six/modules/c++/six.sicd/source/NearestNeighbors_unseq.cpp
+++ b/six/modules/c++/six.sicd/source/NearestNeighbors_unseq.cpp
@@ -809,8 +809,8 @@ struct mdspan_iterator final
     using value_type = TValueType;
     using pointer = TValueType*;
     using const_pointer = const TValueType*;
-    using reference = TValueType&;
-    using const_reference = const TValueType&;
+    using reference = TValueType;
+    using const_reference = const TValueType;
     using iterator_category = std::random_access_iterator_tag;
 
     mdspan_iterator() = default;
@@ -972,8 +972,6 @@ static void finish_nearest_neighbors_unseq(const six::sicd::NearestNeighbors& im
 template<typename ZFloatV, int elements_per_iteration>
 void six::sicd::NearestNeighbors::nearest_neighbors_unseq_(std::span<const zfloat> inputs, std::span<AMP8I_PHS8I> results) const
 {
-    using intv_t = IntV<ZFloatV>;
-
     // View the data as chunks of *elements_per_iteration*.  This allows iterating
     // to go *elements_per_iteration* at a time; and each chunk can be processed
     // using `nearest_neighbors_unseq_T()`, above.
@@ -1001,8 +999,6 @@ void six::sicd::NearestNeighbors::nearest_neighbors_unseq_(std::span<const zfloa
 template<typename ZFloatV, int elements_per_iteration>
 void six::sicd::NearestNeighbors::nearest_neighbors_par_unseq_T(std::span<const zfloat> inputs, std::span<AMP8I_PHS8I> results) const
 {
-    using intv_t = IntV<ZFloatV>;
-
     // View the data as chunks of *elements_per_iteration*.  This allows iterating
     // to go *elements_per_iteration* at a time; and each chunk can be processed
     // using `nearest_neighbors_unseq_T()`, above.

From 26f34a4e62c320ed654b2be53f0d88b2ade19218 Mon Sep 17 00:00:00 2001
From: Dan Smith <dan.smith@maxar.com>
Date: Tue, 23 Jan 2024 13:18:31 -0500
Subject: [PATCH 12/13] remove debug code

---
 .../source/NearestNeighbors_unseq.cpp         | 43 +------------------
 1 file changed, 2 insertions(+), 41 deletions(-)

diff --git a/six/modules/c++/six.sicd/source/NearestNeighbors_unseq.cpp b/six/modules/c++/six.sicd/source/NearestNeighbors_unseq.cpp
index 330ff0d8f..564c5415c 100644
--- a/six/modules/c++/six.sicd/source/NearestNeighbors_unseq.cpp
+++ b/six/modules/c++/six.sicd/source/NearestNeighbors_unseq.cpp
@@ -680,18 +680,7 @@ inline auto lower_bound_(std::span<const float> magnitudes, const FloatV& v)
 template<typename IntV, typename FloatV>
 inline auto lower_bound(std::span<const float> magnitudes, const FloatV& value)
 {
-    auto retval = lower_bound_<IntV>(magnitudes, value);
-
-    #if CODA_OSS_DEBUG
-    for (int i = 0; i < ssize(value); i++)
-    {
-        const auto it = std::lower_bound(magnitudes.begin(), magnitudes.end(), value[i]);
-        const auto result = gsl::narrow<int>(std::distance(magnitudes.begin(), it));
-        assert(retval[i] == result);
-    }
-    #endif
-
-    return retval;
+    return lower_bound_<IntV>(magnitudes, value);
 }
 
 template<typename IntV, typename FloatV>
@@ -911,45 +900,17 @@ struct mdspan_iterator_value final
 template<typename ZFloatV>
 using IntV = decltype(::getPhase(ZFloatV{}, 0.0f));
 
-// The compiler can sometimes do better optimization with fixed-size structures.
 template<typename ZFloatV>
-auto six::sicd::NearestNeighbors::nearest_neighbors_unseq_T(std::span<const zfloat> p) const // TODO: std::span<T, N> ... ?
+auto six::sicd::NearestNeighbors::nearest_neighbors_unseq_T(std::span<const zfloat> p) const // TODO: std::span<T, N> ... ?  The compiler can sometimes do better optimization with fixed-size structures.
 {
     ZFloatV v;
     assert(p.size() == size(v));
-
     copy_from(p, v);
-    #if CODA_OSS_DEBUG
-    for (int i = 0; i < ssize(v); i++)
-    {
-        //const auto z = p[i];
-        //assert(real(v)[i] == z.real());
-        //assert(imag(v)[i] == z.imag());
-    }
-    #endif
 
     using intv_t = IntV<ZFloatV>;
     AMP8I_PHS8I_unseq<intv_t> retval;
-
     retval.phase = ::getPhase(v, converter_.phase_delta());
-    #if CODA_OSS_DEBUG
-    for (int i = 0; i < ssize(retval.phase); i++)
-    {
-        const auto phase_ = converter_.getPhase(p[i]);
-        assert(static_cast<uint8_t>(retval.phase[i]) == phase_);
-    }
-    #endif
-
     retval.amplitude = lookup_and_find_nearest(converter_, retval.phase, v);
-    #if CODA_OSS_DEBUG
-    for (int i = 0; i < ssize(retval.amplitude); i++)
-    {
-        const auto i_ = retval.phase[i];
-        const auto& phase_directions = converter_.get_phase_directions().value;
-        const auto a = find_nearest(phase_directions[i_], p[i]);
-        assert(a == retval.amplitude[i]);
-    }
-    #endif
 
     return retval;
 }

From e6ec7010d34e3ac58a931c733e6cf90c68910b63 Mon Sep 17 00:00:00 2001
From: Dan Smith <dan.smith@maxar.com>
Date: Tue, 23 Jan 2024 13:31:29 -0500
Subject: [PATCH 13/13] reduce duplicate code

---
 .../include/six/sicd/NearestNeighbors.h       | 12 +--
 .../source/NearestNeighbors_unseq.cpp         | 98 ++++++-------------
 2 files changed, 36 insertions(+), 74 deletions(-)

diff --git a/six/modules/c++/six.sicd/include/six/sicd/NearestNeighbors.h b/six/modules/c++/six.sicd/include/six/sicd/NearestNeighbors.h
index 6e53c3b81..e9629551e 100644
--- a/six/modules/c++/six.sicd/include/six/sicd/NearestNeighbors.h
+++ b/six/modules/c++/six.sicd/include/six/sicd/NearestNeighbors.h
@@ -68,15 +68,15 @@ struct NearestNeighbors final
     uint8_t find_nearest(six::zfloat phase_direction, six::zfloat v) const;
     uint8_t getPhase(six::zfloat) const;
 
-#if SIX_sicd_ComplexToAMP8IPHS8I_unseq
+    #if SIX_sicd_ComplexToAMP8IPHS8I_unseq
+    void nearest_neighbors_(execution_policy, std::span<const zfloat> inputs, std::span<AMP8I_PHS8I_t> results) const;
+
     template<typename ZFloatV>
-    auto nearest_neighbors_unseq_T(std::span<const zfloat> p) const; // TODO: std::span<T, N> ... ?
-    template<typename ZFloatV, int elements_per_iteration>
-    void nearest_neighbors_unseq_(std::span<const zfloat> inputs, std::span<AMP8I_PHS8I_t> results) const;
+    auto unseq_nearest_neighbors(std::span<const zfloat> p) const; // TODO: std::span<T, N> ... ?
 
     template<typename ZFloatV, int elements_per_iteration>
-    void nearest_neighbors_par_unseq_T(std::span<const zfloat> inputs, std::span<AMP8I_PHS8I_t> results) const;
-#endif 
+    void nearest_neighbors_T(execution_policy, std::span<const zfloat> inputs, std::span<AMP8I_PHS8I_t> results) const;
+    #endif 
 
 };
 }
diff --git a/six/modules/c++/six.sicd/source/NearestNeighbors_unseq.cpp b/six/modules/c++/six.sicd/source/NearestNeighbors_unseq.cpp
index 564c5415c..0de69a6fc 100644
--- a/six/modules/c++/six.sicd/source/NearestNeighbors_unseq.cpp
+++ b/six/modules/c++/six.sicd/source/NearestNeighbors_unseq.cpp
@@ -901,7 +901,7 @@ template<typename ZFloatV>
 using IntV = decltype(::getPhase(ZFloatV{}, 0.0f));
 
 template<typename ZFloatV>
-auto six::sicd::NearestNeighbors::nearest_neighbors_unseq_T(std::span<const zfloat> p) const // TODO: std::span<T, N> ... ?  The compiler can sometimes do better optimization with fixed-size structures.
+auto six::sicd::NearestNeighbors::unseq_nearest_neighbors(std::span<const zfloat> p) const // TODO: std::span<T, N> ... ?  The compiler can sometimes do better optimization with fixed-size structures.
 {
     ZFloatV v;
     assert(p.size() == size(v));
@@ -931,7 +931,8 @@ static void finish_nearest_neighbors_unseq(const six::sicd::NearestNeighbors& im
 }
 
 template<typename ZFloatV, int elements_per_iteration>
-void six::sicd::NearestNeighbors::nearest_neighbors_unseq_(std::span<const zfloat> inputs, std::span<AMP8I_PHS8I> results) const
+void six::sicd::NearestNeighbors::nearest_neighbors_T(execution_policy policy,
+    std::span<const zfloat> inputs, std::span<AMP8I_PHS8I> results) const
 {
     // View the data as chunks of *elements_per_iteration*.  This allows iterating
     // to go *elements_per_iteration* at a time; and each chunk can be processed
@@ -949,38 +950,22 @@ void six::sicd::NearestNeighbors::nearest_neighbors_unseq_(std::span<const zfloa
 
     const auto func = [&](const auto& v)
     {
-        return nearest_neighbors_unseq_T<ZFloatV>(v);
+        return unseq_nearest_neighbors<ZFloatV>(v);
     };
-    std::transform(/*std::execution::unseq,*/ b, e, d, func);
 
-    // Then finish off anything left
-    finish_nearest_neighbors_unseq<elements_per_iteration>(*this, inputs, results);
-}
-
-template<typename ZFloatV, int elements_per_iteration>
-void six::sicd::NearestNeighbors::nearest_neighbors_par_unseq_T(std::span<const zfloat> inputs, std::span<AMP8I_PHS8I> results) const
-{
-    // View the data as chunks of *elements_per_iteration*.  This allows iterating
-    // to go *elements_per_iteration* at a time; and each chunk can be processed
-    // using `nearest_neighbors_unseq_T()`, above.
-    using extents_t = coda_oss::dextents<size_t, 2>; // two dimensions: M×N
-    const extents_t extents{ inputs.size() / elements_per_iteration, elements_per_iteration };
-    const coda_oss::mdspan<const zfloat, extents_t> md_inputs(inputs.data(), extents);
-    assert(md_inputs.size() <= inputs.size());
-    auto const b = cbegin(md_inputs);
-    auto const e = cend(md_inputs);
-
-    const coda_oss::mdspan<AMP8I_PHS8I, extents_t> md_results(results.data(), extents);
-    assert(md_results.size() <= results.size());
-    auto const d = begin<mdspan_iterator_value>(md_results);
-
-    const auto func = [&](const auto& v)
+    if (policy == execution_policy::unseq)
     {
-        return nearest_neighbors_unseq_T<ZFloatV>(v);
-    };
-    //std::transform(std::execution::par_unseq, b, e, d, func);
-    mt::Transform_par(b, e, d, func);
-
+        std::transform(/*std::execution::unseq,*/ b, e, d, func);
+    }
+    else if (policy == execution_policy::par_unseq)
+    {
+        //std::transform(std::execution::par_unseq, b, e, d, func);
+        mt::Transform_par(b, e, d, func);
+    }
+    else
+    {
+        throw std::logic_error("Unsupported execution_policy");
+    }
 
     // Then finish off anything left
     finish_nearest_neighbors_unseq<elements_per_iteration>(*this, inputs, results);
@@ -1019,74 +1004,51 @@ std::string SIX_SICD_API six_sicd_set_nearest_neighbors_unseq(std::string unseq)
     return retval;
 }
 
-void six::sicd::NearestNeighbors::nearest_neighbors_unseq(std::span<const zfloat> inputs, std::span<AMP8I_PHS8I> results) const
+void six::sicd::NearestNeighbors::nearest_neighbors_(execution_policy policy,
+    std::span<const zfloat> inputs, std::span<AMP8I_PHS8I_t> results) const
 {
     // TODO: there could be more complicated logic here to determine which UNSEQ
     // implementation to use.
 
-
     // This is very simple as it's only used for unit-testing
     const auto& unseq = ::nearest_neighbors_unseq_;
     #if SIX_sicd_has_simd
     if (unseq == unseq_simd)
     {
-        return nearest_neighbors_unseq_<simd_zfloatv, simd_elements_per_iteration>(inputs, results);
+        return nearest_neighbors_T<simd_zfloatv, simd_elements_per_iteration>(policy, inputs, results);
     }
     #endif
     #if SIX_sicd_has_VCL
     if (unseq == unseq_vcl)
     {
-        return nearest_neighbors_unseq_<vcl_zfloatv, vcl_elements_per_iteration>(inputs, results);
+        return nearest_neighbors_T<vcl_zfloatv, vcl_elements_per_iteration>(policy, inputs, results);
     }
     #endif
     #if SIX_sicd_has_valarray
     if (unseq == unseq_valarray)
     {
-        return nearest_neighbors_unseq_<valarray_zfloatv, valarray_elements_per_iteration>(inputs, results);
+        return nearest_neighbors_T<valarray_zfloatv, valarray_elements_per_iteration>(policy, inputs, results);
     }
     #endif
     #if SIX_sicd_has_ximd
     if (unseq == unseq_ximd)
     {
-        return nearest_neighbors_unseq_<ximd_zfloatv, ximd_elements_per_iteration>(inputs, results);
+        return nearest_neighbors_T<ximd_zfloatv, ximd_elements_per_iteration>(policy, inputs, results);
     }
     #endif
 
-    throw std::logic_error("Don't know how to implement nearest_neighbors_unseq() for unseq=" + unseq);
+    throw std::logic_error("Don't know how to implement nearest_neighbors_() for unseq=" + unseq);
+}
+void six::sicd::NearestNeighbors::nearest_neighbors_unseq(std::span<const zfloat> inputs, std::span<AMP8I_PHS8I> results) const
+{
+    // TODO: there could be more complicated logic here to determine which UNSEQ
+    // implementation to use.
+    nearest_neighbors_(execution_policy::unseq, inputs, results);
 }
-
 void six::sicd::NearestNeighbors::nearest_neighbors_par_unseq(std::span<const zfloat> inputs, std::span<AMP8I_PHS8I> results) const
 {
     // TODO: there could be more complicated logic here to determine which UNSEQ
     // implementation to use.
-
-    // This is very simple as it's only used for unit-testing
-    const auto& unseq = ::nearest_neighbors_unseq_;
-    #if SIX_sicd_has_simd
-    if (unseq == unseq_simd)
-    {
-        return nearest_neighbors_par_unseq_T<simd_zfloatv, simd_elements_per_iteration>(inputs, results);
-    }
-    #endif
-    #if SIX_sicd_has_VCL
-    if (unseq == unseq_vcl)
-    {
-        return nearest_neighbors_par_unseq_T<vcl_zfloatv, vcl_elements_per_iteration>(inputs, results);
-    }
-    #endif
-    #if SIX_sicd_has_valarray
-    if (unseq == unseq_valarray)
-    {
-        return nearest_neighbors_par_unseq_T<valarray_zfloatv, valarray_elements_per_iteration>(inputs, results);
-    }
-    #endif
-    #if SIX_sicd_has_ximd
-    if (unseq == unseq_ximd)
-    {
-        return nearest_neighbors_par_unseq_T<ximd_zfloatv, ximd_elements_per_iteration>(inputs, results);
-    }
-    #endif
-
-    throw std::logic_error("Don't know how to implement nearest_neighbors_par_unseq() for unseq=" + unseq);
+    nearest_neighbors_(execution_policy::par_unseq, inputs, results);
 }
 #endif // SIX_sicd_ComplexToAMP8IPHS8I_unseq