Improve #include logic.

default stream 0.

Optimize OperatorJacobiSmoother::Mult by fusing kernels as much
as possible.

general/backends.hpp

@@ -31,6 +31,7 @@
 #endif
 
 #ifdef MFEM_USE_RAJA
+#define CAMP_USE_PLATFORM_DEFAULT_STREAM 1
 #include "RAJA/RAJA.hpp"
 #if defined(RAJA_ENABLE_CUDA) && !defined(MFEM_USE_CUDA)
 #error When RAJA is built with CUDA, MFEM_USE_CUDA=YES is required

general/device.cpp

@@ -11,6 +11,9 @@
 
 #include "forall.hpp"
 #include "occa.hpp"
+#if defined(MFEM_USE_OCCA) && OCCA_CUDA_ENABLED
+#include <occa/modes/cuda/stream.hpp>
+#endif
 #ifdef MFEM_USE_CEED
 #include "../fem/libceed/ceed.hpp"
 #endif
@@ -410,6 +413,10 @@ static void OccaDeviceSetup(const int dev)
 #if OCCA_CUDA_ENABLED
       std::string mode("mode: 'CUDA', device_id : ");
       internal::occaDevice.setup(mode.append(1,'0'+dev));
+      auto def_stream = new occa::cuda::stream(
+         internal::occaDevice.getModeDevice(),
+         occa::properties(), (CUstream)0);
+      internal::occaDevice.setStream(def_stream);
 #else
       MFEM_ABORT("the OCCA CUDA backend requires OCCA built with CUDA!");
 #endif

linalg/solvers.cpp

@@ -162,21 +162,22 @@ void OperatorJacobiSmoother::Mult(const Vector &x, Vector &y) const
    MFEM_ASSERT(x.Size() == N, "invalid input vector");
    MFEM_ASSERT(y.Size() == N, "invalid output vector");
 
+   auto DI = dinv.Read();
+   auto X = x.Read();
    if (iterative_mode && oper)
    {
-      oper->Mult(y, residual);  // r = A x
-      subtract(x, residual, residual); // r = b - A x
+      oper->Mult(y, residual);  // r = A y
+      auto R = residual.Read();
+      auto Y = y.ReadWrite();
+      // y += D^{-1} (x - A y)
+      MFEM_FORALL(i, N, Y[i] += DI[i] * (X[i] - R[i]); );
    }
    else
    {
-      residual = x;
-      y.UseDevice(true);
-      y = 0.0;
+      auto Y = y.Write();
+      // y = D^{-1} x
+      MFEM_FORALL(i, N, Y[i] = DI[i] * X[i]; );
    }
-   auto DI = dinv.Read();
-   auto R = residual.Read();
-   auto Y = y.ReadWrite();
-   MFEM_FORALL(i, N, Y[i] += DI[i] * R[i]; );
 }
 
 OperatorChebyshevSmoother::OperatorChebyshevSmoother(Operator* oper_,