diff --git a/doc/sphinx/src/using.rst b/doc/sphinx/src/using.rst
index 909611ea..40a9d7a5 100644
--- a/doc/sphinx/src/using.rst
+++ b/doc/sphinx/src/using.rst
@@ -220,3 +220,13 @@ data pointer. It accepts anywhere between 1 and 6 sizes.
 
 ``PortableMDArray`` also supports some simple boolean comparitors,
 such as ``==`` and arithmetic such as ``+``, and ``-``.
+
+array.hpp
+^^^^^^^^^
+
+``PortsOfCall::array`` is intended to be a drop-in replacement for ``std::array``, with the
+exception that it works on GPUs.  As of C++17, ``std::array::fill`` and ``std::array::swap`` are
+not yet ``constexpr``, so even with the "relaxed ``constexpr``" compilation mode ``std::array`` is
+not feature-complete on GPUs.  This will change when those member functions become ``constexpr`` in
+C++20.
+
diff --git a/ports-of-call/array.hpp b/ports-of-call/array.hpp
new file mode 100644
index 00000000..8eb90cc0
--- /dev/null
+++ b/ports-of-call/array.hpp
@@ -0,0 +1,270 @@
+#ifndef _PORTS_OF_CALL_ARRAY_HPP_
+#define _PORTS_OF_CALL_ARRAY_HPP_
+
+#include "portability.hpp"
+#include "portable_errors.hpp"
+
+#include <cstddef>
+#include <iterator>
+#include <tuple>
+
+namespace PortsOfCall {
+
+//! The array class is a constexpr (mostly) replacement for std::array
+/*!
+std::array is not fully constexpr until C++20.  Currently we are relying on
+"constexpr"-ness to offload to GPUs.  PortsOfCall::array provides a
+constexpr replacement for std::array with the following notable exceptions:
+    - the reverse iterators are not provided because std::reverse_iterator is
+not constexpr until C++17 (could be provided with modest effort)
+
+    - the `at` function is not provided because it bounds checks and throws an
+exception which we cannot do on a GPU (it could be added and made
+non-constexpr)
+
+    - no std::get (cannot be provided without undefined behavior)
+*/
+
+namespace array_detail {
+
+template <typename T, std::size_t N>
+struct bits {
+  using type = T[N];
+
+  template <typename IntT>
+  PORTABLE_FORCEINLINE_FUNCTION static constexpr T &ref_access(type &array, IntT index) {
+    return array[index];
+  }
+
+  template <typename IntT>
+  PORTABLE_FORCEINLINE_FUNCTION static constexpr T const &ref_access(type const &array,
+                                                                     IntT index) {
+    return array[index];
+  }
+
+  PORTABLE_FORCEINLINE_FUNCTION static constexpr T *ptr_access(type &array) {
+    return static_cast<T *>(array);
+  }
+
+  PORTABLE_FORCEINLINE_FUNCTION static constexpr T const *ptr_access(type const &array) {
+    return static_cast<T const *>(array);
+  }
+};
+
+template <typename T>
+struct bits<T, 0> {
+  using type = T[1];
+
+  template <typename IntT>
+  PORTABLE_FORCEINLINE_FUNCTION static constexpr T &ref_access(type &arr, IntT) {
+    return arr[0];
+  }
+
+  template <typename IntT>
+  PORTABLE_FORCEINLINE_FUNCTION static constexpr T const &ref_access(type const &arr,
+                                                                     IntT) {
+    return arr[0];
+  }
+
+  PORTABLE_FORCEINLINE_FUNCTION static constexpr T *ptr_access(type &arr) {
+    return (T *)&arr;
+  }
+
+  PORTABLE_FORCEINLINE_FUNCTION static constexpr T const *ptr_access(type const &arr) {
+    return (T const *)&arr;
+  }
+};
+} // namespace array_detail
+
+template <typename T, std::size_t N>
+struct array {
+  // member types
+  using value_type = T;
+  using size_type = std::size_t;
+  using difference_type = std::ptrdiff_t;
+  using reference = T &;
+  using const_reference = T const &;
+  using pointer = value_type *;
+  using const_pointer = value_type const *;
+  using iterator = pointer;
+  using const_iterator = const_pointer;
+
+  //! The  array being wrapped by this class
+  typename array_detail::bits<T, N>::type arr{};
+
+  //! Index operator that returns a reference to the desired element.
+  PORTABLE_FUNCTION constexpr reference operator[](size_type index) {
+    return array_detail::bits<T, N>::ref_access(arr, index);
+  }
+
+  //! bounds-checked equivalent to index operator
+  PORTABLE_FUNCTION constexpr reference at(size_type index) {
+    if (index >= size()) {
+      PORTABLE_ALWAYS_THROW_OR_ABORT("invalid index.");
+    }
+    return array_detail::bits<T, N>::ref_access(arr, index);
+  }
+
+  //! Index operator that returns a const reference to the desired element.
+  PORTABLE_FUNCTION constexpr const_reference operator[](size_type index) const {
+    return array_detail::bits<T, N>::ref_access(arr, index);
+  }
+
+  //! bounds-checked equivalent to index operator
+  PORTABLE_FUNCTION constexpr const_reference at(size_type index) const {
+    if (index >= size()) {
+      PORTABLE_ALWAYS_THROW_OR_ABORT("invalid index.");
+    }
+    return array_detail::bits<T, N>::ref_access(arr, index);
+  }
+
+  //! access first element of array
+  PORTABLE_FUNCTION constexpr reference front() {
+    return array_detail::bits<T, N>::ref_access(arr, 0);
+  }
+
+  //! access first element of const array
+  PORTABLE_FUNCTION constexpr const_reference front() const {
+    return array_detail::bits<T, N>::ref_access(arr, 0);
+  }
+
+  //! access last element of array
+  PORTABLE_FUNCTION constexpr reference back() {
+    return N ? array_detail::bits<T, N>::ref_access(arr, N - 1)
+             : array_detail::bits<T, N>::ref_access(arr, 0);
+  }
+
+  //! access first element of const array
+  PORTABLE_FUNCTION constexpr const_reference back() const {
+    return N ? array_detail::bits<T, N>::ref_access(arr, N - 1)
+             : array_detail::bits<T, N>::ref_access(arr, 0);
+  }
+
+  //! get pointer to underlying data of array
+  PORTABLE_FUNCTION constexpr pointer data() {
+    return array_detail::bits<T, N>::ptr_access(arr);
+  }
+
+  //! get pointer to underlying data of const array
+  PORTABLE_FUNCTION constexpr const_pointer data() const {
+    return array_detail::bits<T, N>::ptr_access(arr);
+  }
+
+  //! Provides an iterator pointing to the beginning of the array.
+  PORTABLE_FUNCTION constexpr iterator begin() { return iterator{data()}; }
+
+  //! Provides a const iterator pointing to the beginning of the const array.
+  PORTABLE_FUNCTION constexpr const_iterator begin() const {
+    return const_iterator{data()};
+  }
+
+  //! Provides a const iterator pointing to the beginning of the const array.
+  PORTABLE_FUNCTION constexpr const_iterator cbegin() const {
+    return const_iterator{data()};
+  }
+
+  //! Provides an iterator pointing to the end of the array.
+  PORTABLE_FUNCTION constexpr iterator end() { return iterator{data() + N}; }
+
+  //! Provides a const iterator pointing to the end of the const array.
+  PORTABLE_FUNCTION constexpr const_iterator end() const {
+    return const_iterator{data() + N};
+  }
+
+  //! Provides a const iterator pointing to the end of the const array.
+  PORTABLE_FUNCTION constexpr const_iterator cend() const {
+    return const_iterator{data() + N};
+  }
+
+  //! test if array is empty, i.e., N==0
+  PORTABLE_FUNCTION constexpr bool empty() const { return N == 0; }
+
+  //! return the size of an array
+  PORTABLE_FUNCTION constexpr size_type size() const { return N; }
+
+  //! return the maximum size an array can hold
+  PORTABLE_FUNCTION constexpr size_type max_size() const { return N; }
+
+  //! fill the array with a single value
+  PORTABLE_FUNCTION constexpr void fill(const_reference value) {
+    for (auto &element : *this) {
+      element = value;
+    }
+  }
+
+  //! swap the array contents with another
+  PORTABLE_FUNCTION constexpr void swap(array &other) {
+    using std::swap;
+    for (size_type i = 0; i < N; ++i) {
+      swap(arr[i], other.arr[i]);
+    }
+  }
+};
+
+//! specialization of swap for array
+template <typename T, std::size_t N>
+PORTABLE_FUNCTION constexpr void swap(array<T, N> &left, array<T, N> &right) {
+  left.swap(right);
+}
+
+//! make_array function to return a deduced type array
+// This function was proposed for standardization but superceded by
+// C++17 automatic template parameter deduction.  It can still be
+// useful at C++14 and earlier.
+// * This implementation does not support the special handling of
+//   of reference wrapper
+
+namespace array_detail {
+
+// determine the element type of the array
+// * if explicitly specified D
+// * if not explicitly specified std::common_type_t<ArgTs...>
+
+template <typename D, typename... ArgTs>
+struct array_element {
+  using type = D;
+};
+
+template <typename... ArgTs>
+struct array_element<void, ArgTs...> : std::common_type<ArgTs...> {};
+
+} // namespace array_detail
+
+template <typename D = void, typename... ArgTs>
+PORTABLE_FUNCTION constexpr array<typename array_detail::array_element<D, ArgTs...>::type,
+                                  sizeof...(ArgTs)>
+make_array(ArgTs &&...args) {
+  return {{std::forward<ArgTs>(args)...}};
+}
+
+} // end namespace PortsOfCall
+
+// provide std specializations
+namespace std {
+/* libc++ and libstdc++ define tuple_size/tuple_element as a class and struct,
+ * respectively.  To prevent clang from issuing a warning that these
+ * customization points are mismatched, we are temporarily disabling that
+ * diagnostic. */
+#if defined(__clang__)
+#pragma clang diagnostic push
+#pragma clang diagnostic ignored "-Wmismatched-tags"
+#endif
+
+template <typename T, std::size_t N>
+struct tuple_size<PortsOfCall::array<T, N>>
+    : public std::integral_constant<std::size_t, N> {};
+
+template <std::size_t I, typename T, std::size_t N>
+struct tuple_element<I, PortsOfCall::array<T, N>> {
+ public:
+  static_assert(I < N, "index must be less than number of elements");
+
+  using type = T;
+};
+
+#if defined(__clang__)
+#pragma clang diagnostic pop
+#endif
+} // namespace std
+
+#endif // #ifndef _PORTS_OF_CALL_ARRAY_HPP_
diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
index e59af456..8e83f110 100644
--- a/test/CMakeLists.txt
+++ b/test/CMakeLists.txt
@@ -51,4 +51,4 @@ target_link_libraries(test_portsofcall
 include(Catch)
 catch_discover_tests(test_portsofcall)
 
-target_sources(test_portsofcall PRIVATE test_portability.cpp)
+target_sources(test_portsofcall PRIVATE test_portability.cpp test_array.cpp)
diff --git a/test/test_array.cpp b/test/test_array.cpp
new file mode 100644
index 00000000..3a9423fc
--- /dev/null
+++ b/test/test_array.cpp
@@ -0,0 +1,278 @@
+#include "ports-of-call/array.hpp"
+#include "ports-of-call/portability.hpp"
+
+#ifndef CATCH_CONFIG_FAST_COMPILE
+#define CATCH_CONFIG_FAST_COMPILE
+#include <catch2/catch_test_macros.hpp>
+#endif
+
+#include <iostream>
+
+using namespace PortsOfCall;
+
+TEST_CASE("array nominal element access (GPU)", "[array][GPU]") {
+  // Declare the array
+  constexpr int N = 16;
+  array<double, N> arr;
+
+  // Fill the array
+  for (int i = 0; i < N; ++i) {
+    arr[i] = i + 1;
+  }
+
+  // Can we read from it on the GPU?
+  int count = 0;
+  auto func = PORTABLE_LAMBDA(const int i, int &count) mutable {
+    if (arr[i] = i + 1) {
+      ++count;
+    }
+  };
+  portableReduce("assign_and_check", 0, N, func, count);
+  CHECK(count == N);
+}
+
+TEST_CASE("array nominal element access", "[array]") {
+  // Declare the array
+  constexpr int N = 16;
+  array<double, N> arr;
+
+  // Fill the array
+  for (int i = 0; i < N; ++i) {
+    arr[i] = i + 1;
+  }
+
+  // Can we read from it?
+  for (int i = 0; i < N; ++i) {
+    CHECK(arr[i] == i + 1);
+  }
+}
+
+TEST_CASE("array zero-sized element access", "[array]") {
+  array<double, 0> arr;
+
+  // accessing the zeroth element of a zero-size array is
+  // allowed but undefined behavior.  We will only verify that
+  // the address returned is well defined.
+  CHECK(std::addressof(arr[0]));
+}
+
+TEST_CASE("const array element access", "[array]") {
+  // Declare the array
+  constexpr int N = 5;
+  array<int, N> const arr{{0, 1, 2, 3, 4}};
+
+  // Can we read from it?
+  for (int i = 0; i < N; ++i) {
+    CHECK(arr[i] == i);
+  }
+}
+
+TEST_CASE("const array zero-sized element access", "[array]") {
+  array<double, 0> const arr{};
+
+  // accessing the zeroth element of a zero-size array is
+  // allowed but undefined behavior.  We will only verify that
+  // the address returned is well defined.
+  CHECK(std::addressof(arr[0]));
+}
+
+TEST_CASE("array range-based for loop", "[array]") {
+  SECTION("with non-const array") {
+    // Declare the array
+    constexpr int N = 15;
+    array<int, N> arr;
+
+    // Fill the loop
+    int i = 0;
+    for (auto &x : arr) {
+      x = ++i;
+    }
+
+    // Check the values
+    i = 0;
+    for (auto &x : arr) {
+      ++i;
+      CHECK(x == i);
+    }
+  }
+
+  SECTION("with const array") {
+    // Declare the array
+    constexpr int N = 5;
+    array<int, N> const arr{{1, 2, 3, 4, 5}};
+
+    // Check the values
+    int i = 0;
+    for (auto const &x : arr) {
+      ++i;
+      CHECK(x == i);
+    }
+  }
+}
+
+TEST_CASE("array begins and ends", "[array]") {
+  using std::begin;
+  using std::cbegin;
+  using std::cend;
+  using std::end;
+
+  SECTION("with non-const array") {
+    array<int, 3> arr{{1, 2, 3}};
+
+    CHECK(std::is_same<int *, decltype(begin(arr))>::value);
+    CHECK(std::is_same<int const *, decltype(cbegin(arr))>::value);
+
+    CHECK(std::is_same<int *, decltype(end(arr))>::value);
+    CHECK(std::is_same<int const *, decltype(cend(arr))>::value);
+  }
+
+  SECTION("with const array") {
+    array<int, 3> const arr{{1, 2, 3}};
+
+    CHECK(std::is_same<int const *, decltype(begin(arr))>::value);
+    CHECK(std::is_same<int const *, decltype(cbegin(arr))>::value);
+
+    CHECK(std::is_same<int const *, decltype(end(arr))>::value);
+    CHECK(std::is_same<int const *, decltype(cend(arr))>::value);
+  }
+
+  SECTION("with non-const zero-sized array") {
+    array<int, 0> arr;
+
+    CHECK(begin(arr) == end(arr));
+  }
+
+  SECTION("with const zero-sized array") {
+    array<int, 0> const arr{};
+
+    CHECK(begin(arr) == end(arr));
+  }
+}
+
+TEST_CASE("array front and back", "[array]") {
+  array<int, 3> arr1{{3, 2, 1}};
+
+  CHECK(std::is_same<int &, decltype(arr1.front())>::value);
+  CHECK(arr1.front() == 3);
+
+  CHECK(std::is_same<int &, decltype(arr1.back())>::value);
+  CHECK(arr1.back() == 1);
+
+  array<double, 3> const arr2{{3, 2, 1}};
+
+  CHECK(std::is_same<double const &, decltype(arr2.front())>::value);
+  CHECK(arr2.front() == 3.0);
+
+  CHECK(std::is_same<double const &, decltype(arr2.back())>::value);
+  CHECK(arr2.back() == 1.0);
+}
+
+TEST_CASE("array data", "[array]") {
+  SECTION("with non-const array") {
+    array<int, 3> arr{{3, 2, 1}};
+
+    CHECK(std::is_same<int *, decltype(arr.data())>::value);
+    CHECK(arr.data() == std::addressof(arr[0]));
+  }
+
+  SECTION("with non-const array") {
+    array<int, 3> const arr{{3, 2, 1}};
+
+    CHECK(std::is_same<int const *, decltype(arr.data())>::value);
+    CHECK(arr.data() == std::addressof(arr[0]));
+  }
+}
+
+TEST_CASE("array empty", "[array]") {
+  SECTION("with non-empty array") {
+    array<int, 10> arr;
+
+    CHECK(not arr.empty());
+  }
+
+  SECTION("with empty array") {
+    array<int, 0> arr;
+
+    CHECK(arr.empty());
+  }
+}
+
+TEST_CASE("array sizes", "[array]") {
+  array<int, 42> arr;
+
+  CHECK(arr.size() == arr.max_size());
+}
+
+TEST_CASE("array fill (GPU)", "[array][GPU]") {
+  constexpr std::size_t N = 42;
+  std::size_t count = 0;
+  auto func = PORTABLE_LAMBDA(const int i, std::size_t &count) {
+    constexpr double value = 3.14;
+    array<double, N> arr;
+    arr.fill(value);
+    for (const double x : arr) {
+      count += (x == value) ? 1 : 0;
+    }
+  };
+  constexpr std::size_t M = 5;
+  portableReduce("check", 0, M, func, count);
+  CHECK(count == N * M);
+}
+
+TEST_CASE("array fill", "[array]") {
+  array<double, 42> arr;
+  arr.fill(3.14);
+
+  for (auto const &x : arr) {
+    CHECK(x == 3.14);
+  }
+}
+
+TEST_CASE("array swap", "[array]") {
+  array<int, 10> zeros;
+  zeros.fill(0);
+  array<int, 10> ones;
+  ones.fill(1);
+
+  zeros.swap(ones);
+
+  for (auto const &x : zeros) {
+    CHECK(x == 1);
+  }
+
+  for (auto const &x : ones) {
+    CHECK(x == 0);
+  }
+
+  using std::swap;
+
+  swap(zeros, ones);
+
+  for (auto const &x : zeros) {
+    CHECK(x == 0);
+  }
+
+  for (auto const &x : ones) {
+    CHECK(x == 1);
+  }
+}
+
+TEST_CASE("array tuple_size", "[array]") {
+  array<double, 5> arr;
+
+  CHECK(std::tuple_size<decltype(arr)>::value == 5);
+}
+
+TEST_CASE("array tuple_element", "[array]") {
+  struct Foo {};
+
+  array<Foo, 3> foos;
+
+  CHECK(std::is_same<Foo, std::tuple_element<1, decltype(foos)>::type>::value);
+}
+
+TEST_CASE("make_array function", "[array]") {
+  auto arr = make_array(1.0, 2.0, 3.0);
+
+  CHECK(std::is_same<decltype(arr), array<double, 3>>::value);
+}