Add GPU scan support for Array<T>::PartialSum

Remove explicit template instantiations

examples/ex1p_vo.cpp

@@ -1,348 +0,0 @@
-//                       MFEM Example 1 - Parallel Version
-//
-// Compile with: make ex1p
-//
-// Sample runs:  mpirun -np 4 ex1p -m ../data/square-disc.mesh
-//               mpirun -np 4 ex1p -m ../data/star.mesh
-//               mpirun -np 4 ex1p -m ../data/star-mixed.mesh
-//               mpirun -np 4 ex1p -m ../data/escher.mesh
-//               mpirun -np 4 ex1p -m ../data/fichera.mesh
-//               mpirun -np 4 ex1p -m ../data/fichera-mixed.mesh
-//               mpirun -np 4 ex1p -m ../data/toroid-wedge.mesh
-//               mpirun -np 4 ex1p -m ../data/octahedron.mesh -o 1
-//               mpirun -np 4 ex1p -m ../data/periodic-annulus-sector.msh
-//               mpirun -np 4 ex1p -m ../data/periodic-torus-sector.msh
-//               mpirun -np 4 ex1p -m ../data/square-disc-p2.vtk -o 2
-//               mpirun -np 4 ex1p -m ../data/square-disc-p3.mesh -o 3
-//               mpirun -np 4 ex1p -m ../data/square-disc-nurbs.mesh -o -1
-//               mpirun -np 4 ex1p -m ../data/star-mixed-p2.mesh -o 2
-//               mpirun -np 4 ex1p -m ../data/disc-nurbs.mesh -o -1
-//               mpirun -np 4 ex1p -m ../data/pipe-nurbs.mesh -o -1
-//               mpirun -np 4 ex1p -m ../data/ball-nurbs.mesh -o 2
-//               mpirun -np 4 ex1p -m ../data/fichera-mixed-p2.mesh -o 2
-//               mpirun -np 4 ex1p -m ../data/star-surf.mesh
-//               mpirun -np 4 ex1p -m ../data/square-disc-surf.mesh
-//               mpirun -np 4 ex1p -m ../data/inline-segment.mesh
-//               mpirun -np 4 ex1p -m ../data/amr-quad.mesh
-//               mpirun -np 4 ex1p -m ../data/amr-hex.mesh
-//               mpirun -np 4 ex1p -m ../data/mobius-strip.mesh
-//               mpirun -np 4 ex1p -m ../data/mobius-strip.mesh -o -1 -sc
-//
-// Device sample runs:
-//               mpirun -np 4 ex1p -pa -d cuda
-//               mpirun -np 4 ex1p -pa -d occa-cuda
-//               mpirun -np 4 ex1p -pa -d raja-omp
-//               mpirun -np 4 ex1p -pa -d ceed-cpu
-//               mpirun -np 4 ex1p -pa -d ceed-cpu -o 4 -a
-//             * mpirun -np 4 ex1p -pa -d ceed-cuda
-//             * mpirun -np 4 ex1p -pa -d ceed-hip
-//               mpirun -np 4 ex1p -pa -d ceed-cuda:/gpu/cuda/shared
-//               mpirun -np 4 ex1p -m ../data/beam-tet.mesh -pa -d ceed-cpu
-//
-// Description:  This example code demonstrates the use of MFEM to define a
-//               simple finite element discretization of the Laplace problem
-//               -Delta u = 1 with homogeneous Dirichlet boundary conditions.
-//               Specifically, we discretize using a FE space of the specified
-//               order, or if order < 1 using an isoparametric/isogeometric
-//               space (i.e. quadratic for quadratic curvilinear mesh, NURBS for
-//               NURBS mesh, etc.)
-//
-//               The example highlights the use of mesh refinement, finite
-//               element grid functions, as well as linear and bilinear forms
-//               corresponding to the left-hand side and right-hand side of the
-//               discrete linear system. We also cover the explicit elimination
-//               of essential boundary conditions, static condensation, and the
-//               optional connection to the GLVis tool for visualization.
-
-#include "mfem.hpp"
-#include "linalg/vector_operator.hpp"
-#include <fstream>
-#include <iostream>
-
-using namespace std;
-using namespace mfem;
-
-class CoordCoefficient : public Coefficient
-{
-private:
-   int d;
-
-   mutable Vector x;
-
-public:
-   CoordCoefficient(int d) : d(d), x(3) {}
-
-   double Eval(ElementTransformation &T, const IntegrationPoint &ip)
-   {
-      if (d == -1) { return 1.0; }
-
-      T.Transform(ip, x);
-      return x[d];
-   }
-};
-
-int main(int argc, char *argv[])
-{
-   // 1. Initialize MPI.
-   MPI_Session mpi;
-   int num_procs = mpi.WorldSize();
-   int myid = mpi.WorldRank();
-
-   // 2. Parse command-line options.
-   const char *mesh_file = "../data/star.mesh";
-   int order = 1;
-   bool static_cond = false;
-   bool pa = false;
-   const char *device_config = "cpu";
-   bool visualization = true;
-   bool algebraic_ceed = false;
-
-   OptionsParser args(argc, argv);
-   args.AddOption(&mesh_file, "-m", "--mesh",
-                  "Mesh file to use.");
-   args.AddOption(&order, "-o", "--order",
-                  "Finite element order (polynomial degree) or -1 for"
-                  " isoparametric space.");
-   args.AddOption(&static_cond, "-sc", "--static-condensation", "-no-sc",
-                  "--no-static-condensation", "Enable static condensation.");
-   args.AddOption(&pa, "-pa", "--partial-assembly", "-no-pa",
-                  "--no-partial-assembly", "Enable Partial Assembly.");
-   args.AddOption(&device_config, "-d", "--device",
-                  "Device configuration string, see Device::Configure().");
-#ifdef MFEM_USE_CEED
-   args.AddOption(&algebraic_ceed, "-a", "--algebraic",
-                  "-no-a", "--no-algebraic",
-                  "Use algebraic Ceed solver");
-#endif
-   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
-                  "--no-visualization",
-                  "Enable or disable GLVis visualization.");
-   args.Parse();
-   if (!args.Good())
-   {
-      if (myid == 0)
-      {
-         args.PrintUsage(cout);
-      }
-      return 1;
-   }
-   if (myid == 0)
-   {
-      args.PrintOptions(cout);
-   }
-
-   // 3. Enable hardware devices such as GPUs, and programming models such as
-   //    CUDA, OCCA, RAJA and OpenMP based on command line options.
-   Device device(device_config);
-   if (myid == 0) { device.Print(); }
-
-   // 4. Read the (serial) mesh from the given mesh file on all processors.  We
-   //    can handle triangular, quadrilateral, tetrahedral, hexahedral, surface
-   //    and volume meshes with the same code.
-   Mesh mesh(mesh_file, 1, 1);
-   int dim = mesh.Dimension();
-
-   // 5. Refine the serial mesh on all processors to increase the resolution. In
-   //    this example we do 'ref_levels' of uniform refinement. We choose
-   //    'ref_levels' to be the largest number that gives a final mesh with no
-   //    more than 10,000 elements.
-   {
-      int ref_levels =
-         (int)floor(log(10000./mesh.GetNE())/log(2.)/dim);
-      for (int l = 0; l < ref_levels; l++)
-      {
-         mesh.UniformRefinement();
-      }
-   }
-
-   // 6. Define a parallel mesh by a partitioning of the serial mesh. Refine
-   //    this mesh further in parallel to increase the resolution. Once the
-   //    parallel mesh is defined, the serial mesh can be deleted.
-   ParMesh pmesh(MPI_COMM_WORLD, mesh);
-   mesh.Clear();
-   {
-      int par_ref_levels = 2;
-      for (int l = 0; l < par_ref_levels; l++)
-      {
-         pmesh.UniformRefinement();
-      }
-   }
-
-   // 7. Define a parallel finite element space on the parallel mesh. Here we
-   //    use continuous Lagrange finite elements of the specified order. If
-   //    order < 1, we instead use an isoparametric/isogeometric space.
-   FiniteElementCollection *fec;
-   bool delete_fec;
-   if (order > 0)
-   {
-      fec = new H1_FECollection(order, dim);
-      delete_fec = true;
-   }
-   else if (pmesh.GetNodes())
-   {
-      fec = pmesh.GetNodes()->OwnFEC();
-      delete_fec = false;
-      if (myid == 0)
-      {
-         cout << "Using isoparametric FEs: " << fec->Name() << endl;
-      }
-   }
-   else
-   {
-      fec = new H1_FECollection(order = 1, dim);
-      delete_fec = true;
-   }
-   ParFiniteElementSpace fespace(&pmesh, fec);
-   HYPRE_BigInt size = fespace.GlobalTrueVSize();
-   if (myid == 0)
-   {
-      cout << "Number of finite element unknowns: " << size << endl;
-   }
-
-   // 8. Determine the list of true (i.e. parallel conforming) essential
-   //    boundary dofs. In this example, the boundary conditions are defined
-   //    by marking all the boundary attributes from the mesh as essential
-   //    (Dirichlet) and converting them to a list of true dofs.
-   Array<int> ess_tdof_list;
-   if (pmesh.bdr_attributes.Size())
-   {
-      Array<int> ess_bdr(pmesh.bdr_attributes.Max());
-      ess_bdr = 1;
-      fespace.GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
-   }
-
-   // 9. Set up the parallel linear form b(.) which corresponds to the
-   //    right-hand side of the FEM linear system, which in this case is
-   //    (1,phi_i) where phi_i are the basis functions in fespace.
-   ParLinearForm b(&fespace);
-   ConstantCoefficient one(1.0);
-   b.AddDomainIntegrator(new DomainLFIntegrator(one));
-   b.Assemble();
-
-   // 10. Define the solution vector x as a parallel finite element grid
-   //     function corresponding to fespace. Initialize x with initial guess of
-   //     zero, which satisfies the boundary conditions.
-   ParGridFunction x(&fespace);
-   x = 0.0;
-
-   ParVectorOperator vo(MPI_COMM_WORLD, myid, fespace.TrueVSize(), dim + 1);
-   {
-      for (int d=0; d <= dim; d++)
-      {
-         ParLinearForm bd(&fespace);
-         CoordCoefficient dCoef(d - 1);
-         bd.AddDomainIntegrator(new DomainLFIntegrator(dCoef));
-         bd.Assemble();
-
-         Vector *dv = new Vector(fespace.TrueVSize());
-         bd.ParallelAssemble(*dv);
-
-         vo.SetVector(d, dv, 1.0, true);
-      }
-   }
-
-   // 11. Set up the parallel bilinear form a(.,.) on the finite element space
-   //     corresponding to the Laplacian operator -Delta, by adding the
-   //     Diffusion domain integrator.
-   ParBilinearForm a(&fespace);
-   if (pa) { a.SetAssemblyLevel(AssemblyLevel::PARTIAL); }
-   a.AddDomainIntegrator(new DiffusionIntegrator(one));
-
-   // 12. Assemble the parallel bilinear form and the corresponding linear
-   //     system, applying any necessary transformations such as: parallel
-   //     assembly, eliminating boundary conditions, applying conforming
-   //     constraints for non-conforming AMR, static condensation, etc.
-   if (static_cond) { a.EnableStaticCondensation(); }
-   a.Assemble();
-
-   OperatorPtr A;
-   Vector B, X;
-   a.FormLinearSystem(ess_tdof_list, x, b, A, X, B);
-
-   // 13. Solve the linear system A X = B.
-   //     * With full assembly, use the BoomerAMG preconditioner from hypre.
-   //     * With partial assembly, use Jacobi smoothing, for now.
-   Solver *prec = NULL;
-   if (pa)
-   {
-      if (UsesTensorBasis(fespace))
-      {
-         if (algebraic_ceed)
-         {
-            prec = new ceed::AlgebraicSolver(a, ess_tdof_list);
-         }
-         else
-         {
-            prec = new OperatorJacobiSmoother(a, ess_tdof_list);
-         }
-      }
-   }
-   else
-   {
-      prec = new HypreBoomerAMG;
-   }
-   CGSolver cg(MPI_COMM_WORLD);
-   cg.SetRelTol(1e-12);
-   cg.SetMaxIter(2000);
-   cg.SetPrintLevel(1);
-   if (prec) { cg.SetPreconditioner(*prec); }
-   cg.SetOperator(*A);
-   cg.Mult(B, X);
-   delete prec;
-
-   {
-      Vector com((myid == 0) ? dim+1 : 0);
-      vo.Mult(X, com);
-
-      if (myid == 0)
-      {
-         cout << "Mass: " << com[0] << endl;
-         cout << "Center of mass: (";
-         for (int d=1; d<=dim; d++)
-         {
-            cout << com[d]/com[0];
-            if (d < dim) { cout << " ,"; }
-         }
-         cout << ")" << endl;
-      }
-   }
-
-   // 14. Recover the parallel grid function corresponding to X. This is the
-   //     local finite element solution on each processor.
-   a.RecoverFEMSolution(X, b, x);
-
-   // 15. Save the refined mesh and the solution in parallel. This output can
-   //     be viewed later using GLVis: "glvis -np <np> -m mesh -g sol".
-   {
-      ostringstream mesh_name, sol_name;
-      mesh_name << "mesh." << setfill('0') << setw(6) << myid;
-      sol_name << "sol." << setfill('0') << setw(6) << myid;
-
-      ofstream mesh_ofs(mesh_name.str().c_str());
-      mesh_ofs.precision(8);
-      pmesh.Print(mesh_ofs);
-
-      ofstream sol_ofs(sol_name.str().c_str());
-      sol_ofs.precision(8);
-      x.Save(sol_ofs);
-   }
-
-   // 16. Send the solution by socket to a GLVis server.
-   if (visualization)
-   {
-      char vishost[] = "localhost";
-      int  visport   = 19916;
-      socketstream sol_sock(vishost, visport);
-      sol_sock << "parallel " << num_procs << " " << myid << "\n";
-      sol_sock.precision(8);
-      sol_sock << "solution\n" << pmesh << x << flush;
-   }
-
-   // 17. Free the used memory.
-   if (delete_fec)
-   {
-      delete fec;
-   }
-
-   return 0;
-}

general/array.cpp

@@ -1,214 +0,0 @@
-// Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
-// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
-// LICENSE and NOTICE for details. LLNL-CODE-806117.
-//
-// This file is part of the MFEM library. For more information and source code
-// availability visit https://mfem.org.
-//
-// MFEM is free software; you can redistribute it and/or modify it under the
-// terms of the BSD-3 license. We welcome feedback and contributions, see file
-// CONTRIBUTING.md for details.
-
-
-// Abstract array data type
-
-#include "array.hpp"
-#include "../general/forall.hpp"
-#include <fstream>
-#include <type_traits>
-
-namespace mfem
-{
-
-template <class T>
-void Array<T>::Print(std::ostream &os, int width) const
-{
-   for (int i = 0; i < size; i++)
-   {
-      os << data[i];
-      if ( !((i+1) % width) || i+1 == size )
-      {
-         os << '\n';
-      }
-      else
-      {
-         os << " ";
-      }
-   }
-}
-
-template <class T>
-void Array<T>::Save(std::ostream &os, int fmt) const
-{
-   if (fmt == 0)
-   {
-      os << size << '\n';
-   }
-   for (int i = 0; i < size; i++)
-   {
-      os << operator[](i) << '\n';
-   }
-}
-
-template <class T>
-void Array<T>::Load(std::istream &in, int fmt)
-{
-   if (fmt == 0)
-   {
-      int new_size;
-      in >> new_size;
-      SetSize(new_size);
-   }
-   for (int i = 0; i < size; i++)
-   {
-      in >> operator[](i);
-   }
-}
-
-template <class T>
-T Array<T>::Max() const
-{
-   MFEM_ASSERT(size > 0, "Array is empty with size " << size);
-
-   T max = operator[](0);
-   for (int i = 1; i < size; i++)
-   {
-      if (max < operator[](i))
-      {
-         max = operator[](i);
-      }
-   }
-
-   return max;
-}
-
-template <class T>
-T Array<T>::Min() const
-{
-   MFEM_ASSERT(size > 0, "Array is empty with size " << size);
-
-   T min = operator[](0);
-   for (int i = 1; i < size; i++)
-   {
-      if (operator[](i) < min)
-      {
-         min = operator[](i);
-      }
-   }
-
-   return min;
-}
-
-// Partial Sum
-template <class T>
-void Array<T>::PartialSum()
-{
-   T sum = static_cast<T>(0);
-   for (int i = 0; i < size; i++)
-   {
-      sum+=operator[](i);
-      operator[](i) = sum;
-   }
-}
-
-template <class T>
-void Array<T>::Abs()
-{
-   static_assert(std::is_arithmetic<T>::value, "Use with arithmetic types!");
-   const bool useDevice = UseDevice();
-   const int N = size;
-   auto y = ReadWrite(useDevice);
-   mfem::forall_switch(useDevice, N, [=] MFEM_HOST_DEVICE (int i)
-   {
-      y[i] = std::abs(y[i]);
-   });
-}
-
-// Sum
-template <class T>
-T Array<T>::Sum() const
-{
-   T sum = static_cast<T>(0);
-   for (int i = 0; i < size; i++)
-   {
-      sum+=operator[](i);
-   }
-
-   return sum;
-}
-
-template <class T>
-int Array<T>::IsSorted() const
-{
-   T val_prev = operator[](0), val;
-   for (int i = 1; i < size; i++)
-   {
-      val=operator[](i);
-      if (val < val_prev)
-      {
-         return 0;
-      }
-      val_prev = val;
-   }
-
-   return 1;
-}
-
-template <class T>
-bool Array<T>::IsConstant() const
-{
-   if (size < 2) { return true; }
-   const T v0 = data[0];
-   for (int i = 1; i < size; i++)
-   {
-      if (data[i] != v0)
-      {
-         return false;
-      }
-   }
-
-   return true;
-}
-
-template <class T>
-void Array2D<T>::Load(const char *filename, int fmt)
-{
-   std::ifstream in;
-   in.open(filename, std::ifstream::in);
-   MFEM_VERIFY(in.is_open(), "File " << filename << " does not exist.");
-   Load(in, fmt);
-   in.close();
-}
-
-template <class T>
-void Array2D<T>::Print(std::ostream &os, int width_)
-{
-   int height = this->NumRows();
-   int width  = this->NumCols();
-
-   for (int i = 0; i < height; i++)
-   {
-      os << "[row " << i << "]\n";
-      for (int j = 0; j < width; j++)
-      {
-         os << (*this)(i,j);
-         if ( (j+1) == width_ || (j+1) % width_ == 0 )
-         {
-            os << '\n';
-         }
-         else
-         {
-            os << ' ';
-         }
-      }
-   }
-}
-
-template class Array<char>;
-template class Array<int>;
-template class Array<long long>;
-template class Array<real_t>;
-template class Array2D<int>;
-template class Array2D<real_t>;
-
-} // namespace mfem

general/array.hpp

@@ -16,9 +16,13 @@
 #include "mem_manager.hpp"
 #include "device.hpp"
 #include "error.hpp"
+#include "forall.hpp"
 #include "globals.hpp"
+#include "reducers.hpp"
+#include "scan.hpp"
 
 #include <iostream>
+#include <fstream>
 #include <cstdlib>
 #include <cstring>
 #include <algorithm>
@@ -135,6 +139,8 @@ public:
    /// Return the device flag of the Memory object used by the Array
    bool UseDevice() const { return data.UseDevice(); }
 
+   void UseDevice(bool use_dev) { data.UseDevice(use_dev); }
+
    /// Return true if the data will be deleted by the Array
    inline bool OwnsData() const { return data.OwnsHostPtr(); }
 
@@ -275,11 +281,11 @@ public:
 
    /** @brief Find the maximal element in the array, using the comparison
        operator `<` for class T. */
-   T Max() const;
+   inline T Max() const;
 
    /** @brief Find the minimal element in the array, using the comparison
        operator `<` for class T. */
-   T Min() const;
+   inline T Min() const;
 
    /// Sorts the array in ascending order. This requires operator< to be defined for T.
    void Sort() { std::sort((T*)data, data + size); }
@@ -297,22 +303,22 @@ public:
    }
 
    /// Return 1 if the array is sorted from lowest to highest.  Otherwise return 0.
-   int IsSorted() const;
+   inline int IsSorted() const;
 
    /// Does the Array have Size zero.
    bool IsEmpty() const { return Size() == 0; }
 
    /// Return true if all entries of the array are the same.
-   bool IsConstant() const;
+   inline bool IsConstant() const;
 
    /// Fill the entries of the array with the cumulative sum of the entries.
-   void PartialSum();
+   inline void PartialSum();
 
    /// Replace each entry of the array with its absolute value.
-   void Abs();
+   inline void Abs();
 
    /// Return the sum of all the array entries using the '+'' operator for class 'T'.
-   T Sum() const;
+   inline T Sum() const;
 
    /// Set all entries of the array to the provided constant.
    inline void operator=(const T &a);
@@ -797,8 +803,14 @@ template <typename T> template <typename CT>
 inline Array<T> &Array<T>::operator=(const Array<CT> &src)
 {
    SetSize(src.Size());
-   for (int i = 0; i < size; i++) { (*this)[i] = T(src[i]); }
-   return *this;
+
+   const bool use_dev = UseDevice() || src.UseDevice();
+   const auto x = src.Read(use_dev);
+   auto y = Write(use_dev);
+   mfem::forall_switch(use_dev, size, [=] MFEM_HOST_DEVICE (int i)
+   {
+      y[i] = x[i];
+   });
 }
 
 template <class T>
@@ -1014,19 +1026,24 @@ template <class T>
 inline void Array<T>::GetSubArray(int offset, int sa_size, Array<T> &sa) const
 {
    sa.SetSize(sa_size);
-   for (int i = 0; i < sa_size; i++)
+   const bool use_dev = UseDevice() || sa.UseDevice();
+   const auto x = Read(use_dev);
+   auto y = sa.Write(use_dev);
+   mfem::forall_switch(use_dev, sa_size, [=] MFEM_HOST_DEVICE (int i)
    {
-      sa[i] = (*this)[offset+i];
-   }
+      y[i] = x[offset + i];
+   });
 }
 
 template <class T>
 inline void Array<T>::operator=(const T &a)
 {
-   for (int i = 0; i < size; i++)
+   const bool use_dev = UseDevice();
+   auto x = Write(use_dev);
+   mfem::forall_switch(use_dev, size, [=] MFEM_HOST_DEVICE (int i)
    {
-      data[i] = a;
-   }
+      x[i] = a;
+   });
 }
 
 template <class T>
@@ -1035,6 +1052,153 @@ inline void Array<T>::Assign(const T *p)
    data.CopyFromHost(p, Size());
 }
 
+template <class T>
+inline void Array<T>::Print(std::ostream &os, int width) const
+{
+   for (int i = 0; i < size; i++)
+   {
+      os << data[i];
+      if ( !((i+1) % width) || i+1 == size )
+      {
+         os << '\n';
+      }
+      else
+      {
+         os << " ";
+      }
+   }
+}
+
+template <class T>
+inline void Array<T>::Save(std::ostream &os, int fmt) const
+{
+   if (fmt == 0)
+   {
+      os << size << '\n';
+   }
+   for (int i = 0; i < size; i++)
+   {
+      os << operator[](i) << '\n';
+   }
+}
+
+template <class T>
+void Array<T>::Load(std::istream &in, int fmt)
+{
+   if (fmt == 0)
+   {
+      int new_size;
+      in >> new_size;
+      SetSize(new_size);
+   }
+   for (int i = 0; i < size; i++)
+   {
+      in >> operator[](i);
+   }
+}
+
+template <class T>
+inline T Array<T>::Max() const
+{
+   MFEM_ASSERT(size > 0, "Array is empty with size " << size);
+
+   T max = operator[](0);
+   for (int i = 1; i < size; i++)
+   {
+      if (max < operator[](i))
+      {
+         max = operator[](i);
+      }
+   }
+
+   return max;
+}
+
+template <class T>
+inline T Array<T>::Min() const
+{
+   MFEM_ASSERT(size > 0, "Array is empty with size " << size);
+
+   T min = operator[](0);
+   for (int i = 1; i < size; i++)
+   {
+      if (operator[](i) < min)
+      {
+         min = operator[](i);
+      }
+   }
+
+   return min;
+}
+
+// Partial Sum
+template <class T>
+inline void Array<T>::PartialSum()
+{
+   auto data_ptr = ReadWrite(UseDevice());
+   InclusiveScan(UseDevice(), data_ptr, data_ptr, size);
+}
+
+template <class T>
+inline void Array<T>::Abs()
+{
+   static_assert(std::is_arithmetic<T>::value, "Use with arithmetic types!");
+   const bool useDevice = UseDevice();
+   const int N = size;
+   auto y = ReadWrite(useDevice);
+   mfem::forall_switch(useDevice, N, [=] MFEM_HOST_DEVICE (int i)
+   {
+      y[i] = std::abs(y[i]);
+   });
+}
+
+// Sum
+template <class T>
+inline T Array<T>::Sum() const
+{
+   T sum = static_cast<T>(0);
+   if (size > 0)
+   {
+      const auto m_data = Read(UseDevice());
+      reduce(size, sum, [=] MFEM_HOST_DEVICE(int i, T &r) { r += m_data[i]; },
+      /*  */ SumReducer<T> {}, UseDevice());
+   }
+   return sum;
+}
+
+template <class T>
+inline int Array<T>::IsSorted() const
+{
+   T val_prev = operator[](0), val;
+   for (int i = 1; i < size; i++)
+   {
+      val=operator[](i);
+      if (val < val_prev)
+      {
+         return 0;
+      }
+      val_prev = val;
+   }
+
+   return 1;
+}
+
+template <class T>
+inline bool Array<T>::IsConstant() const
+{
+   if (size < 2) { return true; }
+   const T v0 = data[0];
+   for (int i = 1; i < size; i++)
+   {
+      if (data[i] != v0)
+      {
+         return false;
+      }
+   }
+
+   return true;
+}
+
 
 template <class T>
 inline const T &Array2D<T>::operator()(int i, int j) const
@@ -1074,6 +1238,40 @@ inline T *Array2D<T>::operator[](int i)
    return &array1d[i*N];
 }
 
+template <class T>
+void Array2D<T>::Load(const char *filename, int fmt)
+{
+   std::ifstream in;
+   in.open(filename, std::ifstream::in);
+   MFEM_VERIFY(in.is_open(), "File " << filename << " does not exist.");
+   Load(in, fmt);
+   in.close();
+}
+
+template <class T>
+void Array2D<T>::Print(std::ostream &os, int width_)
+{
+   int height = this->NumRows();
+   int width  = this->NumCols();
+
+   for (int i = 0; i < height; i++)
+   {
+      os << "[row " << i << "]\n";
+      for (int j = 0; j < width; j++)
+      {
+         os << (*this)(i,j);
+         if ( (j+1) == width_ || (j+1) % width_ == 0 )
+         {
+            os << '\n';
+         }
+         else
+         {
+            os << ' ';
+         }
+      }
+   }
+}
+
 
 template <class T>
 inline void Swap(Array2D<T> &a, Array2D<T> &b)

general/reducers.hpp

@@ -12,7 +12,6 @@
 #ifndef MFEM_REDUCERS_HPP
 #define MFEM_REDUCERS_HPP
 
-#include "array.hpp"
 #include "forall.hpp"
 
 #include <cmath>
@@ -514,6 +513,33 @@ template<class B, class R> struct reduction_kernel
       }
    }
 };
+
+template <class T>
+class ReductionWorkspace
+{
+   Memory<T> workspace;
+
+   static ReductionWorkspace &Instance()
+   {
+      static ReductionWorkspace instance;
+      return instance;
+   }
+
+   ~ReductionWorkspace() { workspace.Delete(); }
+
+public:
+   static T *Get(int num_blocks)
+   {
+      ReductionWorkspace &instance = Instance();
+      if (instance.workspace.Capacity() < num_blocks)
+      {
+         instance.workspace.Delete();
+         instance.workspace.New(num_blocks, MemoryType::HOST_PINNED);
+      }
+      return instance.workspace;
+   }
+};
+
 }
 
 /**
@@ -529,8 +555,7 @@ template<class B, class R> struct reduction_kernel
  @tparam T value_type to operate on
  */
 template <class T, class B, class R>
-void reduce(int N, T &res, B &&body, const R &reducer, bool use_dev,
-            Array<T> &workspace)
+void reduce(int N, T &res, B &&body, const R &reducer, bool use_dev)
 {
    if (N == 0)
    {
@@ -567,13 +592,7 @@ void reduce(int N, T &res, B &&body, const R &reducer, bool use_dev,
 
       red_type red{nullptr, std::forward<B>(body), reducer, N, items_per_thread};
       // allocate res to fit block_size entries
-      auto mt = workspace.GetMemory().GetMemoryType();
-      if (mt != MemoryType::HOST_PINNED && mt != MemoryType::MANAGED)
-      {
-         mt = MemoryType::HOST_PINNED;
-      }
-      workspace.SetSize(nblocks, mt);
-      auto work = workspace.HostWrite();
+      auto work = internal::ReductionWorkspace<T>::Get(nblocks);
       red.work = work;
       forall_2D(nblocks, block_size, 1, std::move(red));
       // wait for results

general/scan.hpp

@@ -28,8 +28,37 @@
 
 namespace mfem
 {
-/// Equivalent to InclusiveScan(use_dev, d_in, d_out, num_items, workspace,
-/// std::plus<>{})
+
+namespace internal
+{
+class ScanWorkspace
+{
+   Memory<std::byte> workspace;
+   static ScanWorkspace &Instance()
+   {
+      static ScanWorkspace instance;
+      return instance;
+   }
+   ~ScanWorkspace() { workspace.Delete(); }
+public:
+   static std::byte *Get(int num_bytes)
+   {
+      ScanWorkspace &instance = Instance();
+      if (Size() < num_bytes)
+      {
+         instance.workspace.Delete();
+         instance.workspace.New(num_bytes);
+      }
+      return instance.workspace.Write(MemoryClass::DEVICE, Size());
+   }
+   static int Size()
+   {
+      return Instance().workspace.Capacity();
+   }
+};
+}
+
+/// Equivalent to InclusiveScan(use_dev, d_in, d_out, num_items, std::plus<>{})
 template <class InputIt, class OutputIt>
 void InclusiveScan(bool use_dev, InputIt d_in, OutputIt d_out, size_t num_items)
 {
@@ -37,12 +66,12 @@ void InclusiveScan(bool use_dev, InputIt d_in, OutputIt d_out, size_t num_items)
 #if defined(MFEM_USE_CUDA) || defined(MFEM_USE_HIP)
    if (use_dev && mfem::Device::Allows(Backend::CUDA_MASK | Backend::HIP_MASK))
    {
-      static Array<std::byte> workspace;
-      size_t bytes = workspace.Size();
-      if (bytes)
+      using internal::ScanWorkspace;
+      size_t bytes = ScanWorkspace::Size();
+      if (bytes > 0)
       {
          auto err = MFEM_CUB_NAMESPACE::DeviceScan::InclusiveSum(
-                       workspace.Write(), bytes, d_in, d_out, num_items);
+                       ScanWorkspace::Get(bytes), bytes, d_in, d_out, num_items);
 #if defined(MFEM_USE_CUDA)
          if (err == cudaSuccess)
          {
@@ -57,11 +86,12 @@ void InclusiveScan(bool use_dev, InputIt d_in, OutputIt d_out, size_t num_items)
       }
       // try allocating a larger buffer
       bytes = 0;
+      // get size of buffer
       MFEM_GPU_CHECK(MFEM_CUB_NAMESPACE::DeviceScan::InclusiveSum(
                         nullptr, bytes, d_in, d_out, num_items));
-      workspace.SetSize(bytes);
+      // resize buffer (in ScanWorkspace::Get) and try again
       MFEM_GPU_CHECK(MFEM_CUB_NAMESPACE::DeviceScan::InclusiveSum(
-                        workspace.Write(), bytes, d_in, d_out, num_items));
+                        ScanWorkspace::Get(bytes), bytes, d_in, d_out, num_items));
       return;
    }
 #endif
@@ -101,12 +131,13 @@ void InclusiveScan(bool use_dev, InputIt d_in, OutputIt d_out, size_t num_items,
 #if defined(MFEM_USE_CUDA) || defined(MFEM_USE_HIP)
    if (use_dev && mfem::Device::Allows(Backend::CUDA_MASK | Backend::HIP_MASK))
    {
-      static Array<std::byte> workspace;
-      size_t bytes = workspace.Size();
-      if (bytes)
+      using internal::ScanWorkspace;
+      size_t bytes = ScanWorkspace::Size();
+      if (bytes > 0)
       {
          auto err = MFEM_CUB_NAMESPACE::DeviceScan::InclusiveScan(
-                       workspace.Write(), bytes, d_in, d_out, scan_op, num_items);
+                       ScanWorkspace::Get(bytes), bytes, d_in, d_out, scan_op,
+                       num_items);
 #if defined(MFEM_USE_CUDA)
          if (err == cudaSuccess)
          {
@@ -123,9 +154,9 @@ void InclusiveScan(bool use_dev, InputIt d_in, OutputIt d_out, size_t num_items,
       bytes = 0;
       MFEM_GPU_CHECK(MFEM_CUB_NAMESPACE::DeviceScan::InclusiveScan(
                         nullptr, bytes, d_in, d_out, scan_op, num_items));
-      workspace.SetSize(bytes);
       MFEM_GPU_CHECK(MFEM_CUB_NAMESPACE::DeviceScan::InclusiveScan(
-                        workspace.Write(), bytes, d_in, d_out, scan_op, num_items));
+                        ScanWorkspace::Get(bytes), bytes, d_in, d_out, scan_op,
+                        num_items));
       return;
    }
 #endif
@@ -164,13 +195,13 @@ void ExclusiveScan(bool use_dev, InputIt d_in, OutputIt d_out, size_t num_items,
 #if defined(MFEM_USE_CUDA) || defined(MFEM_USE_HIP)
    if (use_dev && mfem::Device::Allows(Backend::CUDA_MASK | Backend::HIP_MASK))
    {
-      static Array<std::byte> workspace;
-      size_t bytes = workspace.Size();
+      using internal::ScanWorkspace;
+      size_t bytes = ScanWorkspace::Size();
       if (bytes)
       {
          auto err = MFEM_CUB_NAMESPACE::DeviceScan::ExclusiveScan(
-                       workspace.Write(), bytes, d_in, d_out, scan_op, init_value,
-                       num_items);
+                       ScanWorkspace::Get(bytes), bytes, d_in, d_out, scan_op,
+                       init_value, num_items);
 #if defined(MFEM_USE_CUDA)
          if (err == cudaSuccess)
          {
@@ -187,10 +218,9 @@ void ExclusiveScan(bool use_dev, InputIt d_in, OutputIt d_out, size_t num_items,
       bytes = 0;
       MFEM_GPU_CHECK(MFEM_CUB_NAMESPACE::DeviceScan::ExclusiveScan(
                         nullptr, bytes, d_in, d_out, scan_op, init_value, num_items));
-      workspace.SetSize(bytes);
       MFEM_GPU_CHECK(MFEM_CUB_NAMESPACE::DeviceScan::ExclusiveScan(
-                        workspace.Write(), bytes, d_in, d_out, scan_op, init_value,
-                        num_items));
+                        ScanWorkspace::Get(bytes), bytes, d_in, d_out, scan_op,
+                        init_value, num_items));
       return;
    }
 #endif
@@ -213,7 +243,7 @@ void ExclusiveScan(bool use_dev, InputIt d_in, OutputIt d_out, size_t num_items,
 }
 
 /// Equivalent to ExclusiveScan(use_dev, d_in, d_out, num_items, init_value,
-/// workspace, std::plus<>{})
+/// std::plus<>{})
 template <class InputIt, class OutputIt, class T>
 void ExclusiveScan(bool use_dev, InputIt d_in, OutputIt d_out, size_t num_items,
                    T init_value)

linalg/vector.cpp

@@ -92,18 +92,6 @@ struct LpReducer
    }
 };
 
-static Array<real_t>& vector_workspace()
-{
-   static Array<real_t> instance;
-   return instance;
-}
-
-static Array<DevicePair<real_t, real_t>> &Lpvector_workspace()
-{
-   static Array<DevicePair<real_t, real_t>> instance;
-   return instance;
-}
-
 Vector::Vector(const Vector &v)
 {
    const int s = v.Size();
@@ -991,7 +979,7 @@ real_t Vector::Norml2() const
          }
       }
    },
-   L2Reducer{}, UseDevice(), Lpvector_workspace());
+   L2Reducer{}, UseDevice());
    // final answer
    return res.second * sqrt(res.first);
 }
@@ -1006,7 +994,7 @@ real_t Vector::Normlinf() const
    {
       r = fmax(r, fabs(m_data[i]));
    },
-   MaxReducer<real_t> {}, UseDevice(), vector_workspace());
+   MaxReducer<real_t> {}, UseDevice());
    return res;
 }
 
@@ -1020,7 +1008,7 @@ real_t Vector::Norml1() const
    {
       r += fabs(m_data[i]);
    },
-   SumReducer<real_t> {}, UseDevice(), vector_workspace());
+   SumReducer<real_t> {}, UseDevice());
    return res;
 }
 
@@ -1063,7 +1051,7 @@ real_t Vector::Normlp(real_t p) const
             }
          }
       },
-      LpReducer{p}, UseDevice(), Lpvector_workspace());
+      LpReducer{p}, UseDevice());
       // final answer
       return res.second * pow(res.first, 1.0 / p);
    } // end if p < infinity()
@@ -1096,7 +1084,7 @@ real_t Vector::operator*(const Vector &v) const
       {
          r += m_data[i] * v_data[i];
       },
-      SumReducer<real_t> {}, use_dev, vector_workspace());
+      SumReducer<real_t> {}, use_dev);
       return res;
    };
 
@@ -1167,7 +1155,7 @@ real_t Vector::Min() const
       {
          r = fmin(r, m_data[i]);
       },
-      MinReducer<real_t> {}, use_dev, vector_workspace());
+      MinReducer<real_t> {}, use_dev);
       return res;
    };
 
@@ -1213,7 +1201,7 @@ real_t Vector::Max() const
       {
          r = fmax(r, m_data[i]);
       },
-      MaxReducer<real_t> {}, use_dev, vector_workspace());
+      MaxReducer<real_t> {}, use_dev);
       return res;
    };
 
@@ -1248,7 +1236,7 @@ real_t Vector::Sum() const
    {
       r += m_data[i];
    },
-   SumReducer<real_t> {}, UseDevice(), vector_workspace());
+   SumReducer<real_t> {}, UseDevice());
    return res;
 }
 

linalg/vector_operator.cpp

@@ -1,80 +0,0 @@
-// Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
-// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
-// LICENSE and NOTICE for details. LLNL-CODE-806117.
-//
-// This file is part of the MFEM library. For more information and source code
-// availability visit https://mfem.org.
-//
-// MFEM is free software; you can redistribute it and/or modify it under the
-// terms of the BSD-3 license. We welcome feedback and contributions, see file
-// CONTRIBUTING.md for details.
-
-#include "vector_operator.hpp"
-
-namespace mfem
-{
-
-#ifdef MFEM_USE_MPI
-
-ParVectorOperator::ParVectorOperator(MPI_Comm comm,
-                                     int myid,
-                                     int local_vec_size,
-                                     int num_vecs)
-   : Operator((myid == 0) ? num_vecs : 0, local_vec_size),
-     comm(comm),
-     myid(myid),
-     vecs(num_vecs),
-     coefs(num_vecs),
-     owns(num_vecs)
-{
-   vecs = NULL;
-   coefs = 1.0;
-   owns = false;
-}
-
-ParVectorOperator::~ParVectorOperator()
-{
-   for (int i=0; i < vecs.Size(); i++)
-   {
-      if (owns[i]) { delete vecs[i]; }
-      vecs[i] = NULL;
-   }
-}
-
-void ParVectorOperator::SetVector(int idx, Vector *vec,
-                                  double c, bool own_vec)
-{
-   MFEM_VERIFY(idx >= 0 && idx < vecs.Size(),
-               "ParVectorOperator: Index out of range");
-
-   vecs[idx] = vec;
-   coefs[idx] = c;
-   owns[idx] = own_vec;
-}
-
-void ParVectorOperator::Mult(const Vector &x, Vector &y) const
-{
-   for (int i=0; i<vecs.Size(); i++)
-   {
-      double vo = coefs[i] * (*vecs[i] * x);
-      double vi = 0.0;
-      MPI_Reduce(&vo, &vi, 1, MPI_DOUBLE, MPI_SUM, 0, comm);
-      if (myid == 0) { y[i] = vi; }
-   }
-}
-
-/// Action of the transpose operator: `y=A^t(x)`.
-void ParVectorOperator::MultTranspose(const Vector &x, Vector &y) const
-{
-   y = 0.0;
-   for (int i=0; i<vecs.Size(); i++)
-   {
-      double xi = (myid == 0) ? x[i] : 0.0;
-      MPI_Bcast(&xi, 1, MPI_DOUBLE, 0, comm);
-      y.Add(xi * coefs[i], *vecs[i]);
-   }
-}
-
-}
-
-#endif // MFEM_USE_MPI

linalg/vector_operator.hpp

@@ -1,55 +0,0 @@
-// Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
-// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
-// LICENSE and NOTICE for details. LLNL-CODE-806117.
-//
-// This file is part of the MFEM library. For more information and source code
-// availability visit https://mfem.org.
-//
-// MFEM is free software; you can redistribute it and/or modify it under the
-// terms of the BSD-3 license. We welcome feedback and contributions, see file
-// CONTRIBUTING.md for details.
-
-#ifndef MFEM_VECTOR_OPERATOR
-#define MFEM_VECTOR_OPERATOR
-
-#include "operator.hpp"
-#include "vector.hpp"
-
-namespace mfem
-{
-
-#ifdef MFEM_USE_MPI
-
-class ParVectorOperator : public Operator
-{
-private:
-   MPI_Comm comm;
-   int myid;
-
-   Array<Vector*> vecs;
-   Array<double> coefs;
-   Array<bool> owns;
-
-public:
-   ParVectorOperator(MPI_Comm comm,
-                     int myid,
-                     int local_vec_size,
-                     int num_vecs);
-
-   ~ParVectorOperator();
-
-   void SetVector(int idx, Vector *vec,
-                  double c = 1.0, bool own_vec = false);
-
-   /// Operator application: `y=A(x)`.
-   void Mult(const Vector &x, Vector &y) const;
-
-   /// Action of the transpose operator: `y=A^t(x)`.
-   void MultTranspose(const Vector &x, Vector &y) const;
-};
-
-#endif // MFEM_USE_MPI
-
-} // namespace mfem
-
-#endif // MFEM_VECTOR_OPERATOR

tests/unit/general/test_reduction.cpp

@@ -22,7 +22,6 @@ using namespace mfem;
 
 TEST_CASE("Reduce Sum", "[Reduction],[GPU]")
 {
-   Array<int> workspace;
    Array<int> a(1000);
    a.HostReadWrite();
    for (int i = 0; i < a.Size(); ++i)
@@ -36,7 +35,7 @@ TEST_CASE("Reduce Sum", "[Reduction],[GPU]")
       int res = 0;
       mfem::reduce(
       a.Size(), res, [=] MFEM_HOST_DEVICE(int i, int &r) { r += dptr[i]; },
-      SumReducer<int> {}, use_dev, workspace);
+      SumReducer<int> {}, use_dev);
       // correct for even-length summations
       int expected = (AsConst(a)[0] + AsConst(a)[a.Size() - 1]) * a.Size() / 2;
       CAPTURE(use_dev);
@@ -46,7 +45,6 @@ TEST_CASE("Reduce Sum", "[Reduction],[GPU]")
 
 TEST_CASE("Reduce Mult", "[Reduction],[GPU]")
 {
-   Array<long long> workspace;
    Array<long long> a(64);
    a.HostReadWrite();
    for (int i = 0; i < a.Size(); ++i)
@@ -64,7 +62,7 @@ TEST_CASE("Reduce Mult", "[Reduction],[GPU]")
          mfem::reduce(
             a.Size(), res,
          [=] MFEM_HOST_DEVICE(int i, long long &r) { r *= dptr[i]; },
-         MultReducer<long long> {}, use_dev, workspace);
+         MultReducer<long long> {}, use_dev);
          long long expected = 0;
          CAPTURE(use_dev);
          REQUIRE(res == expected);
@@ -76,7 +74,7 @@ TEST_CASE("Reduce Mult", "[Reduction],[GPU]")
          mfem::reduce(
             a.Size(), res,
          [=] MFEM_HOST_DEVICE(int i, long long &r) { r *= dptr[i]; },
-         MultReducer<long long> {}, use_dev, workspace);
+         MultReducer<long long> {}, use_dev);
          long long expected = 21936950640377856;
          CAPTURE(use_dev);
          REQUIRE(res == expected);
@@ -86,7 +84,6 @@ TEST_CASE("Reduce Mult", "[Reduction],[GPU]")
 
 TEST_CASE("Reduce BAnd", "[Reduction],[GPU]")
 {
-   Array<unsigned> workspace;
    Array<unsigned> a(10);
    SECTION("{ Bit unset }")
    {
@@ -108,7 +105,7 @@ TEST_CASE("Reduce BAnd", "[Reduction],[GPU]")
          mfem::reduce(
             a.Size(), res,
          [=] MFEM_HOST_DEVICE(int i, unsigned &r) { r &= dptr[i]; },
-         BAndReducer<unsigned> {}, use_dev, workspace);
+         BAndReducer<unsigned> {}, use_dev);
          CAPTURE(use_dev);
          REQUIRE(res == ((~1u) & ~(1u << unset_bit)));
          REQUIRE((res & (1u << unset_bit)) == 0);
@@ -132,7 +129,7 @@ TEST_CASE("Reduce BAnd", "[Reduction],[GPU]")
          mfem::reduce(
             a.Size(), res,
          [=] MFEM_HOST_DEVICE(int i, unsigned &r) { r &= dptr[i]; },
-         BAndReducer<unsigned> {}, use_dev, workspace);
+         BAndReducer<unsigned> {}, use_dev);
          CAPTURE(use_dev);
          REQUIRE(res == (1u << set_bit));
       }
@@ -141,7 +138,6 @@ TEST_CASE("Reduce BAnd", "[Reduction],[GPU]")
 
 TEST_CASE("Reduce BOr", "[Reduction],[GPU]")
 {
-   Array<unsigned> workspace;
    Array<unsigned> a(0x210);
    a.HostReadWrite();
    for (int i = 0; i < a.Size(); ++i)
@@ -157,7 +153,7 @@ TEST_CASE("Reduce BOr", "[Reduction],[GPU]")
       mfem::reduce(
          a.Size(), res,
       [=] MFEM_HOST_DEVICE(int i, unsigned &r) { r |= dptr[i]; },
-      BOrReducer<unsigned> {}, use_dev, workspace);
+      BOrReducer<unsigned> {}, use_dev);
       CAPTURE(use_dev);
       REQUIRE(res == 0x3ffu);
    }
@@ -165,7 +161,6 @@ TEST_CASE("Reduce BOr", "[Reduction],[GPU]")
 
 TEST_CASE("Reduce Min", "[Reduction],[GPU]")
 {
-   Array<int> workspace;
    Array<int> a(1000);
    auto hptr = a.HostReadWrite();
    for (int i = 0; i < a.Size(); ++i)
@@ -190,7 +185,7 @@ TEST_CASE("Reduce Min", "[Reduction],[GPU]")
             r = dptr[i];
          }
       },
-      MinReducer<int> {}, use_dev, workspace);
+      MinReducer<int> {}, use_dev);
       CAPTURE(use_dev);
       REQUIRE(res == -10);
    }
@@ -198,7 +193,6 @@ TEST_CASE("Reduce Min", "[Reduction],[GPU]")
 
 TEST_CASE("Reduce Max", "[Reduction],[GPU]")
 {
-   Array<int> workspace;
    Array<int> a(1000);
    auto hptr = a.HostReadWrite();
    for (int i = 0; i < a.Size(); ++i)
@@ -223,7 +217,7 @@ TEST_CASE("Reduce Max", "[Reduction],[GPU]")
             r = dptr[i];
          }
       },
-      MaxReducer<int> {}, use_dev, workspace);
+      MaxReducer<int> {}, use_dev);
       CAPTURE(use_dev);
       REQUIRE(res == 999 - 10);
    }
@@ -231,7 +225,6 @@ TEST_CASE("Reduce Max", "[Reduction],[GPU]")
 
 TEST_CASE("Reduce MinMax", "[Reduction],[GPU]")
 {
-   Array<DevicePair<int, int>> workspace;
    Array<int> a(1000);
    auto hptr = a.HostReadWrite();
    for (int i = 0; i < a.Size(); ++i)
@@ -262,7 +255,7 @@ TEST_CASE("Reduce MinMax", "[Reduction],[GPU]")
             r.second = dptr[i];
          }
       },
-      MinMaxReducer<int> {}, use_dev, workspace);
+      MinMaxReducer<int> {}, use_dev);
       CAPTURE(use_dev);
       REQUIRE(res.first == -10);
       REQUIRE(res.second == a.Size() - 11);
@@ -271,7 +264,6 @@ TEST_CASE("Reduce MinMax", "[Reduction],[GPU]")
 
 TEST_CASE("Reduce ArgMin", "[Reduction],[GPU]")
 {
-   Array<DevicePair<double, int>> workspace;
    Array<double> a(1000);
    auto hptr = a.HostReadWrite();
    for (int i = 0; i < a.Size(); ++i)
@@ -297,7 +289,7 @@ TEST_CASE("Reduce ArgMin", "[Reduction],[GPU]")
             r.second = i;
          }
       },
-      ArgMinReducer<double, int> {}, use_dev, workspace);
+      ArgMinReducer<double, int> {}, use_dev);
       CAPTURE(use_dev);
       REQUIRE(res.first == -10);
       REQUIRE(res.second >= 0);
@@ -308,7 +300,6 @@ TEST_CASE("Reduce ArgMin", "[Reduction],[GPU]")
 
 TEST_CASE("Reduce ArgMax", "[Reduction],[GPU]")
 {
-   Array<DevicePair<double, int>> workspace;
    Array<double> a(1000);
 
    auto hptr = a.HostReadWrite();
@@ -337,7 +328,7 @@ TEST_CASE("Reduce ArgMax", "[Reduction],[GPU]")
             r.second = i;
          }
       },
-      ArgMaxReducer<double, int> {}, use_dev, workspace);
+      ArgMaxReducer<double, int> {}, use_dev);
       CAPTURE(use_dev);
       REQUIRE(res.first == a.Size() - 11);
       REQUIRE(res.second >= 0);
@@ -348,7 +339,6 @@ TEST_CASE("Reduce ArgMax", "[Reduction],[GPU]")
 
 TEST_CASE("Reduce ArgMinMax", "[Reduction],[GPU]")
 {
-   Array<MinMaxLocScalar<double, int>> workspace;
    Array<double> a(1000);
    auto hptr = a.HostReadWrite();
    for (int i = 0; i < a.Size(); ++i)
@@ -383,7 +373,7 @@ TEST_CASE("Reduce ArgMinMax", "[Reduction],[GPU]")
             r.max_loc = i;
          }
       },
-      ArgMinMaxReducer<double, int> {}, use_dev, workspace);
+      ArgMinMaxReducer<double, int> {}, use_dev);
       CAPTURE(use_dev);
       REQUIRE(res.min_val == -10);
       REQUIRE(res.min_loc >= 0);