From f29eaf7cc71e635a5242239e0e099ed6444416eb Mon Sep 17 00:00:00 2001
From: Christopher Mauney <mauneyc@pn2305185.lanl.gov>
Date: Mon, 2 Oct 2023 14:45:33 -0600
Subject: [PATCH 01/22] initial commit

---
 ports-of-call/portable_arrays.hpp | 439 +++++++++---------------------
 1 file changed, 135 insertions(+), 304 deletions(-)

diff --git a/ports-of-call/portable_arrays.hpp b/ports-of-call/portable_arrays.hpp
index 81a2ad6d..394e4ec6 100644
--- a/ports-of-call/portable_arrays.hpp
+++ b/ports-of-call/portable_arrays.hpp
@@ -26,73 +26,122 @@
 
 #include "portability.hpp"
 #include <algorithm>
+#include <array>
 #include <assert.h>
 #include <cstddef> // size_t
 #include <cstring> // memset()
+#include <functional>
+#include <numeric>
+#include <type_traits>
 #include <utility> // swap()
 
+constexpr std::size_t MAXDIM = 6;
+using narr = std::array<std::size_t, MAXDIM>;
+namespace detail {
+
+template <std::size_t Ind, typename NX>
+constexpr void set_value(narr &ndat, NX value) {
+  ndat[Ind] = value;
+}
+
+template <std::size_t Ind, typename NX>
+constexpr void set_values(narr &ndat, NX value) {
+  set_value<Ind>(ndat, value);
+}
+
+template <std::size_t Ind, typename NX, typename... NXs>
+constexpr void set_values(narr &ndat, NX value, NXs... nxs) {
+  set_value<Ind>(ndat, value);
+  set_values<Ind - 1>(ndat, nxs...);
+}
+
+template <std::size_t Ind>
+size_t computeIndex(const narr &nd, const size_t index) {
+  return index;
+}
+template <std::size_t Ind, typename... Tail>
+size_t computeIndex(const narr &nd, const size_t index, const Tail... tail) {
+  return index * nd[Ind] + computeIndex<Ind + 1>(nd, tail...);
+}
+
+template <typename NX>
+PORTABLE_INLINE_FUNCTION auto varmul(NX v) {
+  return v;
+}
+
+template <typename NX, typename... NXs>
+PORTABLE_INLINE_FUNCTION auto varmul(NX v, NXs... vs) {
+  return v * varmul(vs...);
+}
+
+} // namespace detail
+
+template <typename... NXs>
+PORTABLE_INLINE_FUNCTION auto nx_arr(NXs... nxs) {
+  narr r;
+  constexpr std::size_t NP = sizeof...(NXs);
+  detail::set_values<NP - 1>(r, nxs...);
+  for (auto i = NP; i < MAXDIM; ++i)
+    r[i] = 1;
+  return r;
+}
+
+template <typename... Indicies>
+size_t compute_index(const narr &nd, const Indicies... idxs) {
+  return 0 + detail::computeIndex<0>(nd, idxs...);
+}
+
+template <typename... NXs>
+PORTABLE_INLINE_FUNCTION auto nx_mul(NXs... nxs) {
+  return detail::varmul(nxs...);
+}
+
 template <typename T>
 class PortableMDArray {
  public:
-  static constexpr int MAXDIM = 6;
+  PORTABLE_FUNCTION PortableMDArray(T *data, narr nxa, std::size_t rank)
+      : pdata_(data), nxs_(nxa), rank_(rank) {}
+  // rank_{std::distance(nxs_.cbegin(), std::find(nxs_.cbegin(), nxs_.cend(),
+  // 1))} {}
+
+  template <typename... NXs>
+  PORTABLE_FUNCTION PortableMDArray(T *p, NXs... nxs) noexcept
+      : PortableMDArray(p, nx_arr(nxs...), sizeof...(NXs)) {}
 
   // ctors
   // default ctor: simply set null PortableMDArray
   PORTABLE_FUNCTION
-  PortableMDArray() noexcept
-      : pdata_(nullptr), nx1_(0), nx2_(0), nx3_(0), nx4_(0), nx5_(0), nx6_(0) {}
-  PORTABLE_FUNCTION PortableMDArray(T *data, int nx1) noexcept
-      : pdata_(data), nx1_(nx1), nx2_(1), nx3_(1), nx4_(1), nx5_(1), nx6_(1) {}
-  PORTABLE_FUNCTION
-  PortableMDArray(T *data, int nx2, int nx1) noexcept
-      : pdata_(data), nx1_(nx1), nx2_(nx2), nx3_(1), nx4_(1), nx5_(1), nx6_(1) {
-  }
-  PORTABLE_FUNCTION
-  PortableMDArray(T *data, int nx3, int nx2, int nx1) noexcept
-      : pdata_(data), nx1_(nx1), nx2_(nx2), nx3_(nx3), nx4_(1), nx5_(1),
-        nx6_(1) {}
-  PORTABLE_FUNCTION
-  PortableMDArray(T *data, int nx4, int nx3, int nx2, int nx1) noexcept
-      : pdata_(data), nx1_(nx1), nx2_(nx2), nx3_(nx3), nx4_(nx4), nx5_(1),
-        nx6_(1) {}
-  PORTABLE_FUNCTION
-  PortableMDArray(T *data, int nx5, int nx4, int nx3, int nx2, int nx1) noexcept
-      : pdata_(data), nx1_(nx1), nx2_(nx2), nx3_(nx3), nx4_(nx4), nx5_(nx5),
-        nx6_(1) {}
-  PORTABLE_FUNCTION
-  PortableMDArray(T *data, int nx6, int nx5, int nx4, int nx3, int nx2,
-                  int nx1) noexcept
-      : pdata_(data), nx1_(nx1), nx2_(nx2), nx3_(nx3), nx4_(nx4), nx5_(nx5),
-        nx6_(nx6) {}
+  PortableMDArray() noexcept : pdata_(nullptr), nxs_{{0}}, rank_{0} {}
 
   // define copy constructor and overload assignment operator so both do deep
   // copies.
   PortableMDArray(const PortableMDArray<T> &t) noexcept;
   PortableMDArray<T> &operator=(const PortableMDArray<T> &t) noexcept;
 
-  // public functions to allocate/deallocate memory for 1D-5D data
-  PORTABLE_FUNCTION void NewPortableMDArray(T *data, int nx1) noexcept;
-  PORTABLE_FUNCTION void NewPortableMDArray(T *data, int nx2, int nx1) noexcept;
-  PORTABLE_FUNCTION void NewPortableMDArray(T *data, int nx3, int nx2,
-                                            int nx1) noexcept;
-  PORTABLE_FUNCTION void NewPortableMDArray(T *data, int nx4, int nx3, int nx2,
-                                            int nx1) noexcept;
-  PORTABLE_FUNCTION void NewPortableMDArray(T *data, int nx5, int nx4, int nx3,
-                                            int nx2, int nx1) noexcept;
-  PORTABLE_FUNCTION void NewPortableMDArray(T *data, int nx6, int nx5, int nx4,
-                                            int nx3, int nx2, int nx1) noexcept;
-
+  template <typename... NXs>
+  PORTABLE_FUNCTION void NewPortableMDArray(T *data, NXs... nxs) noexcept {
+    pdata_ = data;
+    nxs_ = nx_arr(nxs...);
+    rank_ = sizeof...(NXs);
+  }
   // public function to swap underlying data pointers of two equally-sized
   // arrays
   void SwapPortableMDArray(PortableMDArray<T> &array2);
 
   // functions to get array dimensions
-  PORTABLE_FORCEINLINE_FUNCTION int GetDim1() const { return nx1_; }
-  PORTABLE_FORCEINLINE_FUNCTION int GetDim2() const { return nx2_; }
-  PORTABLE_FORCEINLINE_FUNCTION int GetDim3() const { return nx3_; }
-  PORTABLE_FORCEINLINE_FUNCTION int GetDim4() const { return nx4_; }
-  PORTABLE_FORCEINLINE_FUNCTION int GetDim5() const { return nx5_; }
-  PORTABLE_FORCEINLINE_FUNCTION int GetDim6() const { return nx6_; }
+  template <std::size_t I>
+  PORTABLE_FUNCTION constexpr auto GetDim() const {
+    return nxs_[I];
+  }
+  // PORTABLE_FUNCTION constexpr GetSize() const { return computeSize()}
+
+  // legacy API, TODO: deprecate
+  PORTABLE_FORCEINLINE_FUNCTION int GetDim1() const { return GetDim<0>(); }
+  PORTABLE_FORCEINLINE_FUNCTION int GetDim2() const { return GetDim<1>(); }
+  PORTABLE_FORCEINLINE_FUNCTION int GetDim3() const { return GetDim<2>(); }
+  PORTABLE_FORCEINLINE_FUNCTION int GetDim4() const { return GetDim<3>(); }
+  PORTABLE_FORCEINLINE_FUNCTION int GetDim5() const { return GetDim<4>(); }
+  PORTABLE_FORCEINLINE_FUNCTION int GetDim6() const { return GetDim<5>(); }
   PORTABLE_INLINE_FUNCTION int GetDim(size_t i) const {
     // TODO: remove if performance cirtical
     assert(0 < i && i <= 6 && "PortableMDArrays are max 6D");
@@ -113,48 +162,21 @@ class PortableMDArray {
     return -1;
   }
 
-  // a function to get the total size of the array
   PORTABLE_FORCEINLINE_FUNCTION int GetSize() const {
-    return nx1_ * nx2_ * nx3_ * nx4_ * nx5_ * nx6_;
+    return std::accumulate(nxs_.cbegin(), nxs_.cend(), 1,
+                           std::multiplies<std::size_t>());
   }
   PORTABLE_FORCEINLINE_FUNCTION std::size_t GetSizeInBytes() const {
-    return nx1_ * nx2_ * nx3_ * nx4_ * nx5_ * nx6_ * sizeof(T);
-  }
-
-  PORTABLE_INLINE_FUNCTION size_t GetRank() const {
-    for (int i = 6; i >= 1; i--) {
-      if (GetDim(i) > 1) return i;
-    }
-    return 0;
+    return GetSize() * sizeof(T);
   }
 
-  PORTABLE_INLINE_FUNCTION void Reshape(int nx6, int nx5, int nx4, int nx3,
-                                        int nx2, int nx1) {
-    assert(nx6 * nx5 * nx4 * nx3 * nx2 * nx1 == GetSize());
-    nx1_ = nx1;
-    nx2_ = nx2;
-    nx3_ = nx3;
-    nx4_ = nx4;
-    nx5_ = nx5;
-    nx6_ = nx6;
+  PORTABLE_INLINE_FUNCTION size_t GetRank() const { return rank_; }
+  template <typename... NXs>
+  PORTABLE_INLINE_FUNCTION void Reshape(NXs... nxs) {
+    assert(nx_mul(nxs...) == GetSize());
+    nxs_ = nx_arr(nxs...);
+    rank_ = sizeof...(NXs);
   }
-  PORTABLE_INLINE_FUNCTION void Reshape(int nx5, int nx4, int nx3, int nx2,
-                                        int nx1) {
-    Reshape(1, nx5, nx4, nx3, nx2, nx1);
-  }
-  PORTABLE_INLINE_FUNCTION void Reshape(int nx4, int nx3, int nx2, int nx1) {
-    Reshape(1, 1, nx4, nx3, nx2, nx1);
-  }
-  PORTABLE_INLINE_FUNCTION void Reshape(int nx3, int nx2, int nx1) {
-    Reshape(1, 1, 1, nx3, nx2, nx1);
-  }
-  PORTABLE_INLINE_FUNCTION void Reshape(int nx2, int nx1) {
-    Reshape(1, 1, 1, 1, nx2, nx1);
-  }
-  PORTABLE_INLINE_FUNCTION void Reshape(int nx1) {
-    Reshape(1, 1, 1, 1, 1, nx1);
-  }
-
   PORTABLE_FORCEINLINE_FUNCTION bool IsShallowSlice() { return true; }
   PORTABLE_FORCEINLINE_FUNCTION bool IsEmpty() { return GetSize() < 1; }
   // "getter" function to access private data member
@@ -172,61 +194,16 @@ class PortableMDArray {
   // operator[]
 
   // "non-const variants" called for "PortableMDArray<T>()" provide read/write
-  // access via returning by reference, enabling assignment on returned l-value,
-  // e.g.: a(3) = 3.0;
-  PORTABLE_FORCEINLINE_FUNCTION T &operator()() { return pdata_[0]; }
-
-  PORTABLE_FORCEINLINE_FUNCTION T &operator()(const int n) { return pdata_[n]; }
-  // "const variants" called for "const PortableMDArray<T>" returns T by value,
-  // since T is typically a built-in type (versus "const T &" to avoid copying
-  // for general types)
-  PORTABLE_FORCEINLINE_FUNCTION T &operator()() const { return pdata_[0]; }
-  PORTABLE_FORCEINLINE_FUNCTION T &operator()(const int n) const {
-    return pdata_[n];
-  }
-  PORTABLE_FORCEINLINE_FUNCTION T &operator()(const int n, const int i) {
-    return pdata_[i + nx1_ * n];
-  }
-  PORTABLE_FORCEINLINE_FUNCTION T &operator()(const int n, const int i) const {
-    return pdata_[i + nx1_ * n];
-  }
-  PORTABLE_FORCEINLINE_FUNCTION T &operator()(const int n, const int j,
-                                              const int i) {
-    return pdata_[i + nx1_ * (j + nx2_ * n)];
-  }
-  PORTABLE_FORCEINLINE_FUNCTION T &operator()(const int n, const int j,
-                                              const int i) const {
-    return pdata_[i + nx1_ * (j + nx2_ * n)];
-  }
-  PORTABLE_FORCEINLINE_FUNCTION T &operator()(const int n, const int k,
-                                              const int j, const int i) {
-    return pdata_[i + nx1_ * (j + nx2_ * (k + nx3_ * n))];
+  // access via returning by reference, enabling assignment on returned
+  // l-value, e.g.: a(3) = 3.0;
+  template <typename... Is>
+  PORTABLE_FORCEINLINE_FUNCTION T &operator()(const Is... idxs) {
+    return pdata_[compute_index(nxs_, idxs...)];
   }
-  PORTABLE_FORCEINLINE_FUNCTION T &operator()(const int n, const int k,
-                                              const int j, const int i) const {
-    return pdata_[i + nx1_ * (j + nx2_ * (k + nx3_ * n))];
-  }
-  PORTABLE_FORCEINLINE_FUNCTION T &
-  operator()(const int m, const int n, const int k, const int j, const int i) {
-    return pdata_[i + nx1_ * (j + nx2_ * (k + nx3_ * (n + nx4_ * m)))];
-  }
-  PORTABLE_FORCEINLINE_FUNCTION T &operator()(const int m, const int n,
-                                              const int k, const int j,
-                                              const int i) const {
-    return pdata_[i + nx1_ * (j + nx2_ * (k + nx3_ * (n + nx4_ * m)))];
-  }
-  // int l?, int o?
-  PORTABLE_FORCEINLINE_FUNCTION T &operator()(const int p, const int m,
-                                              const int n, const int k,
-                                              const int j, const int i) {
-    return pdata_[i +
-                  nx1_ * (j + nx2_ * (k + nx3_ * (n + nx4_ * (m + nx5_ * p))))];
-  }
-  PORTABLE_FORCEINLINE_FUNCTION T &operator()(const int p, const int m,
-                                              const int n, const int k,
-                                              const int j, const int i) const {
-    return pdata_[i +
-                  nx1_ * (j + nx2_ * (k + nx3_ * (n + nx4_ * (m + nx5_ * p))))];
+
+  template <typename... Is>
+  PORTABLE_FORCEINLINE_FUNCTION T &operator()(const Is... idxs) const {
+    return pdata_[compute_index(nxs_, idxs...)];
   }
 
   PortableMDArray<T> &operator*=(T scale) {
@@ -263,19 +240,16 @@ class PortableMDArray {
 
  private:
   T *pdata_;
-  int nx1_, nx2_, nx3_, nx4_, nx5_, nx6_;
+  narr nxs_;
+  int rank_;
 };
 
 // copy constructor (does a shallow copy)
 
 template <typename T>
 PortableMDArray<T>::PortableMDArray(const PortableMDArray<T> &src) noexcept {
-  nx1_ = src.nx1_;
-  nx2_ = src.nx2_;
-  nx3_ = src.nx3_;
-  nx4_ = src.nx4_;
-  nx5_ = src.nx5_;
-  nx6_ = src.nx6_;
+  nxs_ = src.nxs_;
+  rank_ = src.rank_;
   if (src.pdata_) pdata_ = src.pdata_;
 }
 
@@ -285,12 +259,8 @@ template <typename T>
 PortableMDArray<T> &
 PortableMDArray<T>::operator=(const PortableMDArray<T> &src) noexcept {
   if (this != &src) {
-    nx1_ = src.nx1_;
-    nx2_ = src.nx2_;
-    nx3_ = src.nx3_;
-    nx4_ = src.nx4_;
-    nx5_ = src.nx5_;
-    nx6_ = src.nx6_;
+    nxs_ = src.nxs;
+    rank_ = src.rank_;
     pdata_ = src.pdata_;
   }
   return *this;
@@ -300,19 +270,17 @@ PortableMDArray<T>::operator=(const PortableMDArray<T> &src) noexcept {
 // note this POINTER equivalence, not data equivalence
 template <typename T>
 bool PortableMDArray<T>::operator==(const PortableMDArray<T> &rhs) const {
-  return (pdata_ == rhs.pdata_ && nx1_ == rhs.nx1_ && nx2_ == rhs.nx2_ &&
-          nx3_ == rhs.nx3_ && nx4_ == rhs.nx4_ && nx5_ == rhs.nx5_ &&
-          nx6_ == rhs.nx6_);
+  return (pdata_ == rhs.pdata_ && nxs_ == rhs.nxs_); // NB rank is implied
 }
 
 //----------------------------------------------------------------------------------------
 //! \fn PortableMDArray::InitWithShallowSlice()
-//  \brief shallow copy of nvar elements in dimension dim of an array, starting
-//  at index=indx. Copies pointer to data, but not data itself.
+//  \brief shallow copy of nvar elements in dimension dim of an array,
+//  starting at index=indx. Copies pointer to data, but not data itself.
 
 //  Shallow slice is only able to address the "nvar" range in "dim", and all
-//  entries of the src array for d<dim (cannot access any nx4=2, etc. entries if
-//  dim=3 for example)
+//  entries of the src array for d<dim (cannot access any nx4=2, etc. entries
+//  if dim=3 for example)
 
 template <typename T>
 PORTABLE_FUNCTION void
@@ -320,162 +288,25 @@ PortableMDArray<T>::InitWithShallowSlice(const PortableMDArray<T> &src,
                                          const int dim, const int indx,
                                          const int nvar) {
   pdata_ = src.pdata_;
-  if (dim == 6) {
-    nx6_ = nvar;
-    nx5_ = src.nx5_;
-    nx4_ = src.nx4_;
-    nx3_ = src.nx3_;
-    nx2_ = src.nx2_;
-    nx1_ = src.nx1_;
-    pdata_ += indx * (nx1_ * nx2_ * nx3_ * nx4_ * nx5_);
-  } else if (dim == 5) {
-    nx6_ = 1;
-    nx5_ = nvar;
-    nx4_ = src.nx4_;
-    nx3_ = src.nx3_;
-    nx2_ = src.nx2_;
-    nx1_ = src.nx1_;
-    pdata_ += indx * (nx1_ * nx2_ * nx3_ * nx4_);
-  } else if (dim == 4) {
-    nx6_ = 1;
-    nx5_ = 1;
-    nx4_ = nvar;
-    nx3_ = src.nx3_;
-    nx2_ = src.nx2_;
-    nx1_ = src.nx1_;
-    pdata_ += indx * (nx1_ * nx2_ * nx3_);
-  } else if (dim == 3) {
-    nx6_ = 1;
-    nx5_ = 1;
-    nx4_ = 1;
-    nx3_ = nvar;     // nx3
-    nx2_ = src.nx2_; // nx2
-    nx1_ = src.nx1_; // nx1
-    pdata_ += indx * (nx1_ * nx2_);
-  } else if (dim == 2) {
-    nx6_ = 1;
-    nx5_ = 1;
-    nx4_ = 1;
-    nx3_ = 1;
-    nx2_ = nvar;
-    nx1_ = src.nx1_;
-    pdata_ += indx * (nx1_);
-  } else if (dim == 1) {
-    nx6_ = 1;
-    nx5_ = 1;
-    nx4_ = 1;
-    nx3_ = 1;
-    nx2_ = 1;
-    nx1_ = nvar;
-    pdata_ += indx;
-  }
-  return;
-}
-
-//----------------------------------------------------------------------------------------
-//! \fn PortableMDArray::NewPortableMDArray()
-//  \brief allocate new 1D array with elements initialized to zero.
-
-template <typename T>
-PORTABLE_FUNCTION void
-PortableMDArray<T>::NewPortableMDArray(T *data, int nx1) noexcept {
-  nx1_ = nx1;
-  nx2_ = 1;
-  nx3_ = 1;
-  nx4_ = 1;
-  nx5_ = 1;
-  nx6_ = 1;
-  pdata_ = data;
-}
-
-//----------------------------------------------------------------------------------------
-//! \fn PortableMDArray::NewPortableMDArray()
-//  \brief 2d data allocation
-
-template <typename T>
-PORTABLE_FUNCTION void
-PortableMDArray<T>::NewPortableMDArray(T *data, int nx2, int nx1) noexcept {
-  nx1_ = nx1;
-  nx2_ = nx2;
-  nx3_ = 1;
-  nx4_ = 1;
-  nx5_ = 1;
-  nx6_ = 1;
-  pdata_ = data;
-}
-
-//----------------------------------------------------------------------------------------
-//! \fn PortableMDArray::NewPortableMDArray()
-//  \brief 3d data allocation
+  std::size_t offs = indx;
+  nxs_[dim - 1] = nvar;
 
-template <typename T>
-PORTABLE_FUNCTION void
-PortableMDArray<T>::NewPortableMDArray(T *data, int nx3, int nx2,
-                                       int nx1) noexcept {
-  nx1_ = nx1;
-  nx2_ = nx2;
-  nx3_ = nx3;
-  nx4_ = 1;
-  nx5_ = 1;
-  nx6_ = 1;
-  pdata_ = data;
-}
-
-//----------------------------------------------------------------------------------------
-//! \fn PortableMDArray::NewPortableMDArray()
-//  \brief 4d data allocation
-
-template <typename T>
-PORTABLE_FUNCTION void
-PortableMDArray<T>::NewPortableMDArray(T *data, int nx4, int nx3, int nx2,
-                                       int nx1) noexcept {
-  nx1_ = nx1;
-  nx2_ = nx2;
-  nx3_ = nx3;
-  nx4_ = nx4;
-  nx5_ = 1;
-  nx6_ = 1;
-  pdata_ = data;
-}
-
-//----------------------------------------------------------------------------------------
-//! \fn PortableMDArray::NewPortableMDArray()
-//  \brief 5d data allocation
-
-template <typename T>
-PORTABLE_FUNCTION void
-PortableMDArray<T>::NewPortableMDArray(T *data, int nx5, int nx4, int nx3,
-                                       int nx2, int nx1) noexcept {
-  nx1_ = nx1;
-  nx2_ = nx2;
-  nx3_ = nx3;
-  nx4_ = nx4;
-  nx5_ = nx5;
-  nx6_ = 1;
-  pdata_ = data;
-}
-
-//----------------------------------------------------------------------------------------
-//! \fn PortableMDArray::NewPortableMDArray()
-//  \brief 6d data allocation
+  for (std::size_t i = 0; i < dim - 1; ++i) {
+    nxs_[i] = src.nxs_[i];
+    offs *= nxs_[i];
+  }
+  for (std::size_t i = dim; i < MAXDIM; ++i) {
+    nxs_[i] = 1;
+  }
+  pdata_ += offs;
 
-template <typename T>
-PORTABLE_FUNCTION void
-PortableMDArray<T>::NewPortableMDArray(T *data, int nx6, int nx5, int nx4,
-                                       int nx3, int nx2, int nx1) noexcept {
-  nx1_ = nx1;
-  nx2_ = nx2;
-  nx3_ = nx3;
-  nx4_ = nx4;
-  nx5_ = nx5;
-  nx6_ = nx6;
-  pdata_ = data;
+  return;
 }
 
 //----------------------------------------------------------------------------------------
 //! \fn PortableMDArray::SwapPortableMDArray()
-//  \brief  swap pdata_ pointers of two equally sized PortableMDArrays (shallow
-//  swap)
+//  \brief  swap pdata_ pointers of two equally sized PortableMDArrays
+//  (shallow swap)
 // Does not allocate memory for either array
 
 template <typename T>

From 7a8e259eb0a27258b27e28ce516dfbab9e054c2e Mon Sep 17 00:00:00 2001
From: Christopher Mauney <mauneyc@pn2305185.lanl.gov>
Date: Mon, 2 Oct 2023 15:24:52 -0600
Subject: [PATCH 02/22] spelling + comments

---
 ports-of-call/portable_arrays.hpp | 36 ++++++++++++++++++++++++-------
 1 file changed, 28 insertions(+), 8 deletions(-)

diff --git a/ports-of-call/portable_arrays.hpp b/ports-of-call/portable_arrays.hpp
index 394e4ec6..54be08f8 100644
--- a/ports-of-call/portable_arrays.hpp
+++ b/ports-of-call/portable_arrays.hpp
@@ -35,40 +35,52 @@
 #include <type_traits>
 #include <utility> // swap()
 
+// maximum number of dimensions
 constexpr std::size_t MAXDIM = 6;
+// array type of dimensions/strides
 using narr = std::array<std::size_t, MAXDIM>;
+
 namespace detail {
 
+// set_value end-case
 template <std::size_t Ind, typename NX>
 constexpr void set_value(narr &ndat, NX value) {
   ndat[Ind] = value;
 }
 
+// set_valueS end case, dispatch to set_value
 template <std::size_t Ind, typename NX>
 constexpr void set_values(narr &ndat, NX value) {
   set_value<Ind>(ndat, value);
 }
 
+// set_values general case, set `Ind` position of array (tracking from RIGHT)
+// with head of size list `NX`
 template <std::size_t Ind, typename NX, typename... NXs>
 constexpr void set_values(narr &ndat, NX value, NXs... nxs) {
   set_value<Ind>(ndat, value);
   set_values<Ind - 1>(ndat, nxs...);
 }
 
+// compute_index base case, i.e. fastest moving index
 template <std::size_t Ind>
-size_t computeIndex(const narr &nd, const size_t index) {
+size_t compute_index(const narr &nd, const size_t index) {
   return index;
 }
+
+// compute_index general case, computing slower moving index strides
 template <std::size_t Ind, typename... Tail>
-size_t computeIndex(const narr &nd, const size_t index, const Tail... tail) {
-  return index * nd[Ind] + computeIndex<Ind + 1>(nd, tail...);
+size_t compute_index(const narr &nd, const size_t index, const Tail... tail) {
+  return index * nd[Ind] + compute_index<Ind + 1>(nd, tail...);
 }
 
+// variadic multiplication, recursive end
 template <typename NX>
 PORTABLE_INLINE_FUNCTION auto varmul(NX v) {
   return v;
 }
 
+// variadic value parameter mutliplication
 template <typename NX, typename... NXs>
 PORTABLE_INLINE_FUNCTION auto varmul(NX v, NXs... vs) {
   return v * varmul(vs...);
@@ -76,21 +88,30 @@ PORTABLE_INLINE_FUNCTION auto varmul(NX v, NXs... vs) {
 
 } // namespace detail
 
+// driver of nx array creation.
+// Note we iterate from the RIGHT, to match
+// what was done explicilt prior.
 template <typename... NXs>
 PORTABLE_INLINE_FUNCTION auto nx_arr(NXs... nxs) {
   narr r;
   constexpr std::size_t NP = sizeof...(NXs);
   detail::set_values<NP - 1>(r, nxs...);
+  // fill "empty" dims with 1
   for (auto i = NP; i < MAXDIM; ++i)
     r[i] = 1;
   return r;
 }
 
+// compute index driver.
 template <typename... Indicies>
-size_t compute_index(const narr &nd, const Indicies... idxs) {
-  return 0 + detail::computeIndex<0>(nd, idxs...);
+PORTABLE_INLINE_FUNCTION std::size_t compute_index(const narr &nd,
+                                                   const Indicies... idxs) {
+  // adding `0` if sizeof...(Indicies) == 0
+  return 0 + detail::compute_index<0>(nd, idxs...);
 }
 
+// recursive variadic multiply
+// multiply a set of variadic values
 template <typename... NXs>
 PORTABLE_INLINE_FUNCTION auto nx_mul(NXs... nxs) {
   return detail::varmul(nxs...);
@@ -99,11 +120,11 @@ PORTABLE_INLINE_FUNCTION auto nx_mul(NXs... nxs) {
 template <typename T>
 class PortableMDArray {
  public:
+  // explicit initialization of objects
   PORTABLE_FUNCTION PortableMDArray(T *data, narr nxa, std::size_t rank)
       : pdata_(data), nxs_(nxa), rank_(rank) {}
-  // rank_{std::distance(nxs_.cbegin(), std::find(nxs_.cbegin(), nxs_.cend(),
-  // 1))} {}
 
+  // variadic ctor, dispatch to explicit constructor
   template <typename... NXs>
   PORTABLE_FUNCTION PortableMDArray(T *p, NXs... nxs) noexcept
       : PortableMDArray(p, nx_arr(nxs...), sizeof...(NXs)) {}
@@ -133,7 +154,6 @@ class PortableMDArray {
   PORTABLE_FUNCTION constexpr auto GetDim() const {
     return nxs_[I];
   }
-  // PORTABLE_FUNCTION constexpr GetSize() const { return computeSize()}
 
   // legacy API, TODO: deprecate
   PORTABLE_FORCEINLINE_FUNCTION int GetDim1() const { return GetDim<0>(); }

From 4c13c2c76e51dfe4c523857127d93403d44365a8 Mon Sep 17 00:00:00 2001
From: Christopher Mauney <mauneyc@pn2305185.lanl.gov>
Date: Tue, 3 Oct 2023 10:24:50 -0600
Subject: [PATCH 03/22] added PORTABLE_ function preambles to missing ones

---
 ports-of-call/portable_arrays.hpp | 13 ++++++++-----
 1 file changed, 8 insertions(+), 5 deletions(-)

diff --git a/ports-of-call/portable_arrays.hpp b/ports-of-call/portable_arrays.hpp
index 54be08f8..940dd580 100644
--- a/ports-of-call/portable_arrays.hpp
+++ b/ports-of-call/portable_arrays.hpp
@@ -44,33 +44,36 @@ namespace detail {
 
 // set_value end-case
 template <std::size_t Ind, typename NX>
-constexpr void set_value(narr &ndat, NX value) {
+PORTABLE_INLINE_FUNCTION void set_value(narr &ndat, NX value) {
   ndat[Ind] = value;
 }
 
 // set_valueS end case, dispatch to set_value
 template <std::size_t Ind, typename NX>
-constexpr void set_values(narr &ndat, NX value) {
+PORTABLE_INLINE_FUNCTION void set_values(narr &ndat, NX value) {
   set_value<Ind>(ndat, value);
 }
 
 // set_values general case, set `Ind` position of array (tracking from RIGHT)
 // with head of size list `NX`
 template <std::size_t Ind, typename NX, typename... NXs>
-constexpr void set_values(narr &ndat, NX value, NXs... nxs) {
+PORTABLE_INLINE_FUNCTION void set_values(narr &ndat, NX value, NXs... nxs) {
   set_value<Ind>(ndat, value);
   set_values<Ind - 1>(ndat, nxs...);
 }
 
 // compute_index base case, i.e. fastest moving index
 template <std::size_t Ind>
-size_t compute_index(const narr &nd, const size_t index) {
+PORTABLE_INLINE_FUNCTION size_t compute_index(const narr &nd,
+                                              const size_t index) {
   return index;
 }
 
 // compute_index general case, computing slower moving index strides
 template <std::size_t Ind, typename... Tail>
-size_t compute_index(const narr &nd, const size_t index, const Tail... tail) {
+PORTABLE_INLINE_FUNCTION size_t compute_index(const narr &nd,
+                                              const size_t index,
+                                              const Tail... tail) {
   return index * nd[Ind] + compute_index<Ind + 1>(nd, tail...);
 }
 

From 4742b0b1e4d842533f2497285b2b8594731ab920 Mon Sep 17 00:00:00 2001
From: Christopher Michael Mauney <mauneyc@lanl.gov>
Date: Tue, 3 Oct 2023 18:38:19 -0600
Subject: [PATCH 04/22] getting build/exec on GPUs

---
 CMakeLists.txt                    |  8 +++
 ports-of-call/portable_arrays.hpp | 33 ++++--------
 test/test_portsofcall.cpp         | 83 ++++++++++++++++++++-----------
 3 files changed, 74 insertions(+), 50 deletions(-)

diff --git a/CMakeLists.txt b/CMakeLists.txt
index 35a46e2f..45fe3b46 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -64,6 +64,14 @@ target_compile_features(${POCLIB}
     cxx_std_14
 )
 
+#NOTE this option is done downstream, maybe poc should control it
+set(with_cxx "$<COMPILE_LANGUAGE:CXX>")
+set(ps_kokkos "$<BOOL:${PORTABILITY_STRATEGY_KOKKOS}>")
+
+target_compile_options(${POCLIB}
+INTERFACE
+$<${with_cxx}:$<${ps_kokkos}:--expt-relaxed-constexpr>>)
+
 # TESTING
 # ----------------------------------------
 if(PORTS_OF_CALL_BUILD_TESTING)
diff --git a/ports-of-call/portable_arrays.hpp b/ports-of-call/portable_arrays.hpp
index 940dd580..b11842e1 100644
--- a/ports-of-call/portable_arrays.hpp
+++ b/ports-of-call/portable_arrays.hpp
@@ -77,18 +77,6 @@ PORTABLE_INLINE_FUNCTION size_t compute_index(const narr &nd,
   return index * nd[Ind] + compute_index<Ind + 1>(nd, tail...);
 }
 
-// variadic multiplication, recursive end
-template <typename NX>
-PORTABLE_INLINE_FUNCTION auto varmul(NX v) {
-  return v;
-}
-
-// variadic value parameter mutliplication
-template <typename NX, typename... NXs>
-PORTABLE_INLINE_FUNCTION auto varmul(NX v, NXs... vs) {
-  return v * varmul(vs...);
-}
-
 } // namespace detail
 
 // driver of nx array creation.
@@ -113,11 +101,15 @@ PORTABLE_INLINE_FUNCTION std::size_t compute_index(const narr &nd,
   return 0 + detail::compute_index<0>(nd, idxs...);
 }
 
-// recursive variadic multiply
-// multiply a set of variadic values
-template <typename... NXs>
-PORTABLE_INLINE_FUNCTION auto nx_mul(NXs... nxs) {
-  return detail::varmul(nxs...);
+// multiply reduce array
+// NOTE: we can do product of variadic parameters `vars...` as
+// `arr_mul({vars...})`
+template <typename T, std::size_t N>
+PORTABLE_INLINE_FUNCTION auto arr_mul(const std::array<T, N> &a) {
+  auto r = T{1};
+  for (auto v : a)
+    r *= v;
+  return r;
 }
 
 template <typename T>
@@ -185,10 +177,7 @@ class PortableMDArray {
     return -1;
   }
 
-  PORTABLE_FORCEINLINE_FUNCTION int GetSize() const {
-    return std::accumulate(nxs_.cbegin(), nxs_.cend(), 1,
-                           std::multiplies<std::size_t>());
-  }
+  PORTABLE_FORCEINLINE_FUNCTION int GetSize() const { return arr_mul(nxs_); }
   PORTABLE_FORCEINLINE_FUNCTION std::size_t GetSizeInBytes() const {
     return GetSize() * sizeof(T);
   }
@@ -196,7 +185,7 @@ class PortableMDArray {
   PORTABLE_INLINE_FUNCTION size_t GetRank() const { return rank_; }
   template <typename... NXs>
   PORTABLE_INLINE_FUNCTION void Reshape(NXs... nxs) {
-    assert(nx_mul(nxs...) == GetSize());
+    assert(arr_mul({nxs...}) == GetSize());
     nxs_ = nx_arr(nxs...);
     rank_ = sizeof...(NXs);
   }
diff --git a/test/test_portsofcall.cpp b/test/test_portsofcall.cpp
index 8fb615f4..e450e749 100644
--- a/test/test_portsofcall.cpp
+++ b/test/test_portsofcall.cpp
@@ -11,6 +11,8 @@
 // copies to the public, perform publicly and display publicly, and to
 // permit others to do so.
 
+#include <bits/utility.h>
+#include <functional>
 #include <ports-of-call/portability.hpp>
 #include <ports-of-call/portable_arrays.hpp>
 #include <vector>
@@ -18,20 +20,19 @@
 #define CATCH_CONFIG_RUNNER
 #include "catch2/catch.hpp"
 
-TEST_CASE("EXECUTION_IS_HOST is set correctly", 
-    "[PortsOfCall]") {
+TEST_CASE("EXECUTION_IS_HOST is set correctly", "[PortsOfCall]") {
   // testing this is maybe nontrivial?
   auto isHost = PortsOfCall::EXECUTION_IS_HOST;
 
 #if defined(PORTABILITY_STRATEGY_KOKKOS)
-  auto checkHost = std::is_same<Kokkos::DefaultExecutionSpace, Kokkos::HostSpace::execution_space>::value;
-  REQUIRE( isHost == checkHost );
+  auto checkHost = std::is_same<Kokkos::DefaultExecutionSpace,
+                                Kokkos::HostSpace::execution_space>::value;
+  REQUIRE(isHost == checkHost);
 #elif defined(PORTABILITY_STRATEGY_CUDA)
-  REQUIRE( isHost == false );
+  REQUIRE(isHost == false);
 #else
-  REQUIRE( isHost == true );
+  REQUIRE(isHost == true);
 #endif
-
 }
 
 // this test is lifted directly from `spiner`;
@@ -71,53 +72,79 @@ TEST_CASE("PortableMDArrays can be allocated from a pointer",
     aslc2.InitWithShallowSlice(a, 1, 0, 2);
     REQUIRE(aslc1 == aslc2);
   }
+}
+
+TEST_CASE("PortableMDArrays Sizes are sane", "[PortableMDArray]") {
+  using tape_t = std::vector<int>;
+
+  constexpr std::array<std::size_t, MAXDIM> nxs{2, 5, 10, 10, 5, 2};
+
+  std::vector<int> subsz(MAXDIM);
+  std::vector<tape_t> dats;
+
+  std::partial_sum(nxs.cbegin(), nxs.cend(), subsz.begin(),
+                   std::multiplies<int>());
+  for (auto i = 0; i < MAXDIM; ++i) {
+    dats.push_back(std::vector<int>(subsz[i], nxs[i]));
+  }
+
+  SECTION("Can construct correctly") {
+    auto pmd1 = PortableMDArray<int>(dats[0].data(), nxs[0]);
 
+    REQUIRE(pmd1.GetSize() == subsz[0]);
+    REQUIRE(pmd1.GetRank() == 1);
+
+    auto pmd3 = PortableMDArray<int>(dats[3].data(), nxs[0], nxs[1], nxs[2]);
+
+    REQUIRE(pmd3.GetSize() == subsz[2]);
+    REQUIRE(pmd3.GetRank() == 3);
+    REQUIRE(pmd3.GetDim1() == nxs[2]);
+    REQUIRE(pmd3.GetDim2() == nxs[1]);
+    REQUIRE(pmd3.GetDim3() == nxs[0]);
+  }
 }
 
 PORTABLE_FORCEINLINE_FUNCTION
-Real index_func(size_t i) {
-  return i*i + 2.0*i + 3.0;
-}
+Real index_func(size_t i) { return i * i + 2.0 * i + 3.0; }
 
 TEST_CASE("portableCopy works with all portability strategies",
           "[portableCopy]") {
   // number of elements
   constexpr const size_t N = 32;
   // size in bytes
-  constexpr const size_t Nb = N*sizeof(Real);
+  constexpr const size_t Nb = N * sizeof(Real);
   // vector length N on host of Real
   std::vector<Real> b(N);
   // device pointer
-  Real* a = (Real*)PORTABLE_MALLOC(Nb);
+  Real *a = (Real *)PORTABLE_MALLOC(Nb);
 
   // set device values to 0
-  portableFor("set to 0", 0, N, PORTABLE_LAMBDA(const int& i)
-  {
-    a[i] = 0.0;
-  });
-  
+  portableFor(
+      "set to 0", 0, N, PORTABLE_LAMBDA(const int &i) { a[i] = 0.0; });
+
   // set host values to reference
-  for(size_t i = 0; i < N; ++i) {
+  for (size_t i = 0; i < N; ++i) {
     b[i] = index_func(i);
   }
 
   // copy data to device pointer
   portableCopyToDevice(a, b.data(), Nb);
-  
+
   // check if device values match reference
-  int sum {0};
-  portableReduce("check portableCopy", 0, N, 0, 0, 0, 0,
-		 PORTABLE_LAMBDA(const int& i, const int &j, const int& k, int& isum)
-  {
-    if (a[i] != index_func(i)) {
-      isum += 1;
-    }
-  }, sum);
+  int sum{0};
+  portableReduce(
+      "check portableCopy", 0, N, 0, 0, 0, 0,
+      PORTABLE_LAMBDA(const int &i, const int &j, const int &k, int &isum) {
+        if (a[i] != index_func(i)) {
+          isum += 1;
+        }
+      },
+      sum);
 
   REQUIRE(sum == 0);
 
   // set b to 0
-  for (auto& v : b) {
+  for (auto &v : b) {
     v = 0.0;
   }
 

From 4bf821eb44dc27d94804cd0908024ba92ceff262 Mon Sep 17 00:00:00 2001
From: Christopher Michael Mauney <mauneyc@lanl.gov>
Date: Tue, 3 Oct 2023 20:32:14 -0600
Subject: [PATCH 05/22] fix indexing

---
 ports-of-call/portable_arrays.hpp | 30 +++++++++++++++-------
 test/test_portsofcall.cpp         | 41 +++++++++++++++++++++++++++++++
 2 files changed, 62 insertions(+), 9 deletions(-)

diff --git a/ports-of-call/portable_arrays.hpp b/ports-of-call/portable_arrays.hpp
index b11842e1..54b34cb3 100644
--- a/ports-of-call/portable_arrays.hpp
+++ b/ports-of-call/portable_arrays.hpp
@@ -62,19 +62,29 @@ PORTABLE_INLINE_FUNCTION void set_values(narr &ndat, NX value, NXs... nxs) {
   set_values<Ind - 1>(ndat, nxs...);
 }
 
+template <std::size_t N>
+PORTABLE_INLINE_FUNCTION auto stride_arr(const narr &nd) {
+  std::array<std::size_t, N> sa;
+  sa[0] = 1;
+  for (auto i = 1; i < N; ++i) {
+    sa[i] = sa[i - 1] * nd[i - 1];
+  }
+  return sa;
+}
+
 // compute_index base case, i.e. fastest moving index
-template <std::size_t Ind>
-PORTABLE_INLINE_FUNCTION size_t compute_index(const narr &nd,
-                                              const size_t index) {
+template <std::size_t Ind, std::size_t N>
+PORTABLE_INLINE_FUNCTION size_t
+compute_index(const std::array<std::size_t, N> &sa, const size_t index) {
   return index;
 }
 
 // compute_index general case, computing slower moving index strides
-template <std::size_t Ind, typename... Tail>
-PORTABLE_INLINE_FUNCTION size_t compute_index(const narr &nd,
-                                              const size_t index,
-                                              const Tail... tail) {
-  return index * nd[Ind] + compute_index<Ind + 1>(nd, tail...);
+template <std::size_t Ind, std::size_t N, typename... Tail>
+PORTABLE_INLINE_FUNCTION size_t
+compute_index(const std::array<std::size_t, N> &sa, const size_t index,
+              const Tail... tail) {
+  return index * sa[N - Ind - 1] + compute_index<Ind + 1, N>(sa, tail...);
 }
 
 } // namespace detail
@@ -97,8 +107,10 @@ PORTABLE_INLINE_FUNCTION auto nx_arr(NXs... nxs) {
 template <typename... Indicies>
 PORTABLE_INLINE_FUNCTION std::size_t compute_index(const narr &nd,
                                                    const Indicies... idxs) {
+  constexpr auto N = sizeof...(Indicies);
+  auto strides = detail::stride_arr<N>(nd);
   // adding `0` if sizeof...(Indicies) == 0
-  return 0 + detail::compute_index<0>(nd, idxs...);
+  return 0 + detail::compute_index<0, N>(strides, idxs...);
 }
 
 // multiply reduce array
diff --git a/test/test_portsofcall.cpp b/test/test_portsofcall.cpp
index e450e749..b92db205 100644
--- a/test/test_portsofcall.cpp
+++ b/test/test_portsofcall.cpp
@@ -104,6 +104,47 @@ TEST_CASE("PortableMDArrays Sizes are sane", "[PortableMDArray]") {
   }
 }
 
+template <std::size_t NX, std::size_t NXNY, typename T>
+PORTABLE_FORCEINLINE_FUNCTION Real iflat(T k, T j, T i) {
+  return i + NX * j + NXNY * k;
+}
+
+TEST_CASE("Correct portable indexing", "[PortableMDArray]") {
+  // layout
+  constexpr std::size_t NX = 32, NY = 64, NZ = 4;
+  constexpr std::size_t NC = NX * NY * NZ;
+  constexpr Real scale = 0.1;
+  // size in bytes
+  constexpr const size_t NCb = NC * sizeof(Real);
+
+  // vector length N on host of Real
+  std::vector<Real> tape_ref(NC), tape_buf(NC);
+
+  for (auto n = 0; n < NC; ++n)
+    tape_ref[n] = scale * n;
+
+  // device pointer
+  Real *tape_d = (Real *)PORTABLE_MALLOC(NCb);
+
+  auto view_d = PortableMDArray<Real>(tape_d, NZ, NY, NX);
+
+  // set device values
+  portableFor(
+      "set unique val", 0, NZ, 0, NY, 0, NX,
+      PORTABLE_LAMBDA(const int &k, const int &j, const int &i) {
+        view_d(k, j, i) = scale * iflat<NX, NX * NY>(k, j, i);
+      });
+
+  portableCopyToHost(tape_buf.data(), tape_d, NCb);
+
+  for (auto n = 0; n < NC; ++n) {
+    INFO("REF=" << tape_ref[n] << " BUF=" << tape_buf[n]);
+    REQUIRE_THAT(tape_buf[n], Catch::Matchers::WithinRel(tape_ref[n]));
+  }
+
+  PORTABLE_FREE(tape_d);
+}
+
 PORTABLE_FORCEINLINE_FUNCTION
 Real index_func(size_t i) { return i * i + 2.0 * i + 3.0; }
 

From 4863827f9cc9b2ad95def2269c89a2c5d094fd31 Mon Sep 17 00:00:00 2001
From: Christopher Michael Mauney <mauneyc@lanl.gov>
Date: Wed, 4 Oct 2023 22:17:44 -0600
Subject: [PATCH 06/22] some design changes

---
 ports-of-call/portable_arrays.hpp | 173 +++++++++++++++---------------
 1 file changed, 86 insertions(+), 87 deletions(-)

diff --git a/ports-of-call/portable_arrays.hpp b/ports-of-call/portable_arrays.hpp
index 54b34cb3..e2adab99 100644
--- a/ports-of-call/portable_arrays.hpp
+++ b/ports-of-call/portable_arrays.hpp
@@ -28,6 +28,7 @@
 #include <algorithm>
 #include <array>
 #include <assert.h>
+#include <bits/utility.h>
 #include <cstddef> // size_t
 #include <cstring> // memset()
 #include <functional>
@@ -37,82 +38,13 @@
 
 // maximum number of dimensions
 constexpr std::size_t MAXDIM = 6;
-// array type of dimensions/strides
-using narr = std::array<std::size_t, MAXDIM>;
 
 namespace detail {
+// convert `index_sequence` values to constant
+template <std::size_t I, std::size_t V>
+constexpr std::size_t to_const = V;
 
-// set_value end-case
-template <std::size_t Ind, typename NX>
-PORTABLE_INLINE_FUNCTION void set_value(narr &ndat, NX value) {
-  ndat[Ind] = value;
-}
-
-// set_valueS end case, dispatch to set_value
-template <std::size_t Ind, typename NX>
-PORTABLE_INLINE_FUNCTION void set_values(narr &ndat, NX value) {
-  set_value<Ind>(ndat, value);
-}
-
-// set_values general case, set `Ind` position of array (tracking from RIGHT)
-// with head of size list `NX`
-template <std::size_t Ind, typename NX, typename... NXs>
-PORTABLE_INLINE_FUNCTION void set_values(narr &ndat, NX value, NXs... nxs) {
-  set_value<Ind>(ndat, value);
-  set_values<Ind - 1>(ndat, nxs...);
-}
-
-template <std::size_t N>
-PORTABLE_INLINE_FUNCTION auto stride_arr(const narr &nd) {
-  std::array<std::size_t, N> sa;
-  sa[0] = 1;
-  for (auto i = 1; i < N; ++i) {
-    sa[i] = sa[i - 1] * nd[i - 1];
-  }
-  return sa;
-}
-
-// compute_index base case, i.e. fastest moving index
-template <std::size_t Ind, std::size_t N>
-PORTABLE_INLINE_FUNCTION size_t
-compute_index(const std::array<std::size_t, N> &sa, const size_t index) {
-  return index;
-}
-
-// compute_index general case, computing slower moving index strides
-template <std::size_t Ind, std::size_t N, typename... Tail>
-PORTABLE_INLINE_FUNCTION size_t
-compute_index(const std::array<std::size_t, N> &sa, const size_t index,
-              const Tail... tail) {
-  return index * sa[N - Ind - 1] + compute_index<Ind + 1, N>(sa, tail...);
-}
-
-} // namespace detail
-
-// driver of nx array creation.
-// Note we iterate from the RIGHT, to match
-// what was done explicilt prior.
-template <typename... NXs>
-PORTABLE_INLINE_FUNCTION auto nx_arr(NXs... nxs) {
-  narr r;
-  constexpr std::size_t NP = sizeof...(NXs);
-  detail::set_values<NP - 1>(r, nxs...);
-  // fill "empty" dims with 1
-  for (auto i = NP; i < MAXDIM; ++i)
-    r[i] = 1;
-  return r;
-}
-
-// compute index driver.
-template <typename... Indicies>
-PORTABLE_INLINE_FUNCTION std::size_t compute_index(const narr &nd,
-                                                   const Indicies... idxs) {
-  constexpr auto N = sizeof...(Indicies);
-  auto strides = detail::stride_arr<N>(nd);
-  // adding `0` if sizeof...(Indicies) == 0
-  return 0 + detail::compute_index<0, N>(strides, idxs...);
-}
-
+// array type of dimensions/strides
 // multiply reduce array
 // NOTE: we can do product of variadic parameters `vars...` as
 // `arr_mul({vars...})`
@@ -123,18 +55,24 @@ PORTABLE_INLINE_FUNCTION auto arr_mul(const std::array<T, N> &a) {
     r *= v;
   return r;
 }
+} // namespace detail
 
 template <typename T>
 class PortableMDArray {
  public:
   // explicit initialization of objects
-  PORTABLE_FUNCTION PortableMDArray(T *data, narr nxa, std::size_t rank)
-      : pdata_(data), nxs_(nxa), rank_(rank) {}
+  PORTABLE_FUNCTION PortableMDArray(T *data,
+                                    std::array<std::size_t, MAXDIM> extents,
+                                    std::array<std::size_t, MAXDIM> strides,
+                                    std::size_t rank) noexcept
+      : pdata_(data), nxs_(extents), strides_(strides), rank_(rank) {}
 
   // variadic ctor, dispatch to explicit constructor
-  template <typename... NXs>
-  PORTABLE_FUNCTION PortableMDArray(T *p, NXs... nxs) noexcept
-      : PortableMDArray(p, nx_arr(nxs...), sizeof...(NXs)) {}
+  template <typename... NXs, std::size_t N = sizeof...(NXs)>
+  PORTABLE_FUNCTION PortableMDArray(T *p, NXs... nxs) noexcept {
+    NewPortableMDArray(p, nxs...);
+  } //: PortableMDArray(p, make_nxs_array(nxs...), make_strides_array<N>(), N) {
+  //}
 
   // ctors
   // default ctor: simply set null PortableMDArray
@@ -149,8 +87,7 @@ class PortableMDArray {
   template <typename... NXs>
   PORTABLE_FUNCTION void NewPortableMDArray(T *data, NXs... nxs) noexcept {
     pdata_ = data;
-    nxs_ = nx_arr(nxs...);
-    rank_ = sizeof...(NXs);
+    update_layout(nxs...);
   }
   // public function to swap underlying data pointers of two equally-sized
   // arrays
@@ -189,7 +126,9 @@ class PortableMDArray {
     return -1;
   }
 
-  PORTABLE_FORCEINLINE_FUNCTION int GetSize() const { return arr_mul(nxs_); }
+  PORTABLE_FORCEINLINE_FUNCTION int GetSize() const {
+    return detail::arr_mul(nxs_);
+  }
   PORTABLE_FORCEINLINE_FUNCTION std::size_t GetSizeInBytes() const {
     return GetSize() * sizeof(T);
   }
@@ -197,9 +136,8 @@ class PortableMDArray {
   PORTABLE_INLINE_FUNCTION size_t GetRank() const { return rank_; }
   template <typename... NXs>
   PORTABLE_INLINE_FUNCTION void Reshape(NXs... nxs) {
-    assert(arr_mul({nxs...}) == GetSize());
-    nxs_ = nx_arr(nxs...);
-    rank_ = sizeof...(NXs);
+    assert(detail::arr_mul({nxs...}) == GetSize());
+    update_layout(nxs...);
   }
   PORTABLE_FORCEINLINE_FUNCTION bool IsShallowSlice() { return true; }
   PORTABLE_FORCEINLINE_FUNCTION bool IsEmpty() { return GetSize() < 1; }
@@ -222,12 +160,12 @@ class PortableMDArray {
   // l-value, e.g.: a(3) = 3.0;
   template <typename... Is>
   PORTABLE_FORCEINLINE_FUNCTION T &operator()(const Is... idxs) {
-    return pdata_[compute_index(nxs_, idxs...)];
+    return pdata_[compute_index(idxs...)];
   }
 
   template <typename... Is>
   PORTABLE_FORCEINLINE_FUNCTION T &operator()(const Is... idxs) const {
-    return pdata_[compute_index(nxs_, idxs...)];
+    return pdata_[compute_index(idxs...)];
   }
 
   PortableMDArray<T> &operator*=(T scale) {
@@ -263,8 +201,68 @@ class PortableMDArray {
                             const int indx, const int nvar);
 
  private:
+  template <typename... NX, std::size_t N = sizeof...(NX)>
+  PORTABLE_FORCEINLINE_FUNCTION auto make_nxs_array(NX... nxs) {
+    std::array<std::size_t, MAXDIM> a;
+    std::array<std::size_t, N> t{static_cast<std::size_t>(nxs)...};
+    for (auto i = 0; i < N; ++i) {
+      a[i] = t[N - i - 1];
+    }
+    for (auto i = N; i < MAXDIM; ++i) {
+      a[i] = 1;
+    }
+    return a;
+  }
+
+  template <std::size_t N>
+  PORTABLE_INLINE_FUNCTION auto make_strides_array() {
+    std::array<std::size_t, MAXDIM> a;
+    a[0] = 1;
+    for (auto i = 1; i < N; ++i) {
+      a[i] = a[i - 1] * nxs_[i - 1];
+    }
+    for (auto i = N; i < MAXDIM; ++i)
+      a[i] = 0;
+
+    return a;
+  }
+
+  // driver of nx array creation.
+  // Note we iterate from the RIGHT, to match
+  // what was done explicilt prior.
+  template <typename... NXs, std::size_t N = sizeof...(NXs)>
+  PORTABLE_INLINE_FUNCTION auto update_layout(NXs... nxs) {
+
+    rank_ = N;
+    nxs_ = make_nxs_array(nxs...);
+    strides_ = make_strides_array<N>();
+  }
+
+  // compute_index base case, i.e. fastest moving index
+  template <std::size_t Ind, std::size_t N>
+  PORTABLE_FORCEINLINE_FUNCTION size_t
+  compute_index_impl(const size_t index) const {
+    return index;
+  }
+
+  // compute_index general case, computing slower moving index strides
+  template <std::size_t Ind, std::size_t N, typename... Tail>
+  PORTABLE_FORCEINLINE_FUNCTION size_t
+  compute_index_impl(const size_t index, const Tail... tail) const {
+    return index * strides_[N - Ind - 1] +
+           compute_index_impl<Ind + 1, N>(tail...);
+  }
+  // compute index driver.
+  template <typename... Indicies, std::size_t N = sizeof...(Indicies)>
+  PORTABLE_FORCEINLINE_FUNCTION std::size_t
+  compute_index(const Indicies... idxs) const {
+    // adding `0` if sizeof...(Indicies) == 0
+    return 0 + compute_index_impl<0, N>(idxs...);
+  }
+
   T *pdata_;
-  narr nxs_;
+  std::array<std::size_t, MAXDIM> nxs_;
+  std::array<std::size_t, MAXDIM> strides_;
   int rank_;
 };
 
@@ -274,6 +272,7 @@ template <typename T>
 PortableMDArray<T>::PortableMDArray(const PortableMDArray<T> &src) noexcept {
   nxs_ = src.nxs_;
   rank_ = src.rank_;
+  strides_ = src.strides_;
   if (src.pdata_) pdata_ = src.pdata_;
 }
 

From 1492a988a3b4706e79673d03093b78db8c32c8f7 Mon Sep 17 00:00:00 2001
From: Christopher Mauney <mauneyc@pn2305185.lanl.gov>
Date: Wed, 14 Aug 2024 11:47:45 -0600
Subject: [PATCH 07/22] tests and benchmarks

---
 ports-of-call/portable_arrays.hpp | 13 +++--
 test/CMakeLists.txt               | 35 ++++++------
 test/test_portsofcall.cpp         | 88 ++++++++++++++++++++++++-------
 test/test_utilities.hpp           | 85 +++++++++++++++++++++++++++++
 4 files changed, 179 insertions(+), 42 deletions(-)
 create mode 100644 test/test_utilities.hpp

diff --git a/ports-of-call/portable_arrays.hpp b/ports-of-call/portable_arrays.hpp
index e2adab99..3346760a 100644
--- a/ports-of-call/portable_arrays.hpp
+++ b/ports-of-call/portable_arrays.hpp
@@ -28,7 +28,6 @@
 #include <algorithm>
 #include <array>
 #include <assert.h>
-#include <bits/utility.h>
 #include <cstddef> // size_t
 #include <cstring> // memset()
 #include <functional>
@@ -55,6 +54,14 @@ PORTABLE_INLINE_FUNCTION auto arr_mul(const std::array<T, N> &a) {
     r *= v;
   return r;
 }
+PORTABLE_FORCEINLINE_FUNCTION
+decltype(auto) vp_prod() {
+  return [](auto &&v) {
+    return std::accumulate(v.begin(), v.end(), 1,
+                           std::multiplies<std::size_t>());
+  };
+}
+
 } // namespace detail
 
 template <typename T>
@@ -127,7 +134,7 @@ class PortableMDArray {
   }
 
   PORTABLE_FORCEINLINE_FUNCTION int GetSize() const {
-    return detail::arr_mul(nxs_);
+    return detail::vp_prod()(nxs_);
   }
   PORTABLE_FORCEINLINE_FUNCTION std::size_t GetSizeInBytes() const {
     return GetSize() * sizeof(T);
@@ -136,7 +143,7 @@ class PortableMDArray {
   PORTABLE_INLINE_FUNCTION size_t GetRank() const { return rank_; }
   template <typename... NXs>
   PORTABLE_INLINE_FUNCTION void Reshape(NXs... nxs) {
-    assert(detail::arr_mul({nxs...}) == GetSize());
+    assert(detail::vp_prod()(std::array{nxs...}) == GetSize());
     update_layout(nxs...);
   }
   PORTABLE_FORCEINLINE_FUNCTION bool IsShallowSlice() { return true; }
diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
index cbad4bbc..7fc8a6ba 100644
--- a/test/CMakeLists.txt
+++ b/test/CMakeLists.txt
@@ -1,40 +1,35 @@
-# © 2021. Triad National Security, LLC. All rights reserved.  This
-# program was produced under U.S. Government contract 89233218CNA000001
-# for Los Alamos National Laboratory (LANL), which is operated by Triad
-# National Security, LLC for the U.S.  Department of Energy/National
-# Nuclear Security Administration. All rights in the program are
-# reserved by Triad National Security, LLC, and the U.S. Department of
-# Energy/National Nuclear Security Administration. The Government is
-# granted for itself and others acting on its behalf a nonexclusive,
-# paid-up, irrevocable worldwide license in this material to reproduce,
-# prepare derivative works, distribute copies to the public, perform
+# © 2021. Triad National Security, LLC. All rights reserved.  This program was
+# produced under U.S. Government contract 89233218CNA000001 for Los Alamos
+# National Laboratory (LANL), which is operated by Triad National Security, LLC
+# for the U.S.  Department of Energy/National Nuclear Security Administration.
+# All rights in the program are reserved by Triad National Security, LLC, and
+# the U.S. Department of Energy/National Nuclear Security Administration. The
+# Government is granted for itself and others acting on its behalf a
+# nonexclusive, paid-up, irrevocable worldwide license in this material to
+# reproduce, prepare derivative works, distribute copies to the public, perform
 # publicly and display publicly, and to permit others to do so.
 
 find_package(Catch2 REQUIRED)
 
-# this interface target is to collect
-# compile/link options for the test
+# this interface target is to collect compile/link options for the test
 add_library(portsofcall_iface INTERFACE)
-if (PORTABILITY_STRATEGY_KOKKOS)
+if(PORTABILITY_STRATEGY_KOKKOS)
   if(NOT TARGET Kokkos::kokkos)
     find_package(Kokkos REQUIRED)
   endif()
 
   target_link_libraries(portsofcall_iface INTERFACE Kokkos::kokkos)
   # this comes with ports-of-call target
-  target_compile_definitions(portsofcall_iface INTERFACE PORTABILITY_STRATEGY_KOKKOS)
+  target_compile_definitions(portsofcall_iface
+                             INTERFACE PORTABILITY_STRATEGY_KOKKOS)
 endif()
 
 target_link_libraries(portsofcall_iface INTERFACE Catch2::Catch2)
 # add unit tests
 add_executable(test_portsofcall test_portsofcall.cpp)
 
-target_link_libraries(test_portsofcall
-  PRIVATE
-    ports-of-call::ports-of-call
-    portsofcall_iface
-)
+target_link_libraries(test_portsofcall PRIVATE ports-of-call::ports-of-call
+                                               portsofcall_iface)
 
 include(Catch)
 catch_discover_tests(test_portsofcall)
-
diff --git a/test/test_portsofcall.cpp b/test/test_portsofcall.cpp
index b92db205..feb28ab7 100644
--- a/test/test_portsofcall.cpp
+++ b/test/test_portsofcall.cpp
@@ -11,14 +11,18 @@
 // copies to the public, perform publicly and display publicly, and to
 // permit others to do so.
 
-#include <bits/utility.h>
 #include <functional>
+#include <iostream>
+#include <numeric>
 #include <ports-of-call/portability.hpp>
 #include <ports-of-call/portable_arrays.hpp>
+// #include <sys/_types/_ucontext.h>
+#include <utility>
 #include <vector>
 
 #define CATCH_CONFIG_RUNNER
-#include "catch2/catch.hpp"
+#include "catch2/catch_all.hpp"
+#include "test_utilities.hpp"
 
 TEST_CASE("EXECUTION_IS_HOST is set correctly", "[PortsOfCall]") {
   // testing this is maybe nontrivial?
@@ -104,24 +108,24 @@ TEST_CASE("PortableMDArrays Sizes are sane", "[PortableMDArray]") {
   }
 }
 
-template <std::size_t NX, std::size_t NXNY, typename T>
-PORTABLE_FORCEINLINE_FUNCTION Real iflat(T k, T j, T i) {
-  return i + NX * j + NXNY * k;
-}
-
 TEST_CASE("Correct portable indexing", "[PortableMDArray]") {
   // layout
+  auto iflat = [](auto nx, auto nxny) {
+    return [=](auto k, auto j, auto i) { return i + nx * j + nxny * k; };
+  };
+
   constexpr std::size_t NX = 32, NY = 64, NZ = 4;
   constexpr std::size_t NC = NX * NY * NZ;
   constexpr Real scale = 0.1;
   // size in bytes
-  constexpr const size_t NCb = NC * sizeof(Real);
+  constexpr const std::size_t NCb = NC * sizeof(Real);
 
   // vector length N on host of Real
   std::vector<Real> tape_ref(NC), tape_buf(NC);
 
-  for (auto n = 0; n < NC; ++n)
-    tape_ref[n] = scale * n;
+  for (auto n = 0; n < NC; ++n) {
+    tape_ref[n] = scale * static_cast<Real>(n);
+  }
 
   // device pointer
   Real *tape_d = (Real *)PORTABLE_MALLOC(NCb);
@@ -132,24 +136,22 @@ TEST_CASE("Correct portable indexing", "[PortableMDArray]") {
   portableFor(
       "set unique val", 0, NZ, 0, NY, 0, NX,
       PORTABLE_LAMBDA(const int &k, const int &j, const int &i) {
-        view_d(k, j, i) = scale * iflat<NX, NX * NY>(k, j, i);
+        view_d(k, j, i) = scale * iflat(NX, NX * NY)(k, j, i);
       });
 
   portableCopyToHost(tape_buf.data(), tape_d, NCb);
 
   for (auto n = 0; n < NC; ++n) {
-    INFO("REF=" << tape_ref[n] << " BUF=" << tape_buf[n]);
+    INFO("REF=" << tape_ref[n] << " BUF=" << tape_buf[n] << " n=" << n);
     REQUIRE_THAT(tape_buf[n], Catch::Matchers::WithinRel(tape_ref[n]));
   }
 
   PORTABLE_FREE(tape_d);
 }
 
-PORTABLE_FORCEINLINE_FUNCTION
-Real index_func(size_t i) { return i * i + 2.0 * i + 3.0; }
-
 TEST_CASE("portableCopy works with all portability strategies",
-          "[portableCopy]") {
+          "[PortableMDArray]") {
+  auto index_func = [](auto i) { return i * i + 2.0 * i + 3.0; };
   // number of elements
   constexpr const size_t N = 32;
   // size in bytes
@@ -193,19 +195,67 @@ TEST_CASE("portableCopy works with all portability strategies",
   portableCopyToHost(b.data(), a, Nb);
 
   // count elements that don't match reference
-  sum = 0;
+  int nbad{0};
   for (int i = 0; i < N; ++i) {
     if (b[i] != index_func(i)) {
-      sum += 1;
+      nbad += 1;
     }
   }
   // make sure all elements match
-  REQUIRE(sum == 0);
+  INFO("nbad=" << nbad);
+  REQUIRE(nbad == 0);
 
   // free device memory
   PORTABLE_FREE(a);
 }
 
+// runs benchmarks for indexing.
+// note: to run, execute
+//    ./tests/test_portsofcall "[!benchmark]"
+// (you may need to escape the `!` char depending
+// on your shell)
+TEST_CASE("Benchmark Indexing", "[!benchmark]") {
+  constexpr std::size_t NX = 32, NY = 64, NZ = 16;
+  constexpr std::size_t NU = 256, NV = 256, NW = 128;
+
+  // construct a series of PortableMDArrays
+  // of different rank and sizes, and
+  // iterates through all contiguous indexing.
+
+  // clang-format off
+  SECTION("small") {
+    BENCHMARK("index calc 1D") 
+    { 
+      return testing::idx_contiguous_bm(std::array{NX}); 
+    };
+    BENCHMARK("index calc 2D") 
+    { 
+      return testing::idx_contiguous_bm(std::array{NX, NY}); 
+    };
+    BENCHMARK("index calc 3D") 
+    { 
+      return testing::idx_contiguous_bm(std::array{NX, NY, NZ}); 
+    };
+  }
+
+  SECTION("big") {
+    BENCHMARK("index calc 1D") 
+    { 
+      return testing::idx_contiguous_bm(std::array{NU}); 
+    };
+    BENCHMARK("index calc 2D") 
+    { 
+      return testing::idx_contiguous_bm(std::array{NU, NV}); 
+    };
+    BENCHMARK("index calc 3D") 
+    { 
+      return testing::idx_contiguous_bm(std::array{NU, NV, NW}); 
+    };
+
+  }
+  // clang-format on
+}
+
 int main(int argc, char *argv[]) {
 
 #ifdef PORTABILITY_STRATEGY_KOKKOS
diff --git a/test/test_utilities.hpp b/test/test_utilities.hpp
new file mode 100644
index 00000000..2d5ad3fb
--- /dev/null
+++ b/test/test_utilities.hpp
@@ -0,0 +1,85 @@
+#ifndef _TEST_PORTOFCALL_HPP_
+#define _TEST_PORTOFCALL_HPP_
+
+#include <__utility/integer_sequence.h>
+#include <array>
+#include <functional>
+#include <ports-of-call/portability.hpp>
+#include <ports-of-call/portable_arrays.hpp>
+#include <utility>
+
+namespace testing {
+
+template <typename Array, std::size_t... I>
+constexpr auto array_fix_impl(const Array &a, std::index_sequence<I...> is) {
+  auto mm = std::array<int, is.size() * 2>();
+  for (int i = 0; i < is.size(); ++i) {
+    mm[(i * 2)] = 0;
+    mm[(i * 2) + 1] = a[i];
+  }
+  return mm;
+}
+
+// we want to call
+//    portableFor(..., <0, NX>, <0, NY>, ...)
+// for an arbitrary set of ranges
+// this function takes an array of N sizes, and converts
+// it into an array 2N {0, NX} contiguous pairs.
+// ex:
+//    array{10,12,13} -> array{0,10,0,12,0,13}
+template <typename T, std::size_t N,
+          typename Indx = std::make_index_sequence<N>>
+PORTABLE_FORCEINLINE_FUNCTION constexpr auto
+array_fix(const std::array<T, N> &a) {
+  return array_fix_impl(a, Indx{});
+}
+
+// constructs the portableFor call
+// tparams:
+//    - View : class of data view
+//    - Array: class of container holding sizes
+//    - J... : index sequence
+// params:
+//    - View v: view of data
+//    - array ext_par: array of sizes
+// this was written ad-hoc, and could be more general
+// (e.g. pass label, functor in lambda)
+template <class View, class Array, std::size_t... J>
+PORTABLE_INLINE_FUNCTION constexpr auto pf_invoke(View v, const Array &ext_par,
+                                                  std::index_sequence<J...>) {
+
+  auto stud_arr = array_fix(ext_par);
+  portableFor(
+      "set unique val", stud_arr[J]...,
+      PORTABLE_LAMBDA(auto... is) { v(is...) = 1.0; });
+}
+
+template <class Ptr, class Array, std::size_t... I>
+PORTABLE_FORCEINLINE_FUNCTION constexpr decltype(auto)
+mdview(Ptr *d, const Array &arr, std::index_sequence<I...>) {
+  return PortableMDArray<Real>(d, arr[I]...);
+}
+
+// stand up and strip down a benchmark
+// input is a set of sizes
+template <class T, std::size_t N, typename I1 = std::make_index_sequence<N>,
+          typename I2 = std::make_index_sequence<2 * N>>
+PORTABLE_FUNCTION auto idx_contiguous_bm(const std::array<T, N> &nxa) {
+  auto nc =
+      std::accumulate(std::begin(nxa), std::end(nxa), 1, std::multiplies<T>{});
+
+  Real *tape_d = (Real *)PORTABLE_MALLOC(nc * sizeof(Real));
+
+  auto view_d = mdview(tape_d, nxa, I1{});
+
+  pf_invoke(view_d, nxa, I2{});
+
+  PORTABLE_FREE(tape_d);
+
+  // return here is arbitrary, Catch2 recommends this to avoid
+  // optimizing out a call.
+  return 1;
+}
+
+} // namespace testing
+#endif

From 551910bc1b6cab942d5905cfba0e9692d5867b5b Mon Sep 17 00:00:00 2001
From: Christopher Mauney <mauneyc@pn2305185.lanl.gov>
Date: Thu, 29 Aug 2024 10:09:39 -0600
Subject: [PATCH 08/22] store

---
 ports-of-call/portable_arrays.hpp    |  26 ++-----
 ports-of-call/utility/array.hpp      |  74 ++++++++++++++++++
 ports-of-call/utility/index_algo.hpp | 108 +++++++++++++++++++++++++++
 3 files changed, 187 insertions(+), 21 deletions(-)
 create mode 100644 ports-of-call/utility/array.hpp
 create mode 100644 ports-of-call/utility/index_algo.hpp

diff --git a/ports-of-call/portable_arrays.hpp b/ports-of-call/portable_arrays.hpp
index 3346760a..dc4a7fcf 100644
--- a/ports-of-call/portable_arrays.hpp
+++ b/ports-of-call/portable_arrays.hpp
@@ -25,12 +25,14 @@
 //  NOTE THE TRAILING INDEX INSIDE THE PARENTHESES IS INDEXED FASTEST
 
 #include "portability.hpp"
+#include "utility/array.hpp"
 #include <algorithm>
 #include <array>
 #include <assert.h>
 #include <cstddef> // size_t
 #include <cstring> // memset()
 #include <functional>
+#include <integer_sequence>
 #include <numeric>
 #include <type_traits>
 #include <utility> // swap()
@@ -43,25 +45,6 @@ namespace detail {
 template <std::size_t I, std::size_t V>
 constexpr std::size_t to_const = V;
 
-// array type of dimensions/strides
-// multiply reduce array
-// NOTE: we can do product of variadic parameters `vars...` as
-// `arr_mul({vars...})`
-template <typename T, std::size_t N>
-PORTABLE_INLINE_FUNCTION auto arr_mul(const std::array<T, N> &a) {
-  auto r = T{1};
-  for (auto v : a)
-    r *= v;
-  return r;
-}
-PORTABLE_FORCEINLINE_FUNCTION
-decltype(auto) vp_prod() {
-  return [](auto &&v) {
-    return std::accumulate(v.begin(), v.end(), 1,
-                           std::multiplies<std::size_t>());
-  };
-}
-
 } // namespace detail
 
 template <typename T>
@@ -134,7 +117,7 @@ class PortableMDArray {
   }
 
   PORTABLE_FORCEINLINE_FUNCTION int GetSize() const {
-    return detail::vp_prod()(nxs_);
+    return util::array_reduce(nxs_, std::multiplies<std::size_t>{});
   }
   PORTABLE_FORCEINLINE_FUNCTION std::size_t GetSizeInBytes() const {
     return GetSize() * sizeof(T);
@@ -143,7 +126,8 @@ class PortableMDArray {
   PORTABLE_INLINE_FUNCTION size_t GetRank() const { return rank_; }
   template <typename... NXs>
   PORTABLE_INLINE_FUNCTION void Reshape(NXs... nxs) {
-    assert(detail::vp_prod()(std::array{nxs...}) == GetSize());
+    assert(util::array_reduce(std::array{nxs...},
+                              std::multiplies<std::size_t>{}) == GetSize());
     update_layout(nxs...);
   }
   PORTABLE_FORCEINLINE_FUNCTION bool IsShallowSlice() { return true; }
diff --git a/ports-of-call/utility/array.hpp b/ports-of-call/utility/array.hpp
new file mode 100644
index 00000000..56e2508f
--- /dev/null
+++ b/ports-of-call/utility/array.hpp
@@ -0,0 +1,74 @@
+#ifndef _PORTSOFCALL_UTILITY_ARRAY_HPP_
+#define _PORTSOFCALL_UTILITY_ARRAY_HPP_
+
+#include "../portability.hpp"
+#include <array>
+#include <type_traits>
+
+namespace util {
+
+namespace detail {
+template <class T, std::size_t N, class F, std::size_t... Is>
+PORTABLE_FORCEINLINE_FUNCTION constexpr auto
+array_map_impl(std::array<T, N> const &x, F f, std::index_sequence<Is...>) {
+  return std::array{f(x[Is])...};
+}
+
+template <class T, class U, std::size_t N, class F, std::size_t... Is>
+PORTABLE_FORCEINLINE_FUNCTION constexpr auto
+array_map_impl(std::array<T, N> const &x, std::array<U, N> const &y, F f,
+               std::index_sequence<Is...>) {
+  return std::array{f(x[Is], y[Is])...};
+}
+
+template <std::size_t f, std::size_t l, class T, std::size_t N, class Op>
+PORTABLE_INLINE_FUNCTION constexpr T
+array_reduce_impl(std::array<T, N> const &x, Op op) {
+  if constexpr ((l - f) == 1)
+    return x[f];
+  else {
+    constexpr std::size_t n = l - f;
+    T left_sum = array_reduce_impl<f, f + n / 2>(x, op);
+    T right_sum = array_reduce_impl<f + n / 2, l>(x, op);
+    return op(left_sum, right_sum);
+  }
+}
+
+} // namespace detail
+
+template <class T, std::size_t N, class F>
+PORTABLE_FORCEINLINE_FUNCTION constexpr auto
+array_map(std::array<T, N> const &x, F f) {
+  return detail::array_map_impl(x, f, std::make_index_sequence<N>{});
+}
+
+template <class T, class U, std::size_t N, class F>
+PORTABLE_FORCEINLINE_FUNCTION constexpr auto
+array_map(std::array<T, N> const &x, std::array<U, N> const &y, F f) {
+  return detail::array_map_impl(x, y, f, std::make_index_sequence<N>{});
+}
+
+template <std::size_t I, class T, std::size_t N, class Op>
+PORTABLE_FORCEINLINE_FUNCTION constexpr T
+array_partial_reduce(std::array<T, N> x, T initial_value, Op op) {
+  static_assert(I <= N);
+  if constexpr (I == 0)
+    return initial_value;
+  else
+    return detail::array_reduce_impl<0, I>(x, op);
+}
+
+template <class T, std::size_t N, class Op>
+PORTABLE_FORCEINLINE_FUNCTION constexpr T array_reduce(std::array<T, N> x,
+                                                       T initial_value, Op op) {
+  return array_partial_reduce<N>(x, initial_value, op);
+}
+template <std::size_t... Is>
+PORTABLE_FORCEINLINE_FUNCTION constexpr auto
+as_array(std::index_sequence<sizeof...(Is)>) {
+  return std::array{Is...};
+}
+
+} // namespace util
+
+#endif // _PORTSOFCALL_UTILITY_ARRAY_HPP_
diff --git a/ports-of-call/utility/index_algo.hpp b/ports-of-call/utility/index_algo.hpp
new file mode 100644
index 00000000..31a5d8e8
--- /dev/null
+++ b/ports-of-call/utility/index_algo.hpp
@@ -0,0 +1,108 @@
+#ifndef _PORTSOFCALL_UTILITY_INDEX_ALGO_HPP_
+#define _PORTSOFCALL_UTILITY_INDEX_ALGO_HPP_
+
+#include "../portability.hpp"
+#include "array.hpp"
+#include <array>
+#include <numeric>
+#include <type_traits>
+
+namespace util {
+
+template <auto I, class T, std::size_t N>
+PORTABLE_FORCEINLINE_FUNCTION static constexpr auto
+get_stride(std::array<T, N> const &dim) {
+  static_assert(I < dim.size(), "Dim index is out of bounds");
+
+  // column major
+  return std::accumulate(std::begin(dim), std::begin(dim) + I, std::size_t{1},
+                         std::multiplies<std::size_t>{});
+}
+
+namespace detail {
+template <class T, std::size_t N, std::size_t... Is>
+PORTABLE_FORCEINLINE_FUNCTION static constexpr auto
+get_strides_impl(std::array<T, N> const &dim, std::index_sequence<Is...>) {
+  return std::array{get_stride<Is>(dim)...};
+}
+
+} // namespace detail
+
+template <class T, std::size_t N>
+PORTABLE_FORCEINLINE_FUNCTION static constexpr auto
+get_strides(std::array<T, N> const &dim) {
+  return detail::get_strides_impl(dim, std::make_index_sequence<dim.size()>{});
+}
+
+template <class T, std::size_t N>
+PORTABLE_FORCEINLINE_FUNCTION static constexpr auto
+fast_findex(std::array<T, N> const &ijk, std::array<T, N> const &dim,
+            std::array<T, N> const &stride) {
+  // TODO: assert ijk in bounds
+  return std::inner_product(std::begin(ijk), std::end(ijk), std::begin(stride),
+                            std::size_t{0});
+}
+
+template <class T, std::size_t N>
+PORTABLE_FORCEINLINE_FUNCTION static constexpr auto
+findex(std::array<T, N> const &ijk, std::array<T, N> const &dim) {
+  return fast_findex(ijk, dim, get_strides(dim));
+}
+
+namespace handroll {
+template <auto I, class T, std::size_t N>
+PORTABLE_FORCEINLINE_FUNCTION static constexpr auto
+get_stride(std::array<T, N> const &dim) {
+
+  // column major
+  return array_partial_reduce<I>(dim, T{1}, std::multiplies<std::size_t>{});
+}
+
+namespace detail {
+template <class T, std::size_t N, std::size_t... Is>
+PORTABLE_FORCEINLINE_FUNCTION static constexpr auto
+get_strides_impl(std::array<T, N> const &dim, std::index_sequence<Is...>) {
+  return std::array{get_stride<Is>(dim)...};
+}
+template <class T, std::size_t N>
+PORTABLE_FORCEINLINE_FUNCTION static constexpr auto
+get_strides(std::array<T, N> const &dim) {
+  return detail::get_strides_impl(dim, std::make_index_sequence<dim.size()>{});
+} // namespace handroll::detail
+template <class T, std::size_t N>
+PORTABLE_FORCEINLINE_FUNCTION static constexpr auto
+fast_findex(std::array<T, N> const &ijk, std::array<T, N> const &dim,
+            std::array<T, N> const &stride) {
+  // TODO: assert ijk in bounds
+  return array_reduce(
+      array_map(ijk, stride, [](auto a, auto b) { return a * b; }), T{1},
+      std::plus<std::size_t>{});
+}
+template <class T, std::size_t N>
+PORTABLE_FORCEINLINE_FUNCTION static constexpr auto
+findex(std::array<T, N> const &ijk, std::array<T, N> const &dim) {
+  return fast_findex(ijk, dim, get_strides(dim));
+}
+} // namespace handroll
+
+template <class T, std::size_t N>
+PORTABLE_FORCEINLINE_FUNCTION static constexpr std::array<T, N>
+fast_mindices(std::size_t idx, std::array<T, N> const &dim,
+              std::array<T, N> const &stride) {
+  std::array<T, N> mdidx;
+  for (std::int64_t i = dim.size() - 1; i >= 0; --i) {
+    mdidx[i] = idx / std::size_t(stride[i]);
+    idx -= mdidx[i] * std::size_t(stride[i]);
+  }
+  return mdidx;
+}
+
+template <class Array>
+PORTABLE_FORCEINLINE_FUNCTION static constexpr auto mindices(std::size_t idx,
+                                                             Array dim) {
+  return fast_mindices(idx, dim, get_strides(dim));
+}
+
+} // namespace util
+
+#endif // _PORTSOFCALL_UTILITY_INDEX_ALGO_HPP_

From 0f66ba178bcc89634d8b0bdb537e74adcbb5e170 Mon Sep 17 00:00:00 2001
From: Christopher Mauney <mauneyc@pn2305185.lanl.gov>
Date: Fri, 30 Aug 2024 13:45:17 -0600
Subject: [PATCH 09/22] array and index algs

---
 ports-of-call/array.hpp              |   7 ++
 ports-of-call/portable_arrays.hpp    | 155 ++++++++++++---------------
 ports-of-call/utility/array_algo.hpp |  81 +++++++++-----
 ports-of-call/utility/index_algo.hpp |  61 +++++------
 test/CMakeLists.txt                  |   3 +-
 test/test_mdidx.cpp                  | 144 +++++++++++++++++++++++++
 test/test_utilities.hpp              |  19 ++--
 7 files changed, 310 insertions(+), 160 deletions(-)
 create mode 100644 test/test_mdidx.cpp

diff --git a/ports-of-call/array.hpp b/ports-of-call/array.hpp
index 8eb90cc0..64569f01 100644
--- a/ports-of-call/array.hpp
+++ b/ports-of-call/array.hpp
@@ -3,6 +3,7 @@
 
 #include "portability.hpp"
 #include "portable_errors.hpp"
+#include "ports-of-call/utility/array_algo.hpp"
 
 #include <cstddef>
 #include <iterator>
@@ -192,6 +193,12 @@ struct array {
     }
   }
 
+  PORTABLE_FUNCTION constexpr bool operator==(const array &rhs) const {
+    for (size_type i = 0; i < N; ++i)
+      if (arr[i] != rhs.arr[i]) return false;
+    return true;
+  }
+
   //! swap the array contents with another
   PORTABLE_FUNCTION constexpr void swap(array &other) {
     using std::swap;
diff --git a/ports-of-call/portable_arrays.hpp b/ports-of-call/portable_arrays.hpp
index d96393ac..72ef43f6 100644
--- a/ports-of-call/portable_arrays.hpp
+++ b/ports-of-call/portable_arrays.hpp
@@ -24,15 +24,15 @@
 //  [N4xN3xN2xN1] are accessed as:  A(n,k,j,i) = A[i + N1*(j + N2*(k + N3*n))]
 //  NOTE THE TRAILING INDEX INSIDE THE PARENTHESES IS INDEXED FASTEST
 
+#include "array.hpp"
 #include "portability.hpp"
-#include "utility/array.hpp"
+#include "utility/array_algo.hpp"
 #include <algorithm>
 #include <array>
 #include <assert.h>
 #include <cstddef> // size_t
 #include <cstring> // memset()
 #include <functional>
-#include <integer_sequence>
 #include <numeric>
 #include <type_traits>
 #include <utility> // swap()
@@ -47,12 +47,17 @@ constexpr std::size_t to_const = V;
 
 } // namespace detail
 
-template <typename T>
+template <class T, auto N>
+using Array = PortsOfCall::array<T, N>;
+
+template <auto N>
+using IArray = Array<std::size_t, N>;
+
+template <typename T, std::size_t D = MAXDIM>
 class PortableMDArray {
  public:
   // explicit initialization of objects
-  PORTABLE_FUNCTION PortableMDArray(T *data, std::array<std::size_t, MAXDIM> extents,
-                                    std::array<std::size_t, MAXDIM> strides,
+  PORTABLE_FUNCTION PortableMDArray(T *data, IArray<D> extents, IArray<D> strides,
                                     std::size_t rank) noexcept
       : pdata_(data), nxs_(extents), strides_(strides), rank_(rank) {}
 
@@ -69,8 +74,23 @@ class PortableMDArray {
   PortableMDArray() noexcept : pdata_(nullptr), nxs_{{0}}, rank_{0} {}
   // define copy constructor and overload assignment operator so both do deep
   // copies.
-  PortableMDArray(const PortableMDArray<T> &t) noexcept;
-  PortableMDArray<T> &operator=(const PortableMDArray<T> &t) noexcept;
+  //  PortableMDArray(const PortableMDArray<T> &t) noexcept;
+  //  PortableMDArray<T> &operator=(const PortableMDArray<T> &t) noexcept;
+  PortableMDArray(const PortableMDArray<T, D> &src) noexcept {
+    nxs_ = src.nxs_;
+    rank_ = src.rank_;
+    strides_ = src.strides_;
+    if (src.pdata_) pdata_ = src.pdata_;
+  }
+
+  PortableMDArray<T, D> &operator=(const PortableMDArray<T, D> &src) noexcept {
+    if (this != &src) {
+      nxs_ = src.nxs;
+      rank_ = src.rank_;
+      pdata_ = src.pdata_;
+    }
+    return *this;
+  }
 
   template <typename... NXs>
   PORTABLE_FUNCTION void NewPortableMDArray(T *data, NXs... nxs) noexcept {
@@ -79,7 +99,12 @@ class PortableMDArray {
   }
   // public function to swap underlying data pointers of two equally-sized
   // arrays
-  void SwapPortableMDArray(PortableMDArray<T> &array2);
+  // void SwapPortableMDArray(PortableMDArray<T> &array2);
+
+  void SwapPortableMDArray(PortableMDArray<T, D> &array2) {
+    std::swap(pdata_, array2.pdata_);
+    return;
+  }
 
   // functions to get array dimensions
   template <std::size_t I>
@@ -115,16 +140,16 @@ class PortableMDArray {
   }
 
   PORTABLE_FORCEINLINE_FUNCTION int GetSize() const {
-    return util::array_reduce(nxs_, std::multiplies<std::size_t>{});
+    return util::array_reduce(nxs_, 1, std::multiplies<std::size_t>{});
   }
   PORTABLE_FORCEINLINE_FUNCTION std::size_t GetSizeInBytes() const {
     return GetSize() * sizeof(T);
   }
 
-  PORTABLE_INLINE_FUNCTION size_t GetRank() const { return rank_; }
+  PORTABLE_INLINE_FUNCTION std::size_t GetRank() const { return rank_; }
   template <typename... NXs>
   PORTABLE_INLINE_FUNCTION void Reshape(NXs... nxs) {
-    assert(util::array_reduce(std::array{nxs...}, std::multiplies<std::size_t>{}) ==
+    assert(util::array_reduce(IArray<D>{nxs...}, 1, std::multiplies<std::size_t>{}) ==
            GetSize());
     update_layout(nxs...);
   }
@@ -158,19 +183,19 @@ class PortableMDArray {
     return pdata_[compute_index(idxs...)];
   }
 
-  PortableMDArray<T> &operator*=(T scale) {
+  PortableMDArray<T, D> &operator*=(T scale) {
     std::transform(pdata_, pdata_ + GetSize(), pdata_,
                    [scale](T val) { return scale * val; });
     return *this;
   }
 
-  PortableMDArray<T> &operator+=(const PortableMDArray<T> &other) {
+  PortableMDArray<T, D> &operator+=(const PortableMDArray<T, D> &other) {
     assert(GetSize() == other.GetSize());
     std::transform(pdata_, pdata_ + GetSize(), other.pdata_, pdata_, std::plus<T>());
     return *this;
   }
 
-  PortableMDArray<T> &operator-=(const PortableMDArray<T> &other) {
+  PortableMDArray<T, D> &operator-=(const PortableMDArray<T, D> &other) {
     assert(GetSize() == other.GetSize());
     std::transform(pdata_, pdata_ + GetSize(), other.pdata_, pdata_, std::minus<T>());
     return *this;
@@ -178,23 +203,40 @@ class PortableMDArray {
 
   // Checks that arrays point to same data with same shape
   // note this POINTER equivalence, not data equivalence
-  bool operator==(const PortableMDArray<T> &other) const;
-  bool operator!=(const PortableMDArray<T> &other) const { return !(*this == other); }
+  bool operator==(const PortableMDArray<T, D> &rhs) const {
+    return (pdata_ == rhs.pdata_ && nxs_ == rhs.nxs_); // NB rank is implied
+  }
 
-  // (deferred) initialize an array with slice from another array
-  PORTABLE_FUNCTION
-  void InitWithShallowSlice(const PortableMDArray<T> &src, const int dim, const int indx,
-                            const int nvar);
+  bool operator!=(const PortableMDArray<T, D> &other) const { return !(*this == other); }
+
+  PORTABLE_FUNCTION void InitWithShallowSlice(const PortableMDArray<T, D> &src,
+                                              const int dim, const int indx,
+                                              const int nvar) {
+    pdata_ = src.pdata_;
+    std::size_t offs = indx;
+    nxs_[dim - 1] = nvar;
+
+    for (std::size_t i = 0; i < dim - 1; ++i) {
+      nxs_[i] = src.nxs_[i];
+      offs *= nxs_[i];
+    }
+    for (std::size_t i = dim; i < MAXDIM; ++i) {
+      nxs_[i] = 1;
+    }
+    pdata_ += offs;
+
+    return;
+  }
 
  private:
   template <typename... NX, std::size_t N = sizeof...(NX)>
   PORTABLE_FORCEINLINE_FUNCTION auto make_nxs_array(NX... nxs) {
-    std::array<std::size_t, MAXDIM> a;
-    std::array<std::size_t, N> t{static_cast<std::size_t>(nxs)...};
+    IArray<D> a;
+    IArray<N> t{static_cast<std::size_t>(nxs)...};
     for (auto i = 0; i < N; ++i) {
       a[i] = t[N - i - 1];
     }
-    for (auto i = N; i < MAXDIM; ++i) {
+    for (auto i = N; i < D; ++i) {
       a[i] = 1;
     }
     return a;
@@ -202,12 +244,12 @@ class PortableMDArray {
 
   template <std::size_t N>
   PORTABLE_INLINE_FUNCTION auto make_strides_array() {
-    std::array<std::size_t, MAXDIM> a;
+    IArray<D> a;
     a[0] = 1;
     for (auto i = 1; i < N; ++i) {
       a[i] = a[i - 1] * nxs_[i - 1];
     }
-    for (auto i = N; i < MAXDIM; ++i)
+    for (auto i = N; i < D; ++i)
       a[i] = 0;
 
     return a;
@@ -244,41 +286,12 @@ class PortableMDArray {
   }
 
   T *pdata_;
-  std::array<std::size_t, MAXDIM> nxs_;
-  std::array<std::size_t, MAXDIM> strides_;
+  IArray<D> nxs_;
+  IArray<D> strides_;
   int rank_;
-};
-
-// copy constructor (does a shallow copy)
-
-template <typename T>
-PortableMDArray<T>::PortableMDArray(const PortableMDArray<T> &src) noexcept {
-  nxs_ = src.nxs_;
-  rank_ = src.rank_;
-  strides_ = src.strides_;
-  if (src.pdata_) pdata_ = src.pdata_;
-}
-
-// shallow copy assignment operator
-
-template <typename T>
-PortableMDArray<T> &
-PortableMDArray<T>::operator=(const PortableMDArray<T> &src) noexcept {
-  if (this != &src) {
-    nxs_ = src.nxs;
-    rank_ = src.rank_;
-    pdata_ = src.pdata_;
-  }
-  return *this;
-}
-
-// Checks that arrays point to same data with same shape
-// note this POINTER equivalence, not data equivalence
-template <typename T>
-bool PortableMDArray<T>::operator==(const PortableMDArray<T> &rhs) const {
-  return (pdata_ == rhs.pdata_ && nxs_ == rhs.nxs_); // NB rank is implied
-}
 
+ public:
+};
 //----------------------------------------------------------------------------------------
 //! \fn PortableMDArray::InitWithShallowSlice()
 //  \brief shallow copy of nvar elements in dimension dim of an array,
@@ -288,36 +301,10 @@ bool PortableMDArray<T>::operator==(const PortableMDArray<T> &rhs) const {
 //  entries of the src array for d<dim (cannot access any nx4=2, etc. entries
 //  if dim=3 for example)
 
-template <typename T>
-PORTABLE_FUNCTION void
-PortableMDArray<T>::InitWithShallowSlice(const PortableMDArray<T> &src, const int dim,
-                                         const int indx, const int nvar) {
-  pdata_ = src.pdata_;
-  std::size_t offs = indx;
-  nxs_[dim - 1] = nvar;
-
-  for (std::size_t i = 0; i < dim - 1; ++i) {
-    nxs_[i] = src.nxs_[i];
-    offs *= nxs_[i];
-  }
-  for (std::size_t i = dim; i < MAXDIM; ++i) {
-    nxs_[i] = 1;
-  }
-  pdata_ += offs;
-
-  return;
-}
-
 //----------------------------------------------------------------------------------------
 //! \fn PortableMDArray::SwapPortableMDArray()
 //  \brief  swap pdata_ pointers of two equally sized PortableMDArrays
 //  (shallow swap)
 // Does not allocate memory for either array
 
-template <typename T>
-void PortableMDArray<T>::SwapPortableMDArray(PortableMDArray<T> &array2) {
-  std::swap(pdata_, array2.pdata_);
-  return;
-}
-
 #endif // _PORTABLE_ARRAYS_HPP_
diff --git a/ports-of-call/utility/array_algo.hpp b/ports-of-call/utility/array_algo.hpp
index e4c52944..f10d09c1 100644
--- a/ports-of-call/utility/array_algo.hpp
+++ b/ports-of-call/utility/array_algo.hpp
@@ -7,65 +7,88 @@
 
 namespace util {
 
+template <auto N>
+constexpr auto
+is(std::integral_constant<decltype(N), N>) { // = std::make_index_sequence<N>;
+  return std::make_index_sequence<N>{};
+}
+
+template <class A>
+using value_t = typename A::value_type;
+
+template <class A, class Op>
+using reduction_value_t =
+    decltype(std::declval<Op>()(std::declval<value_t<A>>(), std::declval<value_t<A>>()));
+
 namespace detail {
-template <class T, std::size_t N, class F, std::size_t... Is>
-PORTABLE_FORCEINLINE_FUNCTION constexpr auto
-array_map_impl(std::array<T, N> const &x, F f, std::index_sequence<Is...>) {
+template <class A, class F, std::size_t... Is>
+PORTABLE_FORCEINLINE_FUNCTION constexpr auto array_map_impl(A const &x, F f,
+                                                            std::index_sequence<Is...>) {
   return std::array{f(x[Is])...};
 }
 
-template <class T, class U, std::size_t N, class F, std::size_t... Is>
-PORTABLE_FORCEINLINE_FUNCTION constexpr auto
-array_map_impl(std::array<T, N> const &x, std::array<U, N> const &y, F f,
-               std::index_sequence<Is...>) {
-  return std::array{f(x[Is], y[Is])...};
+template <class A, class B, std::size_t N, class F, std::size_t... Is>
+PORTABLE_FORCEINLINE_FUNCTION constexpr auto array_map_impl(A const &x, B const &y, F f,
+                                                            std::index_sequence<Is...>) {
+  return A{f(x[Is], y[Is])...};
 }
 
-template <std::size_t f, std::size_t l, class T, std::size_t N, class Op>
-PORTABLE_INLINE_FUNCTION constexpr T array_reduce_impl(std::array<T, N> const &x, Op op) {
+template <std::size_t f, std::size_t l, class A, class Op>
+PORTABLE_INLINE_FUNCTION constexpr auto array_reduce_impl(A const &x, Op op) {
   if constexpr ((l - f) == 1)
     return x[f];
   else {
     constexpr std::size_t n = l - f;
-    T left_sum = array_reduce_impl<f, f + n / 2>(x, op);
-    T right_sum = array_reduce_impl<f + n / 2, l>(x, op);
+    auto left_sum = array_reduce_impl<f, f + n / 2>(x, op);
+    auto right_sum = array_reduce_impl<f + n / 2, l>(x, op);
     return op(left_sum, right_sum);
   }
 }
 
 } // namespace detail
 
-template <class T, std::size_t N, class F>
-PORTABLE_FORCEINLINE_FUNCTION constexpr auto array_map(std::array<T, N> const &x, F f) {
-  return detail::array_map_impl(x, f, std::make_index_sequence<N>{});
+template <class A, class F>
+PORTABLE_FORCEINLINE_FUNCTION constexpr auto array_map(A const &x, F f) {
+  return detail::array_map_impl(x, f, is(x.size()));
 }
 
-template <class T, class U, std::size_t N, class F>
-PORTABLE_FORCEINLINE_FUNCTION constexpr auto array_map(std::array<T, N> const &x,
-                                                       std::array<U, N> const &y, F f) {
-  return detail::array_map_impl(x, y, f, std::make_index_sequence<N>{});
+template <class A, class B, class F>
+PORTABLE_FORCEINLINE_FUNCTION constexpr auto array_map(A const &x, B const &y, F f) {
+  return detail::array_map_impl(x, y, f, is(x.size()));
 }
 
-template <std::size_t I, class T, std::size_t N, class Op>
-PORTABLE_FORCEINLINE_FUNCTION constexpr T array_partial_reduce(std::array<T, N> x,
-                                                               T initial_value, Op op) {
-  static_assert(I <= N);
+template <std::size_t I, class A, class Op, class T = reduction_value_t<A, Op>>
+PORTABLE_FORCEINLINE_FUNCTION constexpr T array_partial_reduce(A x, T initial_value,
+                                                               Op op) {
+  static_assert(I <= x.size());
   if constexpr (I == 0)
     return initial_value;
   else
     return detail::array_reduce_impl<0, I>(x, op);
 }
 
-template <class T, std::size_t N, class Op>
-PORTABLE_FORCEINLINE_FUNCTION constexpr T array_reduce(std::array<T, N> x,
-                                                       T initial_value, Op op) {
-  return array_partial_reduce<N>(x, initial_value, op);
+template <class A, class Op, class T = reduction_value_t<A, Op>>
+
+PORTABLE_FORCEINLINE_FUNCTION constexpr T array_reduce(A x, T initial_value, Op op) {
+  return array_partial_reduce<x.size()>(x, initial_value, op);
+}
+/*
+namespace detail {
+template <class A, std::size_t... I>
+PORTABLE_FORCEINLINE_FUNCTION constexpr bool
+arrays_equal_impl(const A &x, const A &y, std::index_sequence<I...>) {
+  return (... && (x[I] == y[I]));
 }
-template <std::size_t... Is>
+} // namespace detail
+template <class A>
+PORTABLE_FORCEINLINE_FUNCTION constexpr bool arrays_equal(const A &x, const A &y) {
+  return detail::arrays_equal_impl(x, y, is(x.size()));
+}*/
+/*template <std::size_t... Is>
 PORTABLE_FORCEINLINE_FUNCTION constexpr auto
 as_array(std::index_sequence<sizeof...(Is)>) {
   return std::array{Is...};
-}
+}*/
 
 } // namespace util
 
diff --git a/ports-of-call/utility/index_algo.hpp b/ports-of-call/utility/index_algo.hpp
index 31a5d8e8..efc4b3b3 100644
--- a/ports-of-call/utility/index_algo.hpp
+++ b/ports-of-call/utility/index_algo.hpp
@@ -2,14 +2,15 @@
 #define _PORTSOFCALL_UTILITY_INDEX_ALGO_HPP_
 
 #include "../portability.hpp"
-#include "array.hpp"
+#include "array_algo.hpp"
 #include <array>
 #include <numeric>
 #include <type_traits>
 
 namespace util {
 
-template <auto I, class T, std::size_t N>
+/*
+ template <auto I, class T, std::size_t N>
 PORTABLE_FORCEINLINE_FUNCTION static constexpr auto
 get_stride(std::array<T, N> const &dim) {
   static_assert(I < dim.size(), "Dim index is out of bounds");
@@ -48,48 +49,41 @@ PORTABLE_FORCEINLINE_FUNCTION static constexpr auto
 findex(std::array<T, N> const &ijk, std::array<T, N> const &dim) {
   return fast_findex(ijk, dim, get_strides(dim));
 }
-
-namespace handroll {
-template <auto I, class T, std::size_t N>
-PORTABLE_FORCEINLINE_FUNCTION static constexpr auto
-get_stride(std::array<T, N> const &dim) {
+*/
+template <auto I, class A>
+PORTABLE_FORCEINLINE_FUNCTION static constexpr auto get_stride(A const &dim) {
 
   // column major
-  return array_partial_reduce<I>(dim, T{1}, std::multiplies<std::size_t>{});
+  return array_partial_reduce<I>(dim, value_t<A>{1}, std::multiplies<std::size_t>{});
 }
 
 namespace detail {
-template <class T, std::size_t N, std::size_t... Is>
+template <class A, std::size_t... Is>
 PORTABLE_FORCEINLINE_FUNCTION static constexpr auto
-get_strides_impl(std::array<T, N> const &dim, std::index_sequence<Is...>) {
-  return std::array{get_stride<Is>(dim)...};
+get_strides_impl(A const &dim, std::index_sequence<Is...>) {
+  return A{get_stride<Is>(dim)...};
 }
-template <class T, std::size_t N>
-PORTABLE_FORCEINLINE_FUNCTION static constexpr auto
-get_strides(std::array<T, N> const &dim) {
-  return detail::get_strides_impl(dim, std::make_index_sequence<dim.size()>{});
-} // namespace handroll::detail
-template <class T, std::size_t N>
+template <class A>
+PORTABLE_FORCEINLINE_FUNCTION static constexpr auto get_strides(A const &dim) {
+  return detail::get_strides_impl(dim, is(dim.size()));
+}
+} // namespace detail
+template <class A>
 PORTABLE_FORCEINLINE_FUNCTION static constexpr auto
-fast_findex(std::array<T, N> const &ijk, std::array<T, N> const &dim,
-            std::array<T, N> const &stride) {
+fast_findex(A const &ijk, A const &dim, A const &stride) {
   // TODO: assert ijk in bounds
-  return array_reduce(
-      array_map(ijk, stride, [](auto a, auto b) { return a * b; }), T{1},
-      std::plus<std::size_t>{});
+  return array_reduce(array_map(ijk, stride, [](auto a, auto b) { return a * b; }),
+                      value_t<A>{1}, std::plus<std::size_t>{});
 }
-template <class T, std::size_t N>
-PORTABLE_FORCEINLINE_FUNCTION static constexpr auto
-findex(std::array<T, N> const &ijk, std::array<T, N> const &dim) {
+template <class A>
+PORTABLE_FORCEINLINE_FUNCTION static constexpr auto findex(A const &ijk, A const &dim) {
   return fast_findex(ijk, dim, get_strides(dim));
 }
-} // namespace handroll
 
-template <class T, std::size_t N>
-PORTABLE_FORCEINLINE_FUNCTION static constexpr std::array<T, N>
-fast_mindices(std::size_t idx, std::array<T, N> const &dim,
-              std::array<T, N> const &stride) {
-  std::array<T, N> mdidx;
+template <class A>
+PORTABLE_FORCEINLINE_FUNCTION static constexpr A
+fast_mindices(std::size_t idx, A const &dim, A const &stride) {
+  A mdidx;
   for (std::int64_t i = dim.size() - 1; i >= 0; --i) {
     mdidx[i] = idx / std::size_t(stride[i]);
     idx -= mdidx[i] * std::size_t(stride[i]);
@@ -97,9 +91,8 @@ fast_mindices(std::size_t idx, std::array<T, N> const &dim,
   return mdidx;
 }
 
-template <class Array>
-PORTABLE_FORCEINLINE_FUNCTION static constexpr auto mindices(std::size_t idx,
-                                                             Array dim) {
+template <class A>
+PORTABLE_FORCEINLINE_FUNCTION static constexpr auto mindices(std::size_t idx, A dim) {
   return fast_mindices(idx, dim, get_strides(dim));
 }
 
diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
index 65d5d205..00e1ed3a 100644
--- a/test/CMakeLists.txt
+++ b/test/CMakeLists.txt
@@ -47,4 +47,5 @@ target_link_libraries(test_portsofcall PRIVATE ports-of-call::ports-of-call
 include(Catch)
 catch_discover_tests(test_portsofcall)
 
-target_sources(test_portsofcall PRIVATE test_portability.cpp test_array.cpp)
+target_sources(test_portsofcall PRIVATE test_portability.cpp test_array.cpp
+                                        test_mdidx.cpp)
diff --git a/test/test_mdidx.cpp b/test/test_mdidx.cpp
new file mode 100644
index 00000000..dbf9b964
--- /dev/null
+++ b/test/test_mdidx.cpp
@@ -0,0 +1,144 @@
+// © (or copyright) 2019-2024. Triad National Security, LLC. All rights
+// reserved.  This program was produced under U.S. Government contract
+// 89233218CNA000001 for Los Alamos National Laboratory (LANL), which is
+// operated by Triad National Security, LLC for the U.S.  Department of
+// Energy/National Nuclear Security Administration. All rights in the
+// program are reserved by Triad National Security, LLC, and the
+// U.S. Department of Energy/National Nuclear Security
+// Administration. The Government is granted for itself and others acting
+// on its behalf a nonexclusive, paid-up, irrevocable worldwide license
+// in this material to reproduce, prepare derivative works, distribute
+// copies to the public, perform publicly and display publicly, and to
+// permit others to do so.
+
+#include <ports-of-call/portability.hpp>
+#include <ports-of-call/portable_arrays.hpp>
+#include <utility>
+#include <vector>
+
+#include "test_utilities.hpp"
+
+#ifndef CATCH_CONFIG_FAST_COMPILE
+#define CATCH_CONFIG_FAST_COMPILE
+#include <catch2/benchmark/catch_benchmark.hpp>
+#include <catch2/catch_test_macros.hpp>
+#include <catch2/matchers/catch_matchers.hpp>
+#include <catch2/matchers/catch_matchers_floating_point.hpp>
+#endif
+
+TEST_CASE("PortableMDArrays Sizes are sane", "[PortableMDArray]") {
+  using tape_t = std::vector<int>;
+
+  constexpr std::array<std::size_t, MAXDIM> nxs{2, 5, 10, 10, 5, 2};
+
+  std::vector<int> subsz(MAXDIM);
+  std::vector<tape_t> dats;
+
+  std::partial_sum(nxs.cbegin(), nxs.cend(), subsz.begin(), std::multiplies<int>());
+  for (auto i = 0; i < MAXDIM; ++i) {
+    dats.push_back(std::vector<int>(subsz[i], nxs[i]));
+  }
+
+  SECTION("Can construct correctly") {
+    auto pmd1 = PortableMDArray<int>(dats[0].data(), nxs[0]);
+
+    REQUIRE(pmd1.GetSize() == subsz[0]);
+    REQUIRE(pmd1.GetRank() == 1);
+
+    auto pmd3 = PortableMDArray<int>(dats[3].data(), nxs[0], nxs[1], nxs[2]);
+
+    REQUIRE(pmd3.GetSize() == subsz[2]);
+    REQUIRE(pmd3.GetRank() == 3);
+    REQUIRE(pmd3.GetDim1() == nxs[2]);
+    REQUIRE(pmd3.GetDim2() == nxs[1]);
+    REQUIRE(pmd3.GetDim3() == nxs[0]);
+  }
+}
+
+TEST_CASE("Correct portable indexing", "[PortableMDArray]") {
+  // layout
+  auto iflat = [](auto nx, auto nxny) {
+    return [=](auto k, auto j, auto i) { return i + nx * j + nxny * k; };
+  };
+
+  constexpr std::size_t NX = 32, NY = 64, NZ = 4;
+  constexpr std::size_t NC = NX * NY * NZ;
+  constexpr Real scale = 0.1;
+  // size in bytes
+  constexpr const std::size_t NCb = NC * sizeof(Real);
+
+  // vector length N on host of Real
+  std::vector<Real> tape_ref(NC), tape_buf(NC);
+
+  for (auto n = 0; n < NC; ++n) {
+    tape_ref[n] = scale * static_cast<Real>(n);
+  }
+
+  // device pointer
+  Real *tape_d = (Real *)PORTABLE_MALLOC(NCb);
+
+  auto view_d = PortableMDArray<Real>(tape_d, NZ, NY, NX);
+
+  // set device values
+  portableFor(
+      "set unique val", 0, NZ, 0, NY, 0, NX,
+      PORTABLE_LAMBDA(const int &k, const int &j, const int &i) {
+        view_d(k, j, i) = scale * iflat(NX, NX * NY)(k, j, i);
+      });
+
+  portableCopyToHost(tape_buf.data(), tape_d, NCb);
+
+  for (auto n = 0; n < NC; ++n) {
+    INFO("REF=" << tape_ref[n] << " BUF=" << tape_buf[n] << " n=" << n);
+    REQUIRE_THAT(tape_buf[n], Catch::Matchers::WithinRel(tape_ref[n]));
+  }
+
+  PORTABLE_FREE(tape_d);
+}
+
+// runs benchmarks for indexing.
+// note: to run, execute
+//    ./tests/test_portsofcall "[!benchmark]"
+// (you may need to escape the `!` char depending
+// on your shell)
+TEST_CASE("Benchmark Indexing", "[!benchmark]") {
+  constexpr std::size_t NX = 32, NY = 64, NZ = 16;
+  constexpr std::size_t NU = 256, NV = 256, NW = 128;
+
+  // construct a series of PortableMDArrays
+  // of different rank and sizes, and
+  // iterates through all contiguous indexing.
+
+  // clang-format off
+  SECTION("small") {
+    BENCHMARK("index calc 1D") 
+    { 
+      return testing::idx_contiguous_bm(std::array{NX}); 
+    };
+    BENCHMARK("index calc 2D") 
+    { 
+      return testing::idx_contiguous_bm(std::array{NX, NY}); 
+    };
+    BENCHMARK("index calc 3D") 
+    { 
+      return testing::idx_contiguous_bm(std::array{NX, NY, NZ}); 
+    };
+  }
+
+  SECTION("big") {
+    BENCHMARK("index calc 1D") 
+    { 
+      return testing::idx_contiguous_bm(std::array{NU}); 
+    };
+    BENCHMARK("index calc 2D") 
+    { 
+      return testing::idx_contiguous_bm(std::array{NU, NV}); 
+    };
+    BENCHMARK("index calc 3D") 
+    { 
+      return testing::idx_contiguous_bm(std::array{NU, NV, NW}); 
+    };
+
+  }
+  // clang-format on
+}
diff --git a/test/test_utilities.hpp b/test/test_utilities.hpp
index 2d5ad3fb..0cb2de41 100644
--- a/test/test_utilities.hpp
+++ b/test/test_utilities.hpp
@@ -1,9 +1,8 @@
 #ifndef _TEST_PORTOFCALL_HPP_
 #define _TEST_PORTOFCALL_HPP_
 
-#include <__utility/integer_sequence.h>
-#include <array>
 #include <functional>
+#include <ports-of-call/array.hpp>
 #include <ports-of-call/portability.hpp>
 #include <ports-of-call/portable_arrays.hpp>
 #include <utility>
@@ -27,10 +26,8 @@ constexpr auto array_fix_impl(const Array &a, std::index_sequence<I...> is) {
 // it into an array 2N {0, NX} contiguous pairs.
 // ex:
 //    array{10,12,13} -> array{0,10,0,12,0,13}
-template <typename T, std::size_t N,
-          typename Indx = std::make_index_sequence<N>>
-PORTABLE_FORCEINLINE_FUNCTION constexpr auto
-array_fix(const std::array<T, N> &a) {
+template <typename T, std::size_t N, typename Indx = std::make_index_sequence<N>>
+PORTABLE_FORCEINLINE_FUNCTION constexpr auto array_fix(const std::array<T, N> &a) {
   return array_fix_impl(a, Indx{});
 }
 
@@ -50,13 +47,12 @@ PORTABLE_INLINE_FUNCTION constexpr auto pf_invoke(View v, const Array &ext_par,
 
   auto stud_arr = array_fix(ext_par);
   portableFor(
-      "set unique val", stud_arr[J]...,
-      PORTABLE_LAMBDA(auto... is) { v(is...) = 1.0; });
+      "set unique val", stud_arr[J]..., PORTABLE_LAMBDA(auto... is) { v(is...) = 1.0; });
 }
 
 template <class Ptr, class Array, std::size_t... I>
-PORTABLE_FORCEINLINE_FUNCTION constexpr decltype(auto)
-mdview(Ptr *d, const Array &arr, std::index_sequence<I...>) {
+PORTABLE_FORCEINLINE_FUNCTION constexpr decltype(auto) mdview(Ptr *d, const Array &arr,
+                                                              std::index_sequence<I...>) {
   return PortableMDArray<Real>(d, arr[I]...);
 }
 
@@ -65,8 +61,7 @@ mdview(Ptr *d, const Array &arr, std::index_sequence<I...>) {
 template <class T, std::size_t N, typename I1 = std::make_index_sequence<N>,
           typename I2 = std::make_index_sequence<2 * N>>
 PORTABLE_FUNCTION auto idx_contiguous_bm(const std::array<T, N> &nxa) {
-  auto nc =
-      std::accumulate(std::begin(nxa), std::end(nxa), 1, std::multiplies<T>{});
+  auto nc = std::accumulate(std::begin(nxa), std::end(nxa), 1, std::multiplies<T>{});
 
   Real *tape_d = (Real *)PORTABLE_MALLOC(nc * sizeof(Real));
 

From ece0a3cf7d0307d76ddb34acfb4b984cca5896c7 Mon Sep 17 00:00:00 2001
From: Christopher Mauney <mauneyc@pn2305185.lanl.gov>
Date: Fri, 30 Aug 2024 13:49:37 -0600
Subject: [PATCH 10/22] comment cleanup and copyright

---
 ports-of-call/array.hpp              | 13 +++++++
 ports-of-call/utility/array_algo.hpp | 13 +++++++
 ports-of-call/utility/index_algo.hpp | 54 +++++++---------------------
 test/test_array.cpp                  | 13 +++++++
 test/test_utilities.hpp              | 17 +++++++--
 5 files changed, 67 insertions(+), 43 deletions(-)

diff --git a/ports-of-call/array.hpp b/ports-of-call/array.hpp
index 64569f01..a5556e4a 100644
--- a/ports-of-call/array.hpp
+++ b/ports-of-call/array.hpp
@@ -1,3 +1,16 @@
+// © (or copyright) 2019-2024. Triad National Security, LLC. All rights
+// reserved.  This program was produced under U.S. Government contract
+// 89233218CNA000001 for Los Alamos National Laboratory (LANL), which is
+// operated by Triad National Security, LLC for the U.S.  Department of
+// Energy/National Nuclear Security Administration. All rights in the
+// program are reserved by Triad National Security, LLC, and the
+// U.S. Department of Energy/National Nuclear Security
+// Administration. The Government is granted for itself and others acting
+// on its behalf a nonexclusive, paid-up, irrevocable worldwide license
+// in this material to reproduce, prepare derivative works, distribute
+// copies to the public, perform publicly and display publicly, and to
+// permit others to do so.
+
 #ifndef _PORTS_OF_CALL_ARRAY_HPP_
 #define _PORTS_OF_CALL_ARRAY_HPP_
 
diff --git a/ports-of-call/utility/array_algo.hpp b/ports-of-call/utility/array_algo.hpp
index f10d09c1..bf15786c 100644
--- a/ports-of-call/utility/array_algo.hpp
+++ b/ports-of-call/utility/array_algo.hpp
@@ -1,3 +1,16 @@
+// © (or copyright) 2019-2024. Triad National Security, LLC. All rights
+// reserved.  This program was produced under U.S. Government contract
+// 89233218CNA000001 for Los Alamos National Laboratory (LANL), which is
+// operated by Triad National Security, LLC for the U.S.  Department of
+// Energy/National Nuclear Security Administration. All rights in the
+// program are reserved by Triad National Security, LLC, and the
+// U.S. Department of Energy/National Nuclear Security
+// Administration. The Government is granted for itself and others acting
+// on its behalf a nonexclusive, paid-up, irrevocable worldwide license
+// in this material to reproduce, prepare derivative works, distribute
+// copies to the public, perform publicly and display publicly, and to
+// permit others to do so.
+
 #ifndef _PORTSOFCALL_UTILITY_ARRAY_ALGO_HPP_
 #define _PORTSOFCALL_UTILITY_ARRAY_ALGO_HPP_
 
diff --git a/ports-of-call/utility/index_algo.hpp b/ports-of-call/utility/index_algo.hpp
index efc4b3b3..3a507620 100644
--- a/ports-of-call/utility/index_algo.hpp
+++ b/ports-of-call/utility/index_algo.hpp
@@ -1,3 +1,16 @@
+// © (or copyright) 2019-2024. Triad National Security, LLC. All rights
+// reserved.  This program was produced under U.S. Government contract
+// 89233218CNA000001 for Los Alamos National Laboratory (LANL), which is
+// operated by Triad National Security, LLC for the U.S.  Department of
+// Energy/National Nuclear Security Administration. All rights in the
+// program are reserved by Triad National Security, LLC, and the
+// U.S. Department of Energy/National Nuclear Security
+// Administration. The Government is granted for itself and others acting
+// on its behalf a nonexclusive, paid-up, irrevocable worldwide license
+// in this material to reproduce, prepare derivative works, distribute
+// copies to the public, perform publicly and display publicly, and to
+// permit others to do so.
+
 #ifndef _PORTSOFCALL_UTILITY_INDEX_ALGO_HPP_
 #define _PORTSOFCALL_UTILITY_INDEX_ALGO_HPP_
 
@@ -9,47 +22,6 @@
 
 namespace util {
 
-/*
- template <auto I, class T, std::size_t N>
-PORTABLE_FORCEINLINE_FUNCTION static constexpr auto
-get_stride(std::array<T, N> const &dim) {
-  static_assert(I < dim.size(), "Dim index is out of bounds");
-
-  // column major
-  return std::accumulate(std::begin(dim), std::begin(dim) + I, std::size_t{1},
-                         std::multiplies<std::size_t>{});
-}
-
-namespace detail {
-template <class T, std::size_t N, std::size_t... Is>
-PORTABLE_FORCEINLINE_FUNCTION static constexpr auto
-get_strides_impl(std::array<T, N> const &dim, std::index_sequence<Is...>) {
-  return std::array{get_stride<Is>(dim)...};
-}
-
-} // namespace detail
-
-template <class T, std::size_t N>
-PORTABLE_FORCEINLINE_FUNCTION static constexpr auto
-get_strides(std::array<T, N> const &dim) {
-  return detail::get_strides_impl(dim, std::make_index_sequence<dim.size()>{});
-}
-
-template <class T, std::size_t N>
-PORTABLE_FORCEINLINE_FUNCTION static constexpr auto
-fast_findex(std::array<T, N> const &ijk, std::array<T, N> const &dim,
-            std::array<T, N> const &stride) {
-  // TODO: assert ijk in bounds
-  return std::inner_product(std::begin(ijk), std::end(ijk), std::begin(stride),
-                            std::size_t{0});
-}
-
-template <class T, std::size_t N>
-PORTABLE_FORCEINLINE_FUNCTION static constexpr auto
-findex(std::array<T, N> const &ijk, std::array<T, N> const &dim) {
-  return fast_findex(ijk, dim, get_strides(dim));
-}
-*/
 template <auto I, class A>
 PORTABLE_FORCEINLINE_FUNCTION static constexpr auto get_stride(A const &dim) {
 
diff --git a/test/test_array.cpp b/test/test_array.cpp
index 3a9423fc..64f0e83b 100644
--- a/test/test_array.cpp
+++ b/test/test_array.cpp
@@ -1,3 +1,16 @@
+// © (or copyright) 2019-2024. Triad National Security, LLC. All rights
+// reserved.  This program was produced under U.S. Government contract
+// 89233218CNA000001 for Los Alamos National Laboratory (LANL), which is
+// operated by Triad National Security, LLC for the U.S.  Department of
+// Energy/National Nuclear Security Administration. All rights in the
+// program are reserved by Triad National Security, LLC, and the
+// U.S. Department of Energy/National Nuclear Security
+// Administration. The Government is granted for itself and others acting
+// on its behalf a nonexclusive, paid-up, irrevocable worldwide license
+// in this material to reproduce, prepare derivative works, distribute
+// copies to the public, perform publicly and display publicly, and to
+// permit others to do so.
+
 #include "ports-of-call/array.hpp"
 #include "ports-of-call/portability.hpp"
 
diff --git a/test/test_utilities.hpp b/test/test_utilities.hpp
index 0cb2de41..30ba2090 100644
--- a/test/test_utilities.hpp
+++ b/test/test_utilities.hpp
@@ -1,5 +1,18 @@
-#ifndef _TEST_PORTOFCALL_HPP_
-#define _TEST_PORTOFCALL_HPP_
+// © (or copyright) 2019-2024. Triad National Security, LLC. All rights
+// reserved.  This program was produced under U.S. Government contract
+// 89233218CNA000001 for Los Alamos National Laboratory (LANL), which is
+// operated by Triad National Security, LLC for the U.S.  Department of
+// Energy/National Nuclear Security Administration. All rights in the
+// program are reserved by Triad National Security, LLC, and the
+// U.S. Department of Energy/National Nuclear Security
+// Administration. The Government is granted for itself and others acting
+// on its behalf a nonexclusive, paid-up, irrevocable worldwide license
+// in this material to reproduce, prepare derivative works, distribute
+// copies to the public, perform publicly and display publicly, and to
+// permit others to do so.
+
+#ifndef _TEST_TEST_UTILITIES_HPP_
+#define _TEST_TEST_UTILITIES_HPP_
 
 #include <functional>
 #include <ports-of-call/array.hpp>

From 66685871b1be63addffd095f6d921effb884d3ea Mon Sep 17 00:00:00 2001
From: Christopher Mauney <mauneyc@pn2305185.lanl.gov>
Date: Fri, 30 Aug 2024 14:07:07 -0600
Subject: [PATCH 11/22] short comments in algos

---
 ports-of-call/portable_arrays.hpp    | 31 +++++++++++++---------------
 ports-of-call/utility/array_algo.hpp | 30 ++++++++++-----------------
 ports-of-call/utility/index_algo.hpp |  9 +++++++-
 3 files changed, 33 insertions(+), 37 deletions(-)

diff --git a/ports-of-call/portable_arrays.hpp b/ports-of-call/portable_arrays.hpp
index 72ef43f6..0938492a 100644
--- a/ports-of-call/portable_arrays.hpp
+++ b/ports-of-call/portable_arrays.hpp
@@ -97,9 +97,11 @@ class PortableMDArray {
     pdata_ = data;
     update_layout(nxs...);
   }
-  // public function to swap underlying data pointers of two equally-sized
-  // arrays
-  // void SwapPortableMDArray(PortableMDArray<T> &array2);
+  //----------------------------------------------------------------------------------------
+  //! \fn PortableMDArray::SwapPortableMDArray()
+  //  \brief  swap pdata_ pointers of two equally sized PortableMDArrays
+  //  (shallow swap)
+  // Does not allocate memory for either array
 
   void SwapPortableMDArray(PortableMDArray<T, D> &array2) {
     std::swap(pdata_, array2.pdata_);
@@ -209,6 +211,15 @@ class PortableMDArray {
 
   bool operator!=(const PortableMDArray<T, D> &other) const { return !(*this == other); }
 
+  //----------------------------------------------------------------------------------------
+  //! \fn PortableMDArray::InitWithShallowSlice()
+  //  \brief shallow copy of nvar elements in dimension dim of an array,
+  //  starting at index=indx. Copies pointer to data, but not data itself.
+
+  //  Shallow slice is only able to address the "nvar" range in "dim", and all
+  //  entries of the src array for d<dim (cannot access any nx4=2, etc. entries
+  //  if dim=3 for example)
+
   PORTABLE_FUNCTION void InitWithShallowSlice(const PortableMDArray<T, D> &src,
                                               const int dim, const int indx,
                                               const int nvar) {
@@ -292,19 +303,5 @@ class PortableMDArray {
 
  public:
 };
-//----------------------------------------------------------------------------------------
-//! \fn PortableMDArray::InitWithShallowSlice()
-//  \brief shallow copy of nvar elements in dimension dim of an array,
-//  starting at index=indx. Copies pointer to data, but not data itself.
-
-//  Shallow slice is only able to address the "nvar" range in "dim", and all
-//  entries of the src array for d<dim (cannot access any nx4=2, etc. entries
-//  if dim=3 for example)
-
-//----------------------------------------------------------------------------------------
-//! \fn PortableMDArray::SwapPortableMDArray()
-//  \brief  swap pdata_ pointers of two equally sized PortableMDArrays
-//  (shallow swap)
-// Does not allocate memory for either array
 
 #endif // _PORTABLE_ARRAYS_HPP_
diff --git a/ports-of-call/utility/array_algo.hpp b/ports-of-call/utility/array_algo.hpp
index bf15786c..3e49a5c7 100644
--- a/ports-of-call/utility/array_algo.hpp
+++ b/ports-of-call/utility/array_algo.hpp
@@ -20,15 +20,17 @@
 
 namespace util {
 
+// shorter make_index_sequence
 template <auto N>
-constexpr auto
-is(std::integral_constant<decltype(N), N>) { // = std::make_index_sequence<N>;
+constexpr auto is(std::integral_constant<decltype(N), N>) {
   return std::make_index_sequence<N>{};
 }
 
+// shorter value_type
 template <class A>
 using value_t = typename A::value_type;
 
+// determines the return type of Op(a,b)
 template <class A, class Op>
 using reduction_value_t =
     decltype(std::declval<Op>()(std::declval<value_t<A>>(), std::declval<value_t<A>>()));
@@ -60,11 +62,16 @@ PORTABLE_INLINE_FUNCTION constexpr auto array_reduce_impl(A const &x, Op op) {
 
 } // namespace detail
 
+// maps an unary function f(x) to each array value, returning an array of results
+// x = {f(a[0]), f(a[1]),..., f(a[N-1])}
 template <class A, class F>
 PORTABLE_FORCEINLINE_FUNCTION constexpr auto array_map(A const &x, F f) {
   return detail::array_map_impl(x, f, is(x.size()));
 }
 
+// maps a binary function to each array value, returning an array of results
+// x = {f(a[0], b[0]), f(a[1], b[1]),..., f(a[N-1], b[N-1])}
+
 template <class A, class B, class F>
 PORTABLE_FORCEINLINE_FUNCTION constexpr auto array_map(A const &x, B const &y, F f) {
   return detail::array_map_impl(x, y, f, is(x.size()));
@@ -82,26 +89,11 @@ PORTABLE_FORCEINLINE_FUNCTION constexpr T array_partial_reduce(A x, T initial_va
 
 template <class A, class Op, class T = reduction_value_t<A, Op>>
 
+// performs a reduction on an array of values
+// e.g. x = sum_i a[i]
 PORTABLE_FORCEINLINE_FUNCTION constexpr T array_reduce(A x, T initial_value, Op op) {
   return array_partial_reduce<x.size()>(x, initial_value, op);
 }
-/*
-namespace detail {
-template <class A, std::size_t... I>
-PORTABLE_FORCEINLINE_FUNCTION constexpr bool
-arrays_equal_impl(const A &x, const A &y, std::index_sequence<I...>) {
-  return (... && (x[I] == y[I]));
-}
-} // namespace detail
-template <class A>
-PORTABLE_FORCEINLINE_FUNCTION constexpr bool arrays_equal(const A &x, const A &y) {
-  return detail::arrays_equal_impl(x, y, is(x.size()));
-}*/
-/*template <std::size_t... Is>
-PORTABLE_FORCEINLINE_FUNCTION constexpr auto
-as_array(std::index_sequence<sizeof...(Is)>) {
-  return std::array{Is...};
-}*/
 
 } // namespace util
 
diff --git a/ports-of-call/utility/index_algo.hpp b/ports-of-call/utility/index_algo.hpp
index 3a507620..168676b9 100644
--- a/ports-of-call/utility/index_algo.hpp
+++ b/ports-of-call/utility/index_algo.hpp
@@ -22,10 +22,11 @@
 
 namespace util {
 
+// returns the stride (offset) needed to move one unit in dimension D
 template <auto I, class A>
 PORTABLE_FORCEINLINE_FUNCTION static constexpr auto get_stride(A const &dim) {
 
-  // column major
+  // NB: column major
   return array_partial_reduce<I>(dim, value_t<A>{1}, std::multiplies<std::size_t>{});
 }
 
@@ -40,6 +41,9 @@ PORTABLE_FORCEINLINE_FUNCTION static constexpr auto get_strides(A const &dim) {
   return detail::get_strides_impl(dim, is(dim.size()));
 }
 } // namespace detail
+
+// returns the flat (1D) index of the md index set {i,j,k}
+// NB: fast because the strides are provided and don't need to be recomputed
 template <class A>
 PORTABLE_FORCEINLINE_FUNCTION static constexpr auto
 fast_findex(A const &ijk, A const &dim, A const &stride) {
@@ -47,11 +51,14 @@ fast_findex(A const &ijk, A const &dim, A const &stride) {
   return array_reduce(array_map(ijk, stride, [](auto a, auto b) { return a * b; }),
                       value_t<A>{1}, std::plus<std::size_t>{});
 }
+
+// same as fast_findex, except the strides are calculated on the fly
 template <class A>
 PORTABLE_FORCEINLINE_FUNCTION static constexpr auto findex(A const &ijk, A const &dim) {
   return fast_findex(ijk, dim, get_strides(dim));
 }
 
+// returns the md index set {i,j,k} of the flat (1D) index
 template <class A>
 PORTABLE_FORCEINLINE_FUNCTION static constexpr A
 fast_mindices(std::size_t idx, A const &dim, A const &stride) {

From 825301a35f1b6dd8184155041ce88d5dc015d98d Mon Sep 17 00:00:00 2001
From: Christopher Mauney <mauneyc@pn2305185.lanl.gov>
Date: Fri, 30 Aug 2024 14:23:19 -0600
Subject: [PATCH 12/22] removed unneeded parameter

---
 ports-of-call/utility/array_algo.hpp | 25 ++++++++++++++++++++++++-
 1 file changed, 24 insertions(+), 1 deletion(-)

diff --git a/ports-of-call/utility/array_algo.hpp b/ports-of-call/utility/array_algo.hpp
index 3e49a5c7..9cb2a516 100644
--- a/ports-of-call/utility/array_algo.hpp
+++ b/ports-of-call/utility/array_algo.hpp
@@ -20,6 +20,29 @@
 
 namespace util {
 
+// NB:
+// the array type is explicitly left unspecified, as it may be
+// desirable to use a type other than std::array.
+// however, it is required that the underlying type must
+// conform to a minimum std::array interface:
+//
+//  1.) constexpr operator[]
+//  2.) constexpr size_type size()
+//  3.) array_type::value_type
+//
+// A casual choice early on was that the array parameter is
+// assumed to be "complete" at instanciation. That is,
+// we use
+//
+//  template<class A>
+//
+// rather than (the more specific)
+//
+//  template<template<class, auto> A, class T, auto N>
+//
+// this imposes some structure that is a little loose, and may
+// be refactored at some point.
+
 // shorter make_index_sequence
 template <auto N>
 constexpr auto is(std::integral_constant<decltype(N), N>) {
@@ -42,7 +65,7 @@ PORTABLE_FORCEINLINE_FUNCTION constexpr auto array_map_impl(A const &x, F f,
   return std::array{f(x[Is])...};
 }
 
-template <class A, class B, std::size_t N, class F, std::size_t... Is>
+template <class A, class B, class F, std::size_t... Is>
 PORTABLE_FORCEINLINE_FUNCTION constexpr auto array_map_impl(A const &x, B const &y, F f,
                                                             std::index_sequence<Is...>) {
   return A{f(x[Is], y[Is])...};

From 703397a56674bd52b0746b42389f8d1e4f6ccca0 Mon Sep 17 00:00:00 2001
From: Christopher Mauney <mauneyc@pn2305185.lanl.gov>
Date: Wed, 4 Sep 2024 13:50:08 -0600
Subject: [PATCH 13/22] some sloppy juggling of types

---
 ports-of-call/portable_arrays.hpp    | 175 +++++++++++++++------------
 ports-of-call/utility/array_algo.hpp |  86 +++++++++++--
 ports-of-call/utility/index_algo.hpp |  13 +-
 test/test_mdidx.cpp                  |   6 +-
 test/test_portability.cpp            |   7 +-
 5 files changed, 185 insertions(+), 102 deletions(-)

diff --git a/ports-of-call/portable_arrays.hpp b/ports-of-call/portable_arrays.hpp
index 0938492a..ca14d3f2 100644
--- a/ports-of-call/portable_arrays.hpp
+++ b/ports-of-call/portable_arrays.hpp
@@ -27,6 +27,7 @@
 #include "array.hpp"
 #include "portability.hpp"
 #include "utility/array_algo.hpp"
+#include "utility/index_algo.hpp"
 #include <algorithm>
 #include <array>
 #include <assert.h>
@@ -47,22 +48,62 @@ constexpr std::size_t to_const = V;
 
 } // namespace detail
 
+#if __cplusplus == 202002L
+template <class T, auto N>
+using Array = std::array<T, N>;
+#else
 template <class T, auto N>
 using Array = PortsOfCall::array<T, N>;
+#endif // __cplusplus
 
 template <auto N>
 using IArray = Array<std::size_t, N>;
 
+/*
+#define ENABLE_CRTP(DerivedType) \
+    constexpr DerivedType &base() noexcept { return static_cast<DerivedType &>(*this); } \
+    constexpr const DerivedType &base() const noexcept { return static_cast<const
+DerivedType &>(*this); }
+
+
+
+template <class Derived>
+struct WithStoredStrides
+{
+  ENABLE_CRTP(Derived)
+  using size_type = typename Derived :: size_type;
+  using index_type = typename Derived :: index_type;
+  using strides_type = decltype(Derived::nx_);
+
+  strides_type strides_;
+
+  template<index_type ...Is>
+  constexpr auto operator() noexcept(Is...is)
+  {
+    return util::findex({is...}, base.nxs_, strides_);
+  }
+};
+*/
+
 template <typename T, std::size_t D = MAXDIM>
 class PortableMDArray {
  public:
+  using element_type = T;
+  using value_type = typename std::remove_cv<T>::type;
+  using index_type = std::ptrdiff_t;
+  using size_type = std::size_t;
+  using pointer = T *;
+  using const_pointer = const T *;
+  using reference = T &;
+  using const_reference = const T &;
+
   // explicit initialization of objects
   PORTABLE_FUNCTION PortableMDArray(T *data, IArray<D> extents, IArray<D> strides,
-                                    std::size_t rank) noexcept
+                                    size_type rank) noexcept
       : pdata_(data), nxs_(extents), strides_(strides), rank_(rank) {}
 
   // variadic ctor, dispatch to explicit constructor
-  template <typename... NXs, std::size_t N = sizeof...(NXs)>
+  template <typename... NXs, size_type N = sizeof...(NXs)>
   PORTABLE_FUNCTION PortableMDArray(T *p, NXs... nxs) noexcept {
     NewPortableMDArray(p, nxs...);
   } //: PortableMDArray(p, make_nxs_array(nxs...), make_strides_array<N>(), N) {
@@ -111,62 +152,81 @@ class PortableMDArray {
   // functions to get array dimensions
   template <std::size_t I>
   PORTABLE_FUNCTION constexpr auto GetDim() const {
-    return nxs_[I];
+    return nxs_[rank_ - I];
   }
 
   // legacy API, TODO: deprecate
-  PORTABLE_FORCEINLINE_FUNCTION int GetDim1() const { return GetDim<0>(); }
-  PORTABLE_FORCEINLINE_FUNCTION int GetDim2() const { return GetDim<1>(); }
-  PORTABLE_FORCEINLINE_FUNCTION int GetDim3() const { return GetDim<2>(); }
-  PORTABLE_FORCEINLINE_FUNCTION int GetDim4() const { return GetDim<3>(); }
-  PORTABLE_FORCEINLINE_FUNCTION int GetDim5() const { return GetDim<4>(); }
-  PORTABLE_FORCEINLINE_FUNCTION int GetDim6() const { return GetDim<5>(); }
+  [[deprecated("Use GetDim<N> instead.")]]
+  PORTABLE_FORCEINLINE_FUNCTION int GetDim1() const {
+    return GetDim<1>();
+  }
+  [[deprecated("Use GetDim<N> instead.")]]
+  PORTABLE_FORCEINLINE_FUNCTION int GetDim2() const {
+    return GetDim<2>();
+  }
+  [[deprecated("Use GetDim<N> instead.")]]
+  PORTABLE_FORCEINLINE_FUNCTION int GetDim3() const {
+    return GetDim<3>();
+  }
+  [[deprecated("Use GetDim<N> instead.")]]
+  PORTABLE_FORCEINLINE_FUNCTION int GetDim4() const {
+    return GetDim<4>();
+  }
+  [[deprecated("Use GetDim<N> instead.")]]
+  PORTABLE_FORCEINLINE_FUNCTION int GetDim5() const {
+    return GetDim<5>();
+  }
+  [[deprecated("Use GetDim<N> instead.")]]
+  PORTABLE_FORCEINLINE_FUNCTION int GetDim6() const {
+    return GetDim<6>();
+  }
   PORTABLE_INLINE_FUNCTION int GetDim(size_t i) const {
     // TODO: remove if performance cirtical
     assert(0 < i && i <= 6 && "PortableMDArrays are max 6D");
     switch (i) {
     case 1:
-      return GetDim1();
+      return GetDim<1>();
     case 2:
-      return GetDim2();
+      return GetDim<2>();
     case 3:
-      return GetDim3();
+      return GetDim<3>();
     case 4:
-      return GetDim4();
+      return GetDim<4>();
     case 5:
-      return GetDim5();
+      return GetDim<5>();
     case 6:
-      return GetDim6();
+      return GetDim<6>();
     }
     return -1;
   }
 
-  PORTABLE_FORCEINLINE_FUNCTION int GetSize() const {
-    return util::array_reduce(nxs_, 1, std::multiplies<std::size_t>{});
+  PORTABLE_FORCEINLINE_FUNCTION int GetSize() const noexcept {
+    return util::array_reduce(nxs_, 1, std::multiplies<size_type>{});
   }
-  PORTABLE_FORCEINLINE_FUNCTION std::size_t GetSizeInBytes() const {
+  PORTABLE_FORCEINLINE_FUNCTION size_type GetSizeInBytes() const noexcept {
     return GetSize() * sizeof(T);
   }
 
-  PORTABLE_INLINE_FUNCTION std::size_t GetRank() const { return rank_; }
+  PORTABLE_INLINE_FUNCTION size_type GetRank() const noexcept { return rank_; }
   template <typename... NXs>
   PORTABLE_INLINE_FUNCTION void Reshape(NXs... nxs) {
-    assert(util::array_reduce(IArray<D>{nxs...}, 1, std::multiplies<std::size_t>{}) ==
-           GetSize());
+    assert(util::array_reduce(IArray<D>{nxs...}, 1, std::multiplies<size_type>{}) ==
+               GetSize() &&
+           "New shape gives different size than existing size");
     update_layout(nxs...);
   }
 
-  PORTABLE_FORCEINLINE_FUNCTION bool IsShallowSlice() { return true; }
-  PORTABLE_FORCEINLINE_FUNCTION bool IsEmpty() { return GetSize() < 1; }
+  PORTABLE_FORCEINLINE_FUNCTION bool IsShallowSlice() const noexcept { return true; }
+  PORTABLE_FORCEINLINE_FUNCTION bool IsEmpty() const noexcept { return GetSize() < 1; }
   // "getter" function to access private data member
   // TODO(felker): Replace this unrestricted "getter" with a limited, safer
   // alternative.
   // TODO(felker): Rename function. Conflicts with "PortableMDArray<> data"
   // OutputData member.
-  PORTABLE_FORCEINLINE_FUNCTION T *data() { return pdata_; }
-  PORTABLE_FORCEINLINE_FUNCTION const T *data() const { return pdata_; }
-  PORTABLE_FORCEINLINE_FUNCTION T *begin() { return pdata_; }
-  PORTABLE_FORCEINLINE_FUNCTION T *end() { return pdata_ + GetSize(); }
+  PORTABLE_FORCEINLINE_FUNCTION pointer data() noexcept { return pdata_; }
+  PORTABLE_FORCEINLINE_FUNCTION const_pointer data() const noexcept { return pdata_; }
+  PORTABLE_FORCEINLINE_FUNCTION pointer begin() noexcept { return pdata_; }
+  PORTABLE_FORCEINLINE_FUNCTION pointer end() noexcept { return pdata_ + GetSize(); }
 
   // overload "function call" operator() to access 1d-5d data
   // provides Fortran-like syntax for multidimensional arrays vs. "subscript"
@@ -176,13 +236,13 @@ class PortableMDArray {
   // access via returning by reference, enabling assignment on returned
   // l-value, e.g.: a(3) = 3.0;
   template <typename... Is>
-  PORTABLE_FORCEINLINE_FUNCTION T &operator()(const Is... idxs) {
-    return pdata_[compute_index(idxs...)];
+  PORTABLE_FORCEINLINE_FUNCTION reference operator()(Is... idxs) noexcept {
+    return pdata_[util::fast_findex({static_cast<size_type>(idxs)...}, nxs_, strides_)];
   }
 
   template <typename... Is>
-  PORTABLE_FORCEINLINE_FUNCTION T &operator()(const Is... idxs) const {
-    return pdata_[compute_index(idxs...)];
+  PORTABLE_FORCEINLINE_FUNCTION T &operator()(Is... idxs) const noexcept {
+    return pdata_[util::fast_findex({static_cast<size_type>(idxs)...}, nxs_, strides_)];
   }
 
   PortableMDArray<T, D> &operator*=(T scale) {
@@ -240,60 +300,15 @@ class PortableMDArray {
   }
 
  private:
-  template <typename... NX, std::size_t N = sizeof...(NX)>
-  PORTABLE_FORCEINLINE_FUNCTION auto make_nxs_array(NX... nxs) {
-    IArray<D> a;
-    IArray<N> t{static_cast<std::size_t>(nxs)...};
-    for (auto i = 0; i < N; ++i) {
-      a[i] = t[N - i - 1];
-    }
-    for (auto i = N; i < D; ++i) {
-      a[i] = 1;
-    }
-    return a;
-  }
-
-  template <std::size_t N>
-  PORTABLE_INLINE_FUNCTION auto make_strides_array() {
-    IArray<D> a;
-    a[0] = 1;
-    for (auto i = 1; i < N; ++i) {
-      a[i] = a[i - 1] * nxs_[i - 1];
-    }
-    for (auto i = N; i < D; ++i)
-      a[i] = 0;
-
-    return a;
-  }
-
   // driver of nx array creation.
   // Note we iterate from the RIGHT, to match
   // what was done explicilt prior.
-  template <typename... NXs, std::size_t N = sizeof...(NXs)>
-  PORTABLE_INLINE_FUNCTION auto update_layout(NXs... nxs) {
+  template <typename... NXs, auto N = sizeof...(NXs)>
+  PORTABLE_INLINE_FUNCTION constexpr auto update_layout(NXs... nxs) noexcept {
 
     rank_ = N;
-    nxs_ = make_nxs_array(nxs...);
-    strides_ = make_strides_array<N>();
-  }
-
-  // compute_index base case, i.e. fastest moving index
-  template <std::size_t Ind, std::size_t N>
-  PORTABLE_FORCEINLINE_FUNCTION size_t compute_index_impl(const size_t index) const {
-    return index;
-  }
-
-  // compute_index general case, computing slower moving index strides
-  template <std::size_t Ind, std::size_t N, typename... Tail>
-  PORTABLE_FORCEINLINE_FUNCTION size_t compute_index_impl(const size_t index,
-                                                          const Tail... tail) const {
-    return index * strides_[N - Ind - 1] + compute_index_impl<Ind + 1, N>(tail...);
-  }
-  // compute index driver.
-  template <typename... Indicies, std::size_t N = sizeof...(Indicies)>
-  PORTABLE_FORCEINLINE_FUNCTION std::size_t compute_index(const Indicies... idxs) const {
-    // adding `0` if sizeof...(Indicies) == 0
-    return 0 + compute_index_impl<0, N>(idxs...);
+    nxs_ = util::make_underfilled_array<D>(util::wrap_vars<IArray<N>>(nxs...));
+    strides_ = util::make_underfilled_array<D, 0>(util::get_strides(nxs_));
   }
 
   T *pdata_;
diff --git a/ports-of-call/utility/array_algo.hpp b/ports-of-call/utility/array_algo.hpp
index 9cb2a516..e1d2e287 100644
--- a/ports-of-call/utility/array_algo.hpp
+++ b/ports-of-call/utility/array_algo.hpp
@@ -17,6 +17,7 @@
 #include "../portability.hpp"
 #include <array>
 #include <type_traits>
+#include <utility>
 
 namespace util {
 
@@ -53,6 +54,16 @@ constexpr auto is(std::integral_constant<decltype(N), N>) {
 template <class A>
 using value_t = typename A::value_type;
 
+template <class A>
+constexpr auto get_size(const A &a) {
+  return a.size();
+}
+
+template <class A, class... B>
+constexpr auto wrap_vars(B... vv) {
+  return A{static_cast<value_t<A>>(vv)...};
+}
+
 // determines the return type of Op(a,b)
 template <class A, class Op>
 using reduction_value_t =
@@ -83,8 +94,67 @@ PORTABLE_INLINE_FUNCTION constexpr auto array_reduce_impl(A const &x, Op op) {
   }
 }
 
+/*
+template <auto, auto V>
+constexpr auto NREP = V;
+
+template <class A, auto Fill, std::size_t... Is, std::size_t... Rs, class... B>
+PORTABLE_FORCEINLINE_FUNCTION constexpr auto
+make_underfilled_reverse_array_impl(std::index_sequence<Is...>,
+                                    std::index_sequence<Rs...>, B &&...vv) {
+    return wrap_vars<A>((std::get<sizeof...(Is) - 1 - Is>(std::tie(vv...)))...,
+                      NREP<Rs, Fill>...);
+}
+template <class A, auto Fill, std::size_t... Is, std::size_t... Rs, class... B>
+PORTABLE_FORCEINLINE_FUNCTION constexpr auto
+make_underfilled_array_impl(std::index_sequence<Is...>, std::index_sequence<Rs...>,
+                            B &&...vv) {
+    return wrap_vars<A>((std::get<Is>(std::tie(vv...)))..., NREP<Rs, Fill>...);
+}
+*/
+
 } // namespace detail
 
+/*
+template <class A, auto Fill = 1, class... B>
+PORTABLE_FORCEINLINE_FUNCTION constexpr auto make_underfilled_reverse_array(B... vv) {
+  constexpr auto D = get_size(A{});
+  return detail::make_underfilled_reverse_array_impl<A, Fill>(
+      std::index_sequence_for<B...>{}, std::make_index_sequence<D - sizeof...(B)>{},
+      vv...);
+}
+*/
+/*
+template <class A, auto Fill = 1, class... B>
+PORTABLE_FORCEINLINE_FUNCTION constexpr auto make_underfilled_array(B... vv) {
+  constexpr auto D = get_size(A{});
+  return detail::make_underfilled_array_impl<A, Fill>(
+      std::index_sequence_for<B...>{}, std::make_index_sequence<D - sizeof...(B)>{},
+      vv...);
+}
+*/
+template <auto P, auto Fill = 1, template <class, auto> class A, class T, auto O>
+PORTABLE_FORCEINLINE_FUNCTION constexpr decltype(auto)
+make_underfilled_array(const A<T, O> &in) {
+  A<T, P> out;
+  for (auto i = 0; i < in.size(); ++i)
+    out[i] = in[i];
+  for (auto i = in.size(); i < out.size(); ++i)
+    out[i] = Fill;
+  return out;
+}
+
+template <auto P, auto Fill = 1, template <class, auto> class A, class T, auto O>
+PORTABLE_FORCEINLINE_FUNCTION constexpr decltype(auto)
+make_underfilled_reversed_array(const A<T, O> &in) {
+  A<T, P> out;
+  for (auto i = in.size() - 1; i >= 0; i--)
+    out[i] = in[i];
+  for (auto i = in.size(); i < out.size(); ++i)
+    out[i] = Fill;
+  return out;
+}
+
 // maps an unary function f(x) to each array value, returning an array of results
 // x = {f(a[0]), f(a[1]),..., f(a[N-1])}
 template <class A, class F>
@@ -97,25 +167,25 @@ PORTABLE_FORCEINLINE_FUNCTION constexpr auto array_map(A const &x, F f) {
 
 template <class A, class B, class F>
 PORTABLE_FORCEINLINE_FUNCTION constexpr auto array_map(A const &x, B const &y, F f) {
-  return detail::array_map_impl(x, y, f, is(x.size()));
+  return detail::array_map_impl(x, y, f, std::make_index_sequence<get_size(A{})>{});
 }
 
-template <std::size_t I, class A, class Op, class T = reduction_value_t<A, Op>>
+template <std::size_t I, std::size_t N, class A, class Op,
+          class T = reduction_value_t<A, Op>>
 PORTABLE_FORCEINLINE_FUNCTION constexpr T array_partial_reduce(A x, T initial_value,
                                                                Op op) {
-  static_assert(I <= x.size());
-  if constexpr (I == 0)
+  static_assert(N <= x.size());
+  if constexpr ((N - I) == 0)
     return initial_value;
   else
-    return detail::array_reduce_impl<0, I>(x, op);
+    return detail::array_reduce_impl<I, N>(x, op);
 }
 
-template <class A, class Op, class T = reduction_value_t<A, Op>>
-
 // performs a reduction on an array of values
 // e.g. x = sum_i a[i]
+template <class A, class Op, class T = reduction_value_t<A, Op>>
 PORTABLE_FORCEINLINE_FUNCTION constexpr T array_reduce(A x, T initial_value, Op op) {
-  return array_partial_reduce<x.size()>(x, initial_value, op);
+  return array_partial_reduce<0, x.size()>(x, initial_value, op);
 }
 
 } // namespace util
diff --git a/ports-of-call/utility/index_algo.hpp b/ports-of-call/utility/index_algo.hpp
index 168676b9..ec5fc0bc 100644
--- a/ports-of-call/utility/index_algo.hpp
+++ b/ports-of-call/utility/index_algo.hpp
@@ -16,9 +16,6 @@
 
 #include "../portability.hpp"
 #include "array_algo.hpp"
-#include <array>
-#include <numeric>
-#include <type_traits>
 
 namespace util {
 
@@ -27,7 +24,8 @@ template <auto I, class A>
 PORTABLE_FORCEINLINE_FUNCTION static constexpr auto get_stride(A const &dim) {
 
   // NB: column major
-  return array_partial_reduce<I>(dim, value_t<A>{1}, std::multiplies<std::size_t>{});
+  return array_partial_reduce<I + 1, get_size(A{})>(dim, value_t<A>{1},
+                                                    std::multiplies<std::size_t>{});
 }
 
 namespace detail {
@@ -36,11 +34,12 @@ PORTABLE_FORCEINLINE_FUNCTION static constexpr auto
 get_strides_impl(A const &dim, std::index_sequence<Is...>) {
   return A{get_stride<Is>(dim)...};
 }
+} // namespace detail
+
 template <class A>
 PORTABLE_FORCEINLINE_FUNCTION static constexpr auto get_strides(A const &dim) {
-  return detail::get_strides_impl(dim, is(dim.size()));
+  return detail::get_strides_impl(dim, std::make_index_sequence<get_size(A{})>{});
 }
-} // namespace detail
 
 // returns the flat (1D) index of the md index set {i,j,k}
 // NB: fast because the strides are provided and don't need to be recomputed
@@ -49,7 +48,7 @@ PORTABLE_FORCEINLINE_FUNCTION static constexpr auto
 fast_findex(A const &ijk, A const &dim, A const &stride) {
   // TODO: assert ijk in bounds
   return array_reduce(array_map(ijk, stride, [](auto a, auto b) { return a * b; }),
-                      value_t<A>{1}, std::plus<std::size_t>{});
+                      value_t<A>{0}, std::plus<std::size_t>{});
 }
 
 // same as fast_findex, except the strides are calculated on the fly
diff --git a/test/test_mdidx.cpp b/test/test_mdidx.cpp
index dbf9b964..c7d9b341 100644
--- a/test/test_mdidx.cpp
+++ b/test/test_mdidx.cpp
@@ -49,9 +49,9 @@ TEST_CASE("PortableMDArrays Sizes are sane", "[PortableMDArray]") {
 
     REQUIRE(pmd3.GetSize() == subsz[2]);
     REQUIRE(pmd3.GetRank() == 3);
-    REQUIRE(pmd3.GetDim1() == nxs[2]);
-    REQUIRE(pmd3.GetDim2() == nxs[1]);
-    REQUIRE(pmd3.GetDim3() == nxs[0]);
+    REQUIRE(pmd3.GetDim<1>() == nxs[2]);
+    REQUIRE(pmd3.GetDim<2>() == nxs[1]);
+    REQUIRE(pmd3.GetDim<3>() == nxs[0]);
   }
 }
 
diff --git a/test/test_portability.cpp b/test/test_portability.cpp
index ac7bd60d..1135426f 100644
--- a/test/test_portability.cpp
+++ b/test/test_portability.cpp
@@ -51,8 +51,8 @@ TEST_CASE("PortableMDArrays can be allocated from a pointer", "[PortableMDArray]
   a.NewPortableMDArray(data.data(), M, N);
 
   SECTION("Shape should be NxM") {
-    REQUIRE(a.GetDim1() == N);
-    REQUIRE(a.GetDim2() == M);
+    REQUIRE(a.GetDim<1>() == N);
+    REQUIRE(a.GetDim<2>() == M);
   }
 
   SECTION("Stride is as set by initialized pointer") {
@@ -87,8 +87,7 @@ TEST_CASE("portableCopy works with all portability strategies", "[portableCopy]"
   Real *a = (Real *)PORTABLE_MALLOC(Nb);
 
   // set device values to 0
-  portableFor(
-      "set to 0", 0, N, PORTABLE_LAMBDA(const int &i) { a[i] = 0.0; });
+  portableFor("set to 0", 0, N, PORTABLE_LAMBDA(const int &i) { a[i] = 0.0; });
 
   // set host values to reference
   for (size_t i = 0; i < N; ++i) {

From 195fecf157f5caeaddd362b2aae4bde17420ed1d Mon Sep 17 00:00:00 2001
From: Christopher Mauney <mauneyc@pn2305185.lanl.gov>
Date: Tue, 8 Oct 2024 13:44:26 -0600
Subject: [PATCH 14/22] working on stuff

---
 ports-of-call/portable_arrays.hpp    |  73 ++++-------
 ports-of-call/utility/array_algo.hpp | 101 ++++++++++------
 ports-of-call/utility/index_algo.hpp |  27 ++++-
 test/test_mdidx.cpp                  | 174 +++++++++++++++++++++++----
 test/test_portability.cpp            |   4 +-
 test/test_utilities.hpp              |  27 ++---
 6 files changed, 272 insertions(+), 134 deletions(-)

diff --git a/ports-of-call/portable_arrays.hpp b/ports-of-call/portable_arrays.hpp
index ca14d3f2..6abab7f4 100644
--- a/ports-of-call/portable_arrays.hpp
+++ b/ports-of-call/portable_arrays.hpp
@@ -29,15 +29,14 @@
 #include "utility/array_algo.hpp"
 #include "utility/index_algo.hpp"
 #include <algorithm>
-#include <array>
 #include <assert.h>
 #include <cstddef> // size_t
 #include <cstring> // memset()
 #include <functional>
-#include <numeric>
 #include <type_traits>
 #include <utility> // swap()
 
+namespace PortsOfCall {
 // maximum number of dimensions
 constexpr std::size_t MAXDIM = 6;
 
@@ -48,6 +47,8 @@ constexpr std::size_t to_const = V;
 
 } // namespace detail
 
+// if C++20, we can use std::array without worry on device
+// else use a hand-rolled array
 #if __cplusplus == 202002L
 template <class T, auto N>
 using Array = std::array<T, N>;
@@ -59,35 +60,10 @@ using Array = PortsOfCall::array<T, N>;
 template <auto N>
 using IArray = Array<std::size_t, N>;
 
-/*
-#define ENABLE_CRTP(DerivedType) \
-    constexpr DerivedType &base() noexcept { return static_cast<DerivedType &>(*this); } \
-    constexpr const DerivedType &base() const noexcept { return static_cast<const
-DerivedType &>(*this); }
-
-
-
-template <class Derived>
-struct WithStoredStrides
-{
-  ENABLE_CRTP(Derived)
-  using size_type = typename Derived :: size_type;
-  using index_type = typename Derived :: index_type;
-  using strides_type = decltype(Derived::nx_);
-
-  strides_type strides_;
-
-  template<index_type ...Is>
-  constexpr auto operator() noexcept(Is...is)
-  {
-    return util::findex({is...}, base.nxs_, strides_);
-  }
-};
-*/
-
 template <typename T, std::size_t D = MAXDIM>
 class PortableMDArray {
  public:
+  using this_type = PortableMDArray<T, D>;
   using element_type = T;
   using value_type = typename std::remove_cv<T>::type;
   using index_type = std::ptrdiff_t;
@@ -103,7 +79,7 @@ class PortableMDArray {
       : pdata_(data), nxs_(extents), strides_(strides), rank_(rank) {}
 
   // variadic ctor, dispatch to explicit constructor
-  template <typename... NXs, size_type N = sizeof...(NXs)>
+  template <class... NXs, size_type N = sizeof...(NXs)>
   PORTABLE_FUNCTION PortableMDArray(T *p, NXs... nxs) noexcept {
     NewPortableMDArray(p, nxs...);
   } //: PortableMDArray(p, make_nxs_array(nxs...), make_strides_array<N>(), N) {
@@ -117,16 +93,13 @@ class PortableMDArray {
   // copies.
   //  PortableMDArray(const PortableMDArray<T> &t) noexcept;
   //  PortableMDArray<T> &operator=(const PortableMDArray<T> &t) noexcept;
-  PortableMDArray(const PortableMDArray<T, D> &src) noexcept {
-    nxs_ = src.nxs_;
-    rank_ = src.rank_;
-    strides_ = src.strides_;
-    if (src.pdata_) pdata_ = src.pdata_;
-  }
+  PortableMDArray(const this_type &src) noexcept
+      : nxs_(src.nxs_), rank_(src.rank_), strides_(src.strides_),
+        pdata_(src.pdata_ ? src.pdata_ : nullptr) {}
 
-  PortableMDArray<T, D> &operator=(const PortableMDArray<T, D> &src) noexcept {
+  PortableMDArray<T, D> &operator=(const this_type &src) noexcept {
     if (this != &src) {
-      nxs_ = src.nxs;
+      nxs_ = src.nxs_;
       rank_ = src.rank_;
       pdata_ = src.pdata_;
     }
@@ -144,7 +117,7 @@ class PortableMDArray {
   //  (shallow swap)
   // Does not allocate memory for either array
 
-  void SwapPortableMDArray(PortableMDArray<T, D> &array2) {
+  void SwapPortableMDArray(this_type &array2) {
     std::swap(pdata_, array2.pdata_);
     return;
   }
@@ -235,29 +208,29 @@ class PortableMDArray {
   // "non-const variants" called for "PortableMDArray<T>()" provide read/write
   // access via returning by reference, enabling assignment on returned
   // l-value, e.g.: a(3) = 3.0;
-  template <typename... Is>
-  PORTABLE_FORCEINLINE_FUNCTION reference operator()(Is... idxs) noexcept {
+  template <typename... Is> //, class = std::enable_if_t<(sizeof...(Is) > 3)>>
+  PORTABLE_FORCEINLINE_FUNCTION constexpr reference operator()(Is... idxs) noexcept {
     return pdata_[util::fast_findex({static_cast<size_type>(idxs)...}, nxs_, strides_)];
   }
 
-  template <typename... Is>
-  PORTABLE_FORCEINLINE_FUNCTION T &operator()(Is... idxs) const noexcept {
+  template <typename... Is> //, class = std::enable_if_t<(sizeof...(Is) > 3)>>
+  PORTABLE_FORCEINLINE_FUNCTION constexpr T &operator()(Is... idxs) const noexcept {
     return pdata_[util::fast_findex({static_cast<size_type>(idxs)...}, nxs_, strides_)];
   }
 
-  PortableMDArray<T, D> &operator*=(T scale) {
+  this_type &operator*=(T scale) {
     std::transform(pdata_, pdata_ + GetSize(), pdata_,
                    [scale](T val) { return scale * val; });
     return *this;
   }
 
-  PortableMDArray<T, D> &operator+=(const PortableMDArray<T, D> &other) {
+  this_type &operator+=(const this_type &other) {
     assert(GetSize() == other.GetSize());
     std::transform(pdata_, pdata_ + GetSize(), other.pdata_, pdata_, std::plus<T>());
     return *this;
   }
 
-  PortableMDArray<T, D> &operator-=(const PortableMDArray<T, D> &other) {
+  this_type &operator-=(const this_type &other) {
     assert(GetSize() == other.GetSize());
     std::transform(pdata_, pdata_ + GetSize(), other.pdata_, pdata_, std::minus<T>());
     return *this;
@@ -265,11 +238,11 @@ class PortableMDArray {
 
   // Checks that arrays point to same data with same shape
   // note this POINTER equivalence, not data equivalence
-  bool operator==(const PortableMDArray<T, D> &rhs) const {
+  bool operator==(const this_type &rhs) const {
     return (pdata_ == rhs.pdata_ && nxs_ == rhs.nxs_); // NB rank is implied
   }
 
-  bool operator!=(const PortableMDArray<T, D> &other) const { return !(*this == other); }
+  bool operator!=(const this_type &other) const { return !(*this == other); }
 
   //----------------------------------------------------------------------------------------
   //! \fn PortableMDArray::InitWithShallowSlice()
@@ -280,9 +253,8 @@ class PortableMDArray {
   //  entries of the src array for d<dim (cannot access any nx4=2, etc. entries
   //  if dim=3 for example)
 
-  PORTABLE_FUNCTION void InitWithShallowSlice(const PortableMDArray<T, D> &src,
-                                              const int dim, const int indx,
-                                              const int nvar) {
+  PORTABLE_FUNCTION void InitWithShallowSlice(const this_type &src, const int dim,
+                                              const int indx, const int nvar) {
     pdata_ = src.pdata_;
     std::size_t offs = indx;
     nxs_[dim - 1] = nvar;
@@ -319,4 +291,5 @@ class PortableMDArray {
  public:
 };
 
+} // namespace PortsOfCall
 #endif // _PORTABLE_ARRAYS_HPP_
diff --git a/ports-of-call/utility/array_algo.hpp b/ports-of-call/utility/array_algo.hpp
index e1d2e287..800f94de 100644
--- a/ports-of-call/utility/array_algo.hpp
+++ b/ports-of-call/utility/array_algo.hpp
@@ -54,13 +54,16 @@ constexpr auto is(std::integral_constant<decltype(N), N>) {
 template <class A>
 using value_t = typename A::value_type;
 
+// constexpr helper to get compile-time size of array A
 template <class A>
-constexpr auto get_size(const A &a) {
+PORTABLE_FORCEINLINE_FUNCTION constexpr auto get_size(const A &a) {
   return a.size();
 }
 
+// wraps a variadic set of variables into an array,
+// making an attempt to cast them.
 template <class A, class... B>
-constexpr auto wrap_vars(B... vv) {
+PORTABLE_FORCEINLINE_FUNCTION constexpr auto wrap_vars(B... vv) {
   return A{static_cast<value_t<A>>(vv)...};
 }
 
@@ -70,31 +73,39 @@ using reduction_value_t =
     decltype(std::declval<Op>()(std::declval<value_t<A>>(), std::declval<value_t<A>>()));
 
 namespace detail {
-template <class A, class F, std::size_t... Is>
+template <class A, class F, auto... Is>
 PORTABLE_FORCEINLINE_FUNCTION constexpr auto array_map_impl(A const &x, F f,
                                                             std::index_sequence<Is...>) {
   return std::array{f(x[Is])...};
 }
 
-template <class A, class B, class F, std::size_t... Is>
+template <class A, class B, class F, auto... Is>
 PORTABLE_FORCEINLINE_FUNCTION constexpr auto array_map_impl(A const &x, B const &y, F f,
                                                             std::index_sequence<Is...>) {
   return A{f(x[Is], y[Is])...};
 }
 
-template <std::size_t f, std::size_t l, class A, class Op>
+template <auto f, auto l, class A, class Op>
 PORTABLE_INLINE_FUNCTION constexpr auto array_reduce_impl(A const &x, Op op) {
   if constexpr ((l - f) == 1)
     return x[f];
   else {
-    constexpr std::size_t n = l - f;
+    constexpr auto n = l - f;
     auto left_sum = array_reduce_impl<f, f + n / 2>(x, op);
     auto right_sum = array_reduce_impl<f + n / 2, l>(x, op);
     return op(left_sum, right_sum);
   }
 }
 
+} // namespace detail
+
+///////////////////////////////////////////
+/// CMM: This verbose code is probably too
+/// messy to keep, but I liked it so I wanted
+/// to keep a record on a commit somehwere if
+/// i ever wanted to lift it.
 /*
+namespace detail{
 template <auto, auto V>
 constexpr auto NREP = V;
 
@@ -111,11 +122,8 @@ make_underfilled_array_impl(std::index_sequence<Is...>, std::index_sequence<Rs..
                             B &&...vv) {
     return wrap_vars<A>((std::get<Is>(std::tie(vv...)))..., NREP<Rs, Fill>...);
 }
-*/
+} //
 
-} // namespace detail
-
-/*
 template <class A, auto Fill = 1, class... B>
 PORTABLE_FORCEINLINE_FUNCTION constexpr auto make_underfilled_reverse_array(B... vv) {
   constexpr auto D = get_size(A{});
@@ -123,8 +131,7 @@ PORTABLE_FORCEINLINE_FUNCTION constexpr auto make_underfilled_reverse_array(B...
       std::index_sequence_for<B...>{}, std::make_index_sequence<D - sizeof...(B)>{},
       vv...);
 }
-*/
-/*
+
 template <class A, auto Fill = 1, class... B>
 PORTABLE_FORCEINLINE_FUNCTION constexpr auto make_underfilled_array(B... vv) {
   constexpr auto D = get_size(A{});
@@ -133,27 +140,7 @@ PORTABLE_FORCEINLINE_FUNCTION constexpr auto make_underfilled_array(B... vv) {
       vv...);
 }
 */
-template <auto P, auto Fill = 1, template <class, auto> class A, class T, auto O>
-PORTABLE_FORCEINLINE_FUNCTION constexpr decltype(auto)
-make_underfilled_array(const A<T, O> &in) {
-  A<T, P> out;
-  for (auto i = 0; i < in.size(); ++i)
-    out[i] = in[i];
-  for (auto i = in.size(); i < out.size(); ++i)
-    out[i] = Fill;
-  return out;
-}
-
-template <auto P, auto Fill = 1, template <class, auto> class A, class T, auto O>
-PORTABLE_FORCEINLINE_FUNCTION constexpr decltype(auto)
-make_underfilled_reversed_array(const A<T, O> &in) {
-  A<T, P> out;
-  for (auto i = in.size() - 1; i >= 0; i--)
-    out[i] = in[i];
-  for (auto i = in.size(); i < out.size(); ++i)
-    out[i] = Fill;
-  return out;
-}
+///////////////////////////////////////////
 
 // maps an unary function f(x) to each array value, returning an array of results
 // x = {f(a[0]), f(a[1]),..., f(a[N-1])}
@@ -164,28 +151,62 @@ PORTABLE_FORCEINLINE_FUNCTION constexpr auto array_map(A const &x, F f) {
 
 // maps a binary function to each array value, returning an array of results
 // x = {f(a[0], b[0]), f(a[1], b[1]),..., f(a[N-1], b[N-1])}
-
 template <class A, class B, class F>
 PORTABLE_FORCEINLINE_FUNCTION constexpr auto array_map(A const &x, B const &y, F f) {
   return detail::array_map_impl(x, y, f, std::make_index_sequence<get_size(A{})>{});
 }
 
-template <std::size_t I, std::size_t N, class A, class Op,
-          class T = reduction_value_t<A, Op>>
+// does a reduction on the sub-array A[I..F)
+// also used to construct full reduction A[..N)
+template <auto I, auto F, class A, class Op, class T = reduction_value_t<A, Op>>
 PORTABLE_FORCEINLINE_FUNCTION constexpr T array_partial_reduce(A x, T initial_value,
                                                                Op op) {
-  static_assert(N <= x.size());
-  if constexpr ((N - I) == 0)
+  static_assert(F <= get_size(A{}));
+  if constexpr ((F - I) == 0)
     return initial_value;
   else
-    return detail::array_reduce_impl<I, N>(x, op);
+    return detail::array_reduce_impl<I, F>(x, op);
 }
 
 // performs a reduction on an array of values
 // e.g. x = sum_i a[i]
+// NB Op is evaluated as a magma, i.e. Op(a,b,c,d) = Op(a,b) + Op(c,d)
+// it executes O(ln n) operations. It's possible this obfuscates
+// potential compiler optimizations for small (2-3) sized arrays
 template <class A, class Op, class T = reduction_value_t<A, Op>>
 PORTABLE_FORCEINLINE_FUNCTION constexpr T array_reduce(A x, T initial_value, Op op) {
-  return array_partial_reduce<0, x.size()>(x, initial_value, op);
+  return array_partial_reduce<0, get_size(A{})>(x, initial_value, op);
+}
+
+// if input values are less than output array length,
+// this function returns an array with "filled" values.
+// e.g.
+//  in: {x, y}
+//  out: {x, y, 1, 1, 1}
+template <auto P, auto Fill = 1, template <class, auto> class A, class T, auto O>
+PORTABLE_FORCEINLINE_FUNCTION constexpr decltype(auto)
+make_underfilled_array(const A<T, O> &in) {
+  A<T, P> out;
+  for (auto i = 0; i < in.size(); ++i)
+    out[i] = in[i];
+  for (auto i = in.size(); i < out.size(); ++i)
+    out[i] = Fill;
+  return out;
+}
+
+// same as above, but reverses input values
+// e.g.
+//  in: {x, y, z}
+//  out: {z, y, x, 1, 1}
+template <auto P, auto Fill = 1, template <class, auto> class A, class T, auto O>
+PORTABLE_FORCEINLINE_FUNCTION constexpr decltype(auto)
+make_underfilled_reversed_array(const A<T, O> &in) {
+  A<T, P> out;
+  for (auto i = in.size() - 1; i >= 0; i--)
+    out[i] = in[i];
+  for (auto i = in.size(); i < out.size(); ++i)
+    out[i] = Fill;
+  return out;
 }
 
 } // namespace util
diff --git a/ports-of-call/utility/index_algo.hpp b/ports-of-call/utility/index_algo.hpp
index ec5fc0bc..fe2796ac 100644
--- a/ports-of-call/utility/index_algo.hpp
+++ b/ports-of-call/utility/index_algo.hpp
@@ -16,18 +16,21 @@
 
 #include "../portability.hpp"
 #include "array_algo.hpp"
+#include <numeric>
 
 namespace util {
 
 // returns the stride (offset) needed to move one unit in dimension D
 template <auto I, class A>
 PORTABLE_FORCEINLINE_FUNCTION static constexpr auto get_stride(A const &dim) {
-
-  // NB: column major
+  // NB: storage order row-major
   return array_partial_reduce<I + 1, get_size(A{})>(dim, value_t<A>{1},
                                                     std::multiplies<std::size_t>{});
 }
 
+// I will rejoice when C++20 arrives and we can
+// template lambdas instead of having to write
+// seperate "detail" and "impls" for index_sequences
 namespace detail {
 template <class A, std::size_t... Is>
 PORTABLE_FORCEINLINE_FUNCTION static constexpr auto
@@ -36,6 +39,12 @@ get_strides_impl(A const &dim, std::index_sequence<Is...>) {
 }
 } // namespace detail
 
+// return an array of strides
+// dim: input array of dimenstions
+// returns array of strides
+// e.g.
+//  dim = {NX, NY, NZ}
+//  return = {1, NX, NX * NY}
 template <class A>
 PORTABLE_FORCEINLINE_FUNCTION static constexpr auto get_strides(A const &dim) {
   return detail::get_strides_impl(dim, std::make_index_sequence<get_size(A{})>{});
@@ -44,11 +53,19 @@ PORTABLE_FORCEINLINE_FUNCTION static constexpr auto get_strides(A const &dim) {
 // returns the flat (1D) index of the md index set {i,j,k}
 // NB: fast because the strides are provided and don't need to be recomputed
 template <class A>
-PORTABLE_FORCEINLINE_FUNCTION static constexpr auto
-fast_findex(A const &ijk, A const &dim, A const &stride) {
+PORTABLE_FUNCTION static constexpr auto fast_findex(A const &ijk, A const &dim,
+                                                    A const &stride) {
   // TODO: assert ijk in bounds
+  constexpr auto N = get_size(A{});
+  if constexpr (N == 1) {
+    return ijk[0];
+  } else if constexpr (N == 2) {
+    return ijk[1] + ijk[0] * stride[0];
+  }
+  // return std::inner_product(std::begin(ijk), std::end(ijk), std::begin(stride),
+  //                           value_t<A>{0});
   return array_reduce(array_map(ijk, stride, [](auto a, auto b) { return a * b; }),
-                      value_t<A>{0}, std::plus<std::size_t>{});
+                      value_t<A>{0}, std::plus<value_t<A>>{});
 }
 
 // same as fast_findex, except the strides are calculated on the fly
diff --git a/test/test_mdidx.cpp b/test/test_mdidx.cpp
index c7d9b341..a5e504b1 100644
--- a/test/test_mdidx.cpp
+++ b/test/test_mdidx.cpp
@@ -13,9 +13,10 @@
 
 #include <ports-of-call/portability.hpp>
 #include <ports-of-call/portable_arrays.hpp>
-#include <utility>
+#include <random>
 #include <vector>
 
+#include "catch2/benchmark/catch_chronometer.hpp"
 #include "test_utilities.hpp"
 
 #ifndef CATCH_CONFIG_FAST_COMPILE
@@ -26,6 +27,8 @@
 #include <catch2/matchers/catch_matchers_floating_point.hpp>
 #endif
 
+using namespace PortsOfCall;
+
 TEST_CASE("PortableMDArrays Sizes are sane", "[PortableMDArray]") {
   using tape_t = std::vector<int>;
 
@@ -61,9 +64,9 @@ TEST_CASE("Correct portable indexing", "[PortableMDArray]") {
     return [=](auto k, auto j, auto i) { return i + nx * j + nxny * k; };
   };
 
-  constexpr std::size_t NX = 32, NY = 64, NZ = 4;
+  constexpr std::size_t NX = 4, NY = 12, NZ = 3;
   constexpr std::size_t NC = NX * NY * NZ;
-  constexpr Real scale = 0.1;
+  constexpr Real scale = 1.0;
   // size in bytes
   constexpr const std::size_t NCb = NC * sizeof(Real);
 
@@ -74,9 +77,7 @@ TEST_CASE("Correct portable indexing", "[PortableMDArray]") {
     tape_ref[n] = scale * static_cast<Real>(n);
   }
 
-  // device pointer
   Real *tape_d = (Real *)PORTABLE_MALLOC(NCb);
-
   auto view_d = PortableMDArray<Real>(tape_d, NZ, NY, NX);
 
   // set device values
@@ -87,12 +88,14 @@ TEST_CASE("Correct portable indexing", "[PortableMDArray]") {
       });
 
   portableCopyToHost(tape_buf.data(), tape_d, NCb);
+  for (auto n = 0; n < NC; ++n) {
+    std::cout << "REF=" << tape_ref[n] << " BUF=" << tape_buf[n] << " n=" << n << "\n";
+  }
 
   for (auto n = 0; n < NC; ++n) {
     INFO("REF=" << tape_ref[n] << " BUF=" << tape_buf[n] << " n=" << n);
     REQUIRE_THAT(tape_buf[n], Catch::Matchers::WithinRel(tape_ref[n]));
   }
-
   PORTABLE_FREE(tape_d);
 }
 
@@ -102,43 +105,172 @@ TEST_CASE("Correct portable indexing", "[PortableMDArray]") {
 // (you may need to escape the `!` char depending
 // on your shell)
 TEST_CASE("Benchmark Indexing", "[!benchmark]") {
-  constexpr std::size_t NX = 32, NY = 64, NZ = 16;
+  constexpr std::size_t NX = 16, NY = 16, NZ = 78;
   constexpr std::size_t NU = 256, NV = 256, NW = 128;
 
   // construct a series of PortableMDArrays
   // of different rank and sizes, and
   // iterates through all contiguous indexing.
 
+  Real a = 1.0;
   // clang-format off
-  SECTION("small") {
-    BENCHMARK("index calc 1D") 
+  SECTION("small 1D") {
+    BENCHMARK_ADVANCED("baseline continguous")(Catch::Benchmark::Chronometer meter)
+    { 
+      Real *tape_d = (Real *)PORTABLE_MALLOC(sizeof(Real)*NX);
+
+      tape_d[0] = a;
+      meter.measure([=] {portableFor(
+        "fill 1D",
+        1, NX, 
+        PORTABLE_LAMBDA(auto i){tape_d[i]=a+tape_d[i-1];}); return 0;});
+
+      PORTABLE_FREE(tape_d);
+      return a;
+    };
+
+
+    BENCHMARK_ADVANCED("PMD continguous (space optimized)")(Catch::Benchmark::Chronometer meter)
     { 
-      return testing::idx_contiguous_bm(std::array{NX}); 
+      Real *tape_d = (Real *)PORTABLE_MALLOC(sizeof(Real)*NX);
+      auto view_d = PortableMDArray<Real, 1>(tape_d, NX);
+
+      view_d(0) = a;
+      meter.measure([=] {portableFor(
+        "fill 1D",
+        1, NX, 
+        PORTABLE_LAMBDA(auto i){view_d(i)=a+view_d(i-1);});return 0;});
+
+      PORTABLE_FREE(tape_d);
     };
-    BENCHMARK("index calc 2D") 
+
+    BENCHMARK_ADVANCED("PMD continguous (default size)")(Catch::Benchmark::Chronometer meter)
     { 
-      return testing::idx_contiguous_bm(std::array{NX, NY}); 
+      Real *tape_d = (Real *)PORTABLE_MALLOC(sizeof(Real)*NX);
+      auto view_d = PortableMDArray<Real>(tape_d, NX);
+
+      view_d(0) = a;
+      meter.measure([=] {portableFor(
+        "fill 1D",
+        1, NX, 
+        PORTABLE_LAMBDA(auto i){view_d(i)=a+view_d(i-1);});return 0;});
+
+      auto b = view_d(0);
+      PORTABLE_FREE(tape_d);
     };
-    BENCHMARK("index calc 3D") 
+
+    BENCHMARK_ADVANCED("baseline shuffle")(Catch::Benchmark::Chronometer meter)
     { 
-      return testing::idx_contiguous_bm(std::array{NX, NY, NZ}); 
+      Real *tape_d = (Real *)PORTABLE_MALLOC(sizeof(Real)*NX);
+      std::array<std::size_t, NX> offs;
+      std::iota(std::begin(offs), std::end(offs), 0);
+      std::shuffle(std::begin(offs), std::end(offs), std::mt19937{std::random_device{}()});
+
+      meter.measure([=] {portableFor(
+        "fill 1D",
+        0, NX, 
+        PORTABLE_LAMBDA(auto i){tape_d[offs[i]]=i;}); return 0;});
+
+      PORTABLE_FREE(tape_d);
     };
+
+    BENCHMARK_ADVANCED("PMD shuffle (space optimized)")(Catch::Benchmark::Chronometer meter)
+    { 
+      Real *tape_d = (Real *)PORTABLE_MALLOC(sizeof(Real)*NX);
+      auto view_d = PortableMDArray<Real, 1>(tape_d, NX);
+
+      std::array<std::size_t, NX> offs;
+      std::iota(std::begin(offs), std::end(offs), 0);
+      std::shuffle(std::begin(offs), std::end(offs), std::mt19937{std::random_device{}()});
+
+      meter.measure([=] {portableFor(
+        "fill 1D",
+        1, NX, 
+        PORTABLE_LAMBDA(auto i){view_d(offs[i])=i;});return 0;});
+
+      PORTABLE_FREE(tape_d);
+    };
+
+
   }
+  SECTION("small 2D") {
+    BENCHMARK_ADVANCED("baseline continguous")(Catch::Benchmark::Chronometer meter)
+    { 
+      Real *tape_d = (Real *)PORTABLE_MALLOC(sizeof(Real)*NX*NY);
+
+      portableFor("prep 2d", 0, NY, 0, 1, PORTABLE_LAMBDA(auto j, auto i) { tape_d[i+j*NX] = a;});
+      meter.measure([=] {portableFor(
+        "fill 1D",
+        0, NY, 1, NX,
+        PORTABLE_LAMBDA(auto j, auto i){tape_d[i + j * NX]=a + tape_d[(i-1) + j * NX];}); return 0;});
 
-  SECTION("big") {
-    BENCHMARK("index calc 1D") 
+      PORTABLE_FREE(tape_d);
+    };
+    BENCHMARK_ADVANCED("PMD contiguous (space optimized)")(Catch::Benchmark::Chronometer meter)
     { 
-      return testing::idx_contiguous_bm(std::array{NU}); 
+      Real *tape_d = (Real *)PORTABLE_MALLOC(sizeof(Real)*NX*NY);
+      auto view_d = PortableMDArray<Real, 2>(tape_d, NY, NX);
+
+      portableFor("prep 2d", 0, NY, 0, 1, PORTABLE_LAMBDA(auto j, auto i) { view_d(j, i) = a;});
+      view_d(0,0) = a;
+
+      meter.measure([=] {portableFor(
+        "fill 2D",
+        0, NY, 1, NX,
+        PORTABLE_LAMBDA(auto j, auto i){view_d(j,i)=a + view_d(j, i-1);});return 0;});
+
+      PORTABLE_FREE(tape_d);
     };
-    BENCHMARK("index calc 2D") 
+
+    BENCHMARK_ADVANCED("PMD continguous (default size)")(Catch::Benchmark::Chronometer meter)
     { 
-      return testing::idx_contiguous_bm(std::array{NU, NV}); 
+      Real *tape_d = (Real *)PORTABLE_MALLOC(sizeof(Real)*NX*NY);
+      auto view_d = PortableMDArray<Real>(tape_d, NY, NX);
+
+      portableFor("prep 2d", 0, NY, 0, 1, PORTABLE_LAMBDA(auto j, auto i) { view_d(j, i) = a;});
+      view_d(0,0) = a;
+
+      meter.measure([=] {portableFor(
+        "fill 2D",
+        0, NY, 1, NX,
+        PORTABLE_LAMBDA(auto j, auto i){view_d(j,i)=a + view_d(j, i-1);});return 0;});
+
+      PORTABLE_FREE(tape_d);
     };
-    BENCHMARK("index calc 3D") 
+  }
+SECTION("small 5D"){
+
+    BENCHMARK_ADVANCED("baseline continguous")(Catch::Benchmark::Chronometer meter)
     { 
-      return testing::idx_contiguous_bm(std::array{NU, NV, NW}); 
+      Real *tape_d = (Real *)PORTABLE_MALLOC(sizeof(Real)*NX*NX*NY*NY*NZ);
+
+      portableFor(
+        "prep 5D",
+        0, NZ, 0, NY, 0, NY, 0, NX, 0, 1,
+        PORTABLE_LAMBDA(auto m, auto n, auto k, auto j, auto i){tape_d[i + j * NX + k * NX * NX + n * NX * NX * NY + m * NY * NY * NX * NX ]=a;});
+
+      meter.measure([=] {portableFor(
+        "fill 5D",
+        0, NZ, 0, NY, 0, NY, 0, NX, 1, NX,
+        PORTABLE_LAMBDA(auto m, auto n, auto k, auto j, auto i){tape_d[i + j * NX + k * NX * NX + n * NX * NX * NY + m * NY * NY * NX * NX ]=a + tape_d[i + j * NX + k * NX * NX + n * NX * NX * NY + m * NY * NY * NX * NX ];});return 0;});
+
+      PORTABLE_FREE(tape_d);
     };
 
+    BENCHMARK_ADVANCED("PMD continguous")(Catch::Benchmark::Chronometer meter)
+    { 
+      Real *tape_d = (Real *)PORTABLE_MALLOC(sizeof(Real)*NX*NX*NY*NY*NZ);
+      auto view_d = PortableMDArray<Real, 5>(tape_d, NZ, NY, NY, NX, NX);
+
+      portableFor("prep 5d", 0, NZ, 0, NY, 0, NY, 0, NX, 0, 1, PORTABLE_LAMBDA(auto m, auto n, auto k, auto j, auto i) { view_d(m,n,k,j, i) = a;});
+      meter.measure([=] {portableFor(
+        "fill 5D",
+        0, NZ, 0, NY, 0, NY, 0, NX, 1, NX,
+        PORTABLE_LAMBDA(auto m, auto n, auto k, auto j, auto i){view_d(m,n,k,j,i)=a + view_d(m,n,k,j,i-1);});return 0;});
+
+      PORTABLE_FREE(tape_d);
+    };
   }
+
   // clang-format on
 }
diff --git a/test/test_portability.cpp b/test/test_portability.cpp
index 1135426f..19947768 100644
--- a/test/test_portability.cpp
+++ b/test/test_portability.cpp
@@ -42,7 +42,7 @@ TEST_CASE("PortableMDArrays can be allocated from a pointer", "[PortableMDArray]
   constexpr int N = 2;
   constexpr int M = 3;
   std::vector<int> data(N * M);
-  PortableMDArray<int> a;
+  PortsOfCall::PortableMDArray<int> a;
   int tot = 0;
   for (int i = 0; i < N * M; i++) {
     data[i] = tot;
@@ -66,7 +66,7 @@ TEST_CASE("PortableMDArrays can be allocated from a pointer", "[PortableMDArray]
   }
 
   SECTION("Identical slices of the same data should compare equal") {
-    PortableMDArray<int> aslc1, aslc2;
+    PortsOfCall::PortableMDArray<int> aslc1, aslc2;
     aslc1.InitWithShallowSlice(a, 1, 0, 2);
     aslc2.InitWithShallowSlice(a, 1, 0, 2);
     REQUIRE(aslc1 == aslc2);
diff --git a/test/test_utilities.hpp b/test/test_utilities.hpp
index 30ba2090..7d2023a1 100644
--- a/test/test_utilities.hpp
+++ b/test/test_utilities.hpp
@@ -14,10 +14,10 @@
 #ifndef _TEST_TEST_UTILITIES_HPP_
 #define _TEST_TEST_UTILITIES_HPP_
 
-#include <functional>
 #include <ports-of-call/array.hpp>
 #include <ports-of-call/portability.hpp>
 #include <ports-of-call/portable_arrays.hpp>
+#include <type_traits>
 #include <utility>
 
 namespace testing {
@@ -66,27 +66,22 @@ PORTABLE_INLINE_FUNCTION constexpr auto pf_invoke(View v, const Array &ext_par,
 template <class Ptr, class Array, std::size_t... I>
 PORTABLE_FORCEINLINE_FUNCTION constexpr decltype(auto) mdview(Ptr *d, const Array &arr,
                                                               std::index_sequence<I...>) {
-  return PortableMDArray<Real>(d, arr[I]...);
+  return PortsOfCall::PortableMDArray<Real>(d, arr[I]...);
 }
 
-// stand up and strip down a benchmark
-// input is a set of sizes
-template <class T, std::size_t N, typename I1 = std::make_index_sequence<N>,
-          typename I2 = std::make_index_sequence<2 * N>>
-PORTABLE_FUNCTION auto idx_contiguous_bm(const std::array<T, N> &nxa) {
-  auto nc = std::accumulate(std::begin(nxa), std::end(nxa), 1, std::multiplies<T>{});
+// stand up a benchmark
+template <class... Ns, std::size_t N = sizeof...(Ns),
+          typename I1 = std::make_index_sequence<N>,
+          typename I2 = std::make_index_sequence<2 * N>,
+          class = std::enable_if_t<std::conjunction_v<std::is_integral<Ns>...>>>
+PORTABLE_FUNCTION auto alloc_tape(Ns... ns) {
+  const auto nc = (ns * ... * std::size_t{1});
 
   Real *tape_d = (Real *)PORTABLE_MALLOC(nc * sizeof(Real));
 
-  auto view_d = mdview(tape_d, nxa, I1{});
+  auto view_d = mdview(tape_d, std::array{ns...}, I1{});
 
-  pf_invoke(view_d, nxa, I2{});
-
-  PORTABLE_FREE(tape_d);
-
-  // return here is arbitrary, Catch2 recommends this to avoid
-  // optimizing out a call.
-  return 1;
+  return std::tuple(tape_d, view_d);
 }
 
 } // namespace testing

From f2431e13894243638e49a69e17322ff03c8016f2 Mon Sep 17 00:00:00 2001
From: Christopher Mauney <mauneyc@pn2305185.lanl.gov>
Date: Fri, 15 Nov 2024 09:23:00 -0700
Subject: [PATCH 15/22] add span

---
 ports-of-call/portable_arrays.hpp | 5 ++---
 1 file changed, 2 insertions(+), 3 deletions(-)

diff --git a/ports-of-call/portable_arrays.hpp b/ports-of-call/portable_arrays.hpp
index 6abab7f4..c0c3ae0c 100644
--- a/ports-of-call/portable_arrays.hpp
+++ b/ports-of-call/portable_arrays.hpp
@@ -26,6 +26,7 @@
 
 #include "array.hpp"
 #include "portability.hpp"
+#include "span.hpp"
 #include "utility/array_algo.hpp"
 #include "utility/index_algo.hpp"
 #include <algorithm>
@@ -283,12 +284,10 @@ class PortableMDArray {
     strides_ = util::make_underfilled_array<D, 0>(util::get_strides(nxs_));
   }
 
-  T *pdata_;
+  span<T> pdata_;
   IArray<D> nxs_;
   IArray<D> strides_;
   int rank_;
-
- public:
 };
 
 } // namespace PortsOfCall

From 3c769f2da4760cf46110fbaae7a30cec7ff56cf8 Mon Sep 17 00:00:00 2001
From: Christopher Mauney <mauneyc@pn2305185.lanl.gov>
Date: Tue, 19 Nov 2024 09:00:03 -0700
Subject: [PATCH 16/22] undo span changes

---
 ports-of-call/portable_arrays.hpp | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/ports-of-call/portable_arrays.hpp b/ports-of-call/portable_arrays.hpp
index c0c3ae0c..ae2837ef 100644
--- a/ports-of-call/portable_arrays.hpp
+++ b/ports-of-call/portable_arrays.hpp
@@ -284,7 +284,7 @@ class PortableMDArray {
     strides_ = util::make_underfilled_array<D, 0>(util::get_strides(nxs_));
   }
 
-  span<T> pdata_;
+  T *pdata_;
   IArray<D> nxs_;
   IArray<D> strides_;
   int rank_;

From 9733f8044ce56308299097346d9bcf0d548f1b21 Mon Sep 17 00:00:00 2001
From: Christopher Michael Mauney <mauneyc@lanl.gov>
Date: Tue, 19 Nov 2024 10:56:19 -0700
Subject: [PATCH 17/22] removed benchmark testing; fix some tests for gpu

---
 ports-of-call/portable_arrays.hpp |  32 +++---
 test/CMakeLists.txt               |   1 +
 test/test_mdidx.cpp               | 181 +-----------------------------
 test/test_utilities.hpp           |   2 +-
 4 files changed, 21 insertions(+), 195 deletions(-)

diff --git a/ports-of-call/portable_arrays.hpp b/ports-of-call/portable_arrays.hpp
index ae2837ef..f3ba64c5 100644
--- a/ports-of-call/portable_arrays.hpp
+++ b/ports-of-call/portable_arrays.hpp
@@ -130,28 +130,28 @@ class PortableMDArray {
   }
 
   // legacy API, TODO: deprecate
-  [[deprecated("Use GetDim<N> instead.")]]
-  PORTABLE_FORCEINLINE_FUNCTION int GetDim1() const {
+  [[deprecated("Use GetDim<N> instead.")]] PORTABLE_FORCEINLINE_FUNCTION int
+  GetDim1() const {
     return GetDim<1>();
   }
-  [[deprecated("Use GetDim<N> instead.")]]
-  PORTABLE_FORCEINLINE_FUNCTION int GetDim2() const {
+  [[deprecated("Use GetDim<N> instead.")]] PORTABLE_FORCEINLINE_FUNCTION int
+  GetDim2() const {
     return GetDim<2>();
   }
-  [[deprecated("Use GetDim<N> instead.")]]
-  PORTABLE_FORCEINLINE_FUNCTION int GetDim3() const {
+  [[deprecated("Use GetDim<N> instead.")]] PORTABLE_FORCEINLINE_FUNCTION int
+  GetDim3() const {
     return GetDim<3>();
   }
-  [[deprecated("Use GetDim<N> instead.")]]
-  PORTABLE_FORCEINLINE_FUNCTION int GetDim4() const {
+  [[deprecated("Use GetDim<N> instead.")]] PORTABLE_FORCEINLINE_FUNCTION int
+  GetDim4() const {
     return GetDim<4>();
   }
-  [[deprecated("Use GetDim<N> instead.")]]
-  PORTABLE_FORCEINLINE_FUNCTION int GetDim5() const {
+  [[deprecated("Use GetDim<N> instead.")]] PORTABLE_FORCEINLINE_FUNCTION int
+  GetDim5() const {
     return GetDim<5>();
   }
-  [[deprecated("Use GetDim<N> instead.")]]
-  PORTABLE_FORCEINLINE_FUNCTION int GetDim6() const {
+  [[deprecated("Use GetDim<N> instead.")]] PORTABLE_FORCEINLINE_FUNCTION int
+  GetDim6() const {
     return GetDim<6>();
   }
   PORTABLE_INLINE_FUNCTION int GetDim(size_t i) const {
@@ -280,8 +280,12 @@ class PortableMDArray {
   PORTABLE_INLINE_FUNCTION constexpr auto update_layout(NXs... nxs) noexcept {
 
     rank_ = N;
-    nxs_ = util::make_underfilled_array<D>(util::wrap_vars<IArray<N>>(nxs...));
-    strides_ = util::make_underfilled_array<D, 0>(util::get_strides(nxs_));
+    // NB: the functions here have trouble with template argument deduction
+    // on some compilers, so they are hand-written.
+    auto szs = util::wrap_vars<IArray<N>>(nxs...);
+    nxs_ = util::make_underfilled_array<D, 1, Array, size_type, N>(szs);
+    auto sds = util::get_strides(nxs_);
+    strides_ = util::make_underfilled_array<D, 0, Array, size_type, D>(sds);
   }
 
   T *pdata_;
diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
index 22445319..39dfba58 100644
--- a/test/CMakeLists.txt
+++ b/test/CMakeLists.txt
@@ -68,5 +68,6 @@ target_sources(test_portsofcall
     test_span.cpp
     test_static_vector.cpp
     test_static_vector_iterator.cpp
+    test_mdidx.cpp
     "$<$<BOOL:${PORTABILITY_STRATEGY_KOKKOS}>:test_static_vector_kokkos.cpp>"
 )
diff --git a/test/test_mdidx.cpp b/test/test_mdidx.cpp
index a5e504b1..108f6a20 100644
--- a/test/test_mdidx.cpp
+++ b/test/test_mdidx.cpp
@@ -60,9 +60,6 @@ TEST_CASE("PortableMDArrays Sizes are sane", "[PortableMDArray]") {
 
 TEST_CASE("Correct portable indexing", "[PortableMDArray]") {
   // layout
-  auto iflat = [](auto nx, auto nxny) {
-    return [=](auto k, auto j, auto i) { return i + nx * j + nxny * k; };
-  };
 
   constexpr std::size_t NX = 4, NY = 12, NZ = 3;
   constexpr std::size_t NC = NX * NY * NZ;
@@ -84,7 +81,7 @@ TEST_CASE("Correct portable indexing", "[PortableMDArray]") {
   portableFor(
       "set unique val", 0, NZ, 0, NY, 0, NX,
       PORTABLE_LAMBDA(const int &k, const int &j, const int &i) {
-        view_d(k, j, i) = scale * iflat(NX, NX * NY)(k, j, i);
+        view_d(k, j, i) = scale * (i + NX * j + NX * NY * k);
       });
 
   portableCopyToHost(tape_buf.data(), tape_d, NCb);
@@ -98,179 +95,3 @@ TEST_CASE("Correct portable indexing", "[PortableMDArray]") {
   }
   PORTABLE_FREE(tape_d);
 }
-
-// runs benchmarks for indexing.
-// note: to run, execute
-//    ./tests/test_portsofcall "[!benchmark]"
-// (you may need to escape the `!` char depending
-// on your shell)
-TEST_CASE("Benchmark Indexing", "[!benchmark]") {
-  constexpr std::size_t NX = 16, NY = 16, NZ = 78;
-  constexpr std::size_t NU = 256, NV = 256, NW = 128;
-
-  // construct a series of PortableMDArrays
-  // of different rank and sizes, and
-  // iterates through all contiguous indexing.
-
-  Real a = 1.0;
-  // clang-format off
-  SECTION("small 1D") {
-    BENCHMARK_ADVANCED("baseline continguous")(Catch::Benchmark::Chronometer meter)
-    { 
-      Real *tape_d = (Real *)PORTABLE_MALLOC(sizeof(Real)*NX);
-
-      tape_d[0] = a;
-      meter.measure([=] {portableFor(
-        "fill 1D",
-        1, NX, 
-        PORTABLE_LAMBDA(auto i){tape_d[i]=a+tape_d[i-1];}); return 0;});
-
-      PORTABLE_FREE(tape_d);
-      return a;
-    };
-
-
-    BENCHMARK_ADVANCED("PMD continguous (space optimized)")(Catch::Benchmark::Chronometer meter)
-    { 
-      Real *tape_d = (Real *)PORTABLE_MALLOC(sizeof(Real)*NX);
-      auto view_d = PortableMDArray<Real, 1>(tape_d, NX);
-
-      view_d(0) = a;
-      meter.measure([=] {portableFor(
-        "fill 1D",
-        1, NX, 
-        PORTABLE_LAMBDA(auto i){view_d(i)=a+view_d(i-1);});return 0;});
-
-      PORTABLE_FREE(tape_d);
-    };
-
-    BENCHMARK_ADVANCED("PMD continguous (default size)")(Catch::Benchmark::Chronometer meter)
-    { 
-      Real *tape_d = (Real *)PORTABLE_MALLOC(sizeof(Real)*NX);
-      auto view_d = PortableMDArray<Real>(tape_d, NX);
-
-      view_d(0) = a;
-      meter.measure([=] {portableFor(
-        "fill 1D",
-        1, NX, 
-        PORTABLE_LAMBDA(auto i){view_d(i)=a+view_d(i-1);});return 0;});
-
-      auto b = view_d(0);
-      PORTABLE_FREE(tape_d);
-    };
-
-    BENCHMARK_ADVANCED("baseline shuffle")(Catch::Benchmark::Chronometer meter)
-    { 
-      Real *tape_d = (Real *)PORTABLE_MALLOC(sizeof(Real)*NX);
-      std::array<std::size_t, NX> offs;
-      std::iota(std::begin(offs), std::end(offs), 0);
-      std::shuffle(std::begin(offs), std::end(offs), std::mt19937{std::random_device{}()});
-
-      meter.measure([=] {portableFor(
-        "fill 1D",
-        0, NX, 
-        PORTABLE_LAMBDA(auto i){tape_d[offs[i]]=i;}); return 0;});
-
-      PORTABLE_FREE(tape_d);
-    };
-
-    BENCHMARK_ADVANCED("PMD shuffle (space optimized)")(Catch::Benchmark::Chronometer meter)
-    { 
-      Real *tape_d = (Real *)PORTABLE_MALLOC(sizeof(Real)*NX);
-      auto view_d = PortableMDArray<Real, 1>(tape_d, NX);
-
-      std::array<std::size_t, NX> offs;
-      std::iota(std::begin(offs), std::end(offs), 0);
-      std::shuffle(std::begin(offs), std::end(offs), std::mt19937{std::random_device{}()});
-
-      meter.measure([=] {portableFor(
-        "fill 1D",
-        1, NX, 
-        PORTABLE_LAMBDA(auto i){view_d(offs[i])=i;});return 0;});
-
-      PORTABLE_FREE(tape_d);
-    };
-
-
-  }
-  SECTION("small 2D") {
-    BENCHMARK_ADVANCED("baseline continguous")(Catch::Benchmark::Chronometer meter)
-    { 
-      Real *tape_d = (Real *)PORTABLE_MALLOC(sizeof(Real)*NX*NY);
-
-      portableFor("prep 2d", 0, NY, 0, 1, PORTABLE_LAMBDA(auto j, auto i) { tape_d[i+j*NX] = a;});
-      meter.measure([=] {portableFor(
-        "fill 1D",
-        0, NY, 1, NX,
-        PORTABLE_LAMBDA(auto j, auto i){tape_d[i + j * NX]=a + tape_d[(i-1) + j * NX];}); return 0;});
-
-      PORTABLE_FREE(tape_d);
-    };
-    BENCHMARK_ADVANCED("PMD contiguous (space optimized)")(Catch::Benchmark::Chronometer meter)
-    { 
-      Real *tape_d = (Real *)PORTABLE_MALLOC(sizeof(Real)*NX*NY);
-      auto view_d = PortableMDArray<Real, 2>(tape_d, NY, NX);
-
-      portableFor("prep 2d", 0, NY, 0, 1, PORTABLE_LAMBDA(auto j, auto i) { view_d(j, i) = a;});
-      view_d(0,0) = a;
-
-      meter.measure([=] {portableFor(
-        "fill 2D",
-        0, NY, 1, NX,
-        PORTABLE_LAMBDA(auto j, auto i){view_d(j,i)=a + view_d(j, i-1);});return 0;});
-
-      PORTABLE_FREE(tape_d);
-    };
-
-    BENCHMARK_ADVANCED("PMD continguous (default size)")(Catch::Benchmark::Chronometer meter)
-    { 
-      Real *tape_d = (Real *)PORTABLE_MALLOC(sizeof(Real)*NX*NY);
-      auto view_d = PortableMDArray<Real>(tape_d, NY, NX);
-
-      portableFor("prep 2d", 0, NY, 0, 1, PORTABLE_LAMBDA(auto j, auto i) { view_d(j, i) = a;});
-      view_d(0,0) = a;
-
-      meter.measure([=] {portableFor(
-        "fill 2D",
-        0, NY, 1, NX,
-        PORTABLE_LAMBDA(auto j, auto i){view_d(j,i)=a + view_d(j, i-1);});return 0;});
-
-      PORTABLE_FREE(tape_d);
-    };
-  }
-SECTION("small 5D"){
-
-    BENCHMARK_ADVANCED("baseline continguous")(Catch::Benchmark::Chronometer meter)
-    { 
-      Real *tape_d = (Real *)PORTABLE_MALLOC(sizeof(Real)*NX*NX*NY*NY*NZ);
-
-      portableFor(
-        "prep 5D",
-        0, NZ, 0, NY, 0, NY, 0, NX, 0, 1,
-        PORTABLE_LAMBDA(auto m, auto n, auto k, auto j, auto i){tape_d[i + j * NX + k * NX * NX + n * NX * NX * NY + m * NY * NY * NX * NX ]=a;});
-
-      meter.measure([=] {portableFor(
-        "fill 5D",
-        0, NZ, 0, NY, 0, NY, 0, NX, 1, NX,
-        PORTABLE_LAMBDA(auto m, auto n, auto k, auto j, auto i){tape_d[i + j * NX + k * NX * NX + n * NX * NX * NY + m * NY * NY * NX * NX ]=a + tape_d[i + j * NX + k * NX * NX + n * NX * NX * NY + m * NY * NY * NX * NX ];});return 0;});
-
-      PORTABLE_FREE(tape_d);
-    };
-
-    BENCHMARK_ADVANCED("PMD continguous")(Catch::Benchmark::Chronometer meter)
-    { 
-      Real *tape_d = (Real *)PORTABLE_MALLOC(sizeof(Real)*NX*NX*NY*NY*NZ);
-      auto view_d = PortableMDArray<Real, 5>(tape_d, NZ, NY, NY, NX, NX);
-
-      portableFor("prep 5d", 0, NZ, 0, NY, 0, NY, 0, NX, 0, 1, PORTABLE_LAMBDA(auto m, auto n, auto k, auto j, auto i) { view_d(m,n,k,j, i) = a;});
-      meter.measure([=] {portableFor(
-        "fill 5D",
-        0, NZ, 0, NY, 0, NY, 0, NX, 1, NX,
-        PORTABLE_LAMBDA(auto m, auto n, auto k, auto j, auto i){view_d(m,n,k,j,i)=a + view_d(m,n,k,j,i-1);});return 0;});
-
-      PORTABLE_FREE(tape_d);
-    };
-  }
-
-  // clang-format on
-}
diff --git a/test/test_utilities.hpp b/test/test_utilities.hpp
index 7d2023a1..c9e8f1cf 100644
--- a/test/test_utilities.hpp
+++ b/test/test_utilities.hpp
@@ -55,7 +55,7 @@ PORTABLE_FORCEINLINE_FUNCTION constexpr auto array_fix(const std::array<T, N> &a
 // this was written ad-hoc, and could be more general
 // (e.g. pass label, functor in lambda)
 template <class View, class Array, std::size_t... J>
-PORTABLE_INLINE_FUNCTION constexpr auto pf_invoke(View v, const Array &ext_par,
+PORTABLE_INLINE_FUNCTION constexpr void pf_invoke(View v, const Array &ext_par,
                                                   std::index_sequence<J...>) {
 
   auto stud_arr = array_fix(ext_par);

From 68bae3b3b8d740af45eb537071dfeca382ec6bfd Mon Sep 17 00:00:00 2001
From: Jonah Miller <jonahm@lanl.gov>
Date: Fri, 13 Dec 2024 09:36:16 -0700
Subject: [PATCH 18/22] fix compiler errors.

---
 ports-of-call/portable_arrays.hpp    | 18 ++++++++++--------
 ports-of-call/utility/array_algo.hpp | 10 ++++++----
 ports-of-call/utility/index_algo.hpp |  4 ++--
 test/test_mdidx.cpp                  |  9 +++++----
 test/test_portability.cpp            |  3 ++-
 5 files changed, 25 insertions(+), 19 deletions(-)

diff --git a/ports-of-call/portable_arrays.hpp b/ports-of-call/portable_arrays.hpp
index f3ba64c5..547b15ce 100644
--- a/ports-of-call/portable_arrays.hpp
+++ b/ports-of-call/portable_arrays.hpp
@@ -39,7 +39,7 @@
 
 namespace PortsOfCall {
 // maximum number of dimensions
-constexpr std::size_t MAXDIM = 6;
+constexpr std::size_t DEFAULT_MAXDIM = 6;
 
 namespace detail {
 // convert `index_sequence` values to constant
@@ -61,9 +61,10 @@ using Array = PortsOfCall::array<T, N>;
 template <auto N>
 using IArray = Array<std::size_t, N>;
 
-template <typename T, std::size_t D = MAXDIM>
+template <typename T, std::size_t D = DEFAULT_MAXDIM>
 class PortableMDArray {
  public:
+  static constexpr const std::size_t MAXDIM = D;
   using this_type = PortableMDArray<T, D>;
   using element_type = T;
   using value_type = typename std::remove_cv<T>::type;
@@ -89,14 +90,14 @@ class PortableMDArray {
   // ctors
   // default ctor: simply set null PortableMDArray
   PORTABLE_FUNCTION
-  PortableMDArray() noexcept : pdata_(nullptr), nxs_{{0}}, rank_{0} {}
+  PortableMDArray() noexcept : pdata_(nullptr), nxs_{{0}}, strides_{{0}}, rank_{0} {}
   // define copy constructor and overload assignment operator so both do deep
   // copies.
   //  PortableMDArray(const PortableMDArray<T> &t) noexcept;
   //  PortableMDArray<T> &operator=(const PortableMDArray<T> &t) noexcept;
   PortableMDArray(const this_type &src) noexcept
-      : nxs_(src.nxs_), rank_(src.rank_), strides_(src.strides_),
-        pdata_(src.pdata_ ? src.pdata_ : nullptr) {}
+      : pdata_(src.pdata_ ? src.pdata_ : nullptr), nxs_(src.nxs_), strides_(src.strides_),
+        rank_(src.rank_) {}
 
   PortableMDArray<T, D> &operator=(const this_type &src) noexcept {
     if (this != &src) {
@@ -254,17 +255,18 @@ class PortableMDArray {
   //  entries of the src array for d<dim (cannot access any nx4=2, etc. entries
   //  if dim=3 for example)
 
-  PORTABLE_FUNCTION void InitWithShallowSlice(const this_type &src, const int dim,
+  PORTABLE_FUNCTION void InitWithShallowSlice(const this_type &src, const size_t dim,
                                               const int indx, const int nvar) {
     pdata_ = src.pdata_;
     std::size_t offs = indx;
     nxs_[dim - 1] = nvar;
 
-    for (std::size_t i = 0; i < dim - 1; ++i) {
+    // this one needs to be an int because compared to input int
+    for (size_t i = 0; i < dim - 1; ++i) {
       nxs_[i] = src.nxs_[i];
       offs *= nxs_[i];
     }
-    for (std::size_t i = dim; i < MAXDIM; ++i) {
+    for (size_t i = dim; i < MAXDIM; ++i) {
       nxs_[i] = 1;
     }
     pdata_ += offs;
diff --git a/ports-of-call/utility/array_algo.hpp b/ports-of-call/utility/array_algo.hpp
index 800f94de..98945b4a 100644
--- a/ports-of-call/utility/array_algo.hpp
+++ b/ports-of-call/utility/array_algo.hpp
@@ -186,10 +186,11 @@ PORTABLE_FORCEINLINE_FUNCTION constexpr T array_reduce(A x, T initial_value, Op
 template <auto P, auto Fill = 1, template <class, auto> class A, class T, auto O>
 PORTABLE_FORCEINLINE_FUNCTION constexpr decltype(auto)
 make_underfilled_array(const A<T, O> &in) {
+  using i_t = typename A<T, 0>::size_type;
   A<T, P> out;
-  for (auto i = 0; i < in.size(); ++i)
+  for (i_t i = 0; i < in.size(); ++i)
     out[i] = in[i];
-  for (auto i = in.size(); i < out.size(); ++i)
+  for (i_t i = in.size(); i < out.size(); ++i)
     out[i] = Fill;
   return out;
 }
@@ -201,10 +202,11 @@ make_underfilled_array(const A<T, O> &in) {
 template <auto P, auto Fill = 1, template <class, auto> class A, class T, auto O>
 PORTABLE_FORCEINLINE_FUNCTION constexpr decltype(auto)
 make_underfilled_reversed_array(const A<T, O> &in) {
+  using i_t = typename A<T, 0>::size_type;
   A<T, P> out;
-  for (auto i = in.size() - 1; i >= 0; i--)
+  for (i_t i = in.size() - 1; i >= 0; i--)
     out[i] = in[i];
-  for (auto i = in.size(); i < out.size(); ++i)
+  for (i_t i = in.size(); i < out.size(); ++i)
     out[i] = Fill;
   return out;
 }
diff --git a/ports-of-call/utility/index_algo.hpp b/ports-of-call/utility/index_algo.hpp
index fe2796ac..cd5fbd9b 100644
--- a/ports-of-call/utility/index_algo.hpp
+++ b/ports-of-call/utility/index_algo.hpp
@@ -53,8 +53,8 @@ PORTABLE_FORCEINLINE_FUNCTION static constexpr auto get_strides(A const &dim) {
 // returns the flat (1D) index of the md index set {i,j,k}
 // NB: fast because the strides are provided and don't need to be recomputed
 template <class A>
-PORTABLE_FUNCTION static constexpr auto fast_findex(A const &ijk, A const &dim,
-                                                    A const &stride) {
+PORTABLE_FUNCTION static constexpr auto
+fast_findex(A const &ijk, [[maybe_unused]] A const &dim, A const &stride) {
   // TODO: assert ijk in bounds
   constexpr auto N = get_size(A{});
   if constexpr (N == 1) {
diff --git a/test/test_mdidx.cpp b/test/test_mdidx.cpp
index 108f6a20..7ef0d15b 100644
--- a/test/test_mdidx.cpp
+++ b/test/test_mdidx.cpp
@@ -28,6 +28,7 @@
 #endif
 
 using namespace PortsOfCall;
+constexpr const std::size_t MAXDIM = DEFAULT_MAXDIM;
 
 TEST_CASE("PortableMDArrays Sizes are sane", "[PortableMDArray]") {
   using tape_t = std::vector<int>;
@@ -38,7 +39,7 @@ TEST_CASE("PortableMDArrays Sizes are sane", "[PortableMDArray]") {
   std::vector<tape_t> dats;
 
   std::partial_sum(nxs.cbegin(), nxs.cend(), subsz.begin(), std::multiplies<int>());
-  for (auto i = 0; i < MAXDIM; ++i) {
+  for (std::size_t i = 0; i < MAXDIM; ++i) {
     dats.push_back(std::vector<int>(subsz[i], nxs[i]));
   }
 
@@ -70,7 +71,7 @@ TEST_CASE("Correct portable indexing", "[PortableMDArray]") {
   // vector length N on host of Real
   std::vector<Real> tape_ref(NC), tape_buf(NC);
 
-  for (auto n = 0; n < NC; ++n) {
+  for (std::size_t n = 0; n < NC; ++n) {
     tape_ref[n] = scale * static_cast<Real>(n);
   }
 
@@ -85,11 +86,11 @@ TEST_CASE("Correct portable indexing", "[PortableMDArray]") {
       });
 
   portableCopyToHost(tape_buf.data(), tape_d, NCb);
-  for (auto n = 0; n < NC; ++n) {
+  for (std::size_t n = 0; n < NC; ++n) {
     std::cout << "REF=" << tape_ref[n] << " BUF=" << tape_buf[n] << " n=" << n << "\n";
   }
 
-  for (auto n = 0; n < NC; ++n) {
+  for (std::size_t n = 0; n < NC; ++n) {
     INFO("REF=" << tape_ref[n] << " BUF=" << tape_buf[n] << " n=" << n);
     REQUIRE_THAT(tape_buf[n], Catch::Matchers::WithinRel(tape_ref[n]));
   }
diff --git a/test/test_portability.cpp b/test/test_portability.cpp
index 811f6322..b2542b4c 100644
--- a/test/test_portability.cpp
+++ b/test/test_portability.cpp
@@ -87,7 +87,8 @@ TEST_CASE("portableCopy works with all portability strategies", "[portableCopy]"
   Real *a = (Real *)PORTABLE_MALLOC(Nb);
 
   // set device values to 0
-  portableFor("set to 0", 0, N, PORTABLE_LAMBDA(const int &i) { a[i] = 0.0; });
+  portableFor(
+      "set to 0", 0, N, PORTABLE_LAMBDA(const int &i) { a[i] = 0.0; });
 
   // set host values to reference
   for (size_t i = 0; i < N; ++i) {

From d47b0756ebe05577b111dc1d54c43da5897ff811 Mon Sep 17 00:00:00 2001
From: Jonah Miller <jonahm@lanl.gov>
Date: Fri, 13 Dec 2024 09:41:36 -0700
Subject: [PATCH 19/22] return deprecated GetDim calls. Been long enough and
 Chris's change is already API breaking

---
 ports-of-call/portable_arrays.hpp | 25 -------------------------
 1 file changed, 25 deletions(-)

diff --git a/ports-of-call/portable_arrays.hpp b/ports-of-call/portable_arrays.hpp
index 547b15ce..73539608 100644
--- a/ports-of-call/portable_arrays.hpp
+++ b/ports-of-call/portable_arrays.hpp
@@ -130,31 +130,6 @@ class PortableMDArray {
     return nxs_[rank_ - I];
   }
 
-  // legacy API, TODO: deprecate
-  [[deprecated("Use GetDim<N> instead.")]] PORTABLE_FORCEINLINE_FUNCTION int
-  GetDim1() const {
-    return GetDim<1>();
-  }
-  [[deprecated("Use GetDim<N> instead.")]] PORTABLE_FORCEINLINE_FUNCTION int
-  GetDim2() const {
-    return GetDim<2>();
-  }
-  [[deprecated("Use GetDim<N> instead.")]] PORTABLE_FORCEINLINE_FUNCTION int
-  GetDim3() const {
-    return GetDim<3>();
-  }
-  [[deprecated("Use GetDim<N> instead.")]] PORTABLE_FORCEINLINE_FUNCTION int
-  GetDim4() const {
-    return GetDim<4>();
-  }
-  [[deprecated("Use GetDim<N> instead.")]] PORTABLE_FORCEINLINE_FUNCTION int
-  GetDim5() const {
-    return GetDim<5>();
-  }
-  [[deprecated("Use GetDim<N> instead.")]] PORTABLE_FORCEINLINE_FUNCTION int
-  GetDim6() const {
-    return GetDim<6>();
-  }
   PORTABLE_INLINE_FUNCTION int GetDim(size_t i) const {
     // TODO: remove if performance cirtical
     assert(0 < i && i <= 6 && "PortableMDArrays are max 6D");

From e180050b57f27bc38d8e43244df761225565dab0 Mon Sep 17 00:00:00 2001
From: Jonah Miller <jonahm@lanl.gov>
Date: Fri, 13 Dec 2024 10:01:24 -0700
Subject: [PATCH 20/22] fast index type punning

---
 ports-of-call/portable_arrays.hpp    | 20 ++++++++++++++++++--
 ports-of-call/utility/index_algo.hpp |  7 ++++++-
 2 files changed, 24 insertions(+), 3 deletions(-)

diff --git a/ports-of-call/portable_arrays.hpp b/ports-of-call/portable_arrays.hpp
index 73539608..ce9a550c 100644
--- a/ports-of-call/portable_arrays.hpp
+++ b/ports-of-call/portable_arrays.hpp
@@ -187,12 +187,28 @@ class PortableMDArray {
   // l-value, e.g.: a(3) = 3.0;
   template <typename... Is> //, class = std::enable_if_t<(sizeof...(Is) > 3)>>
   PORTABLE_FORCEINLINE_FUNCTION constexpr reference operator()(Is... idxs) noexcept {
-    return pdata_[util::fast_findex({static_cast<size_type>(idxs)...}, nxs_, strides_)];
+    index_type idx;
+    if constexpr (sizeof...(Is) == 0) {
+      idx = 0;
+    } else if constexpr (sizeof...(Is) == 1) {
+      idx = get_first(idxs...);
+    } else {
+      idx = util::fast_findex({static_cast<size_type>(idxs)...}, nxs_, strides_);
+    }
+    return pdata_[idx];
   }
 
   template <typename... Is> //, class = std::enable_if_t<(sizeof...(Is) > 3)>>
   PORTABLE_FORCEINLINE_FUNCTION constexpr T &operator()(Is... idxs) const noexcept {
-    return pdata_[util::fast_findex({static_cast<size_type>(idxs)...}, nxs_, strides_)];
+    index_type idx;
+    if constexpr (sizeof...(Is) == 0) {
+      idx = 0;
+    } else if constexpr (sizeof...(Is) == 1) {
+      idx = get_first(idxs...);
+    } else {
+      idx = util::fast_findex({static_cast<size_type>(idxs)...}, nxs_, strides_);
+    }
+    return pdata_[idx];
   }
 
   this_type &operator*=(T scale) {
diff --git a/ports-of-call/utility/index_algo.hpp b/ports-of-call/utility/index_algo.hpp
index cd5fbd9b..74b9b9da 100644
--- a/ports-of-call/utility/index_algo.hpp
+++ b/ports-of-call/utility/index_algo.hpp
@@ -50,10 +50,15 @@ PORTABLE_FORCEINLINE_FUNCTION static constexpr auto get_strides(A const &dim) {
   return detail::get_strides_impl(dim, std::make_index_sequence<get_size(A{})>{});
 }
 
+template <typename I, typename... Is>
+PORTABLE_FORCEINLINE_FUNCTION static constexpr auto get_first(I head, Is... tail) {
+  return head;
+}
+
 // returns the flat (1D) index of the md index set {i,j,k}
 // NB: fast because the strides are provided and don't need to be recomputed
 template <class A>
-PORTABLE_FUNCTION static constexpr auto
+PORTABLE_FORCEINLINE_FUNCTION static constexpr auto
 fast_findex(A const &ijk, [[maybe_unused]] A const &dim, A const &stride) {
   // TODO: assert ijk in bounds
   constexpr auto N = get_size(A{});

From b3aba8684fde0427ee4e0d53afb48659ab951dec Mon Sep 17 00:00:00 2001
From: Jonah Miller <jonahm@lanl.gov>
Date: Fri, 13 Dec 2024 10:03:39 -0700
Subject: [PATCH 21/22] namespace

---
 ports-of-call/portable_arrays.hpp | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/ports-of-call/portable_arrays.hpp b/ports-of-call/portable_arrays.hpp
index ce9a550c..1dec9312 100644
--- a/ports-of-call/portable_arrays.hpp
+++ b/ports-of-call/portable_arrays.hpp
@@ -191,7 +191,7 @@ class PortableMDArray {
     if constexpr (sizeof...(Is) == 0) {
       idx = 0;
     } else if constexpr (sizeof...(Is) == 1) {
-      idx = get_first(idxs...);
+      idx = util::get_first(idxs...);
     } else {
       idx = util::fast_findex({static_cast<size_type>(idxs)...}, nxs_, strides_);
     }
@@ -204,7 +204,7 @@ class PortableMDArray {
     if constexpr (sizeof...(Is) == 0) {
       idx = 0;
     } else if constexpr (sizeof...(Is) == 1) {
-      idx = get_first(idxs...);
+      idx = util::get_first(idxs...);
     } else {
       idx = util::fast_findex({static_cast<size_type>(idxs)...}, nxs_, strides_);
     }

From 3b5d9bae3a427704d91c8a04f3e6b78f3c9a1db3 Mon Sep 17 00:00:00 2001
From: Christopher Mauney <mauneyc@pn2305185.lanl.gov>
Date: Wed, 18 Dec 2024 07:15:33 -0700
Subject: [PATCH 22/22] null case for index

---
 ports-of-call/utility/index_algo.hpp | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

diff --git a/ports-of-call/utility/index_algo.hpp b/ports-of-call/utility/index_algo.hpp
index 74b9b9da..8b6d120d 100644
--- a/ports-of-call/utility/index_algo.hpp
+++ b/ports-of-call/utility/index_algo.hpp
@@ -62,7 +62,9 @@ PORTABLE_FORCEINLINE_FUNCTION static constexpr auto
 fast_findex(A const &ijk, [[maybe_unused]] A const &dim, A const &stride) {
   // TODO: assert ijk in bounds
   constexpr auto N = get_size(A{});
-  if constexpr (N == 1) {
+  if constexpr (N == 0) {
+    return 0;
+  } else if constexpr (N == 1) {
     return ijk[0];
   } else if constexpr (N == 2) {
     return ijk[1] + ijk[0] * stride[0];