message: port jacobi2 and jacobi3 to device
diff: @@ -2461,52 +2461,60 @@ void SparseMatrix::DiagScale(const Vector &b, Vector &x, double sc) const
  });
  }
  -void SparseMatrix::Jacobi2(const Vector &b, const Vector &x0, Vector &x1,
  - double sc) const
  +template <bool useFabs>
  +static void JacobiDispatch(const Vector &b, const Vector &x0, Vector &x1,
  + const Memory<int> &I, const Memory<int> &J,
  + const Memory<double> &A, const int height,
  + const double sc)
  {
  - MFEM_VERIFY(Finalized(), "Matrix must be finalized.");
  + const bool useDevice = b.UseDevice() || x0.UseDevice() || x1.UseDevice();
  - for (int i = 0; i < height; i++)
  + const auto bp = b.Read(useDevice);
  + const auto x0p = x0.Read(useDevice);
  + auto x1p = x1.Write(useDevice);
  +
  + const auto Ip = Read(I, height+1, useDevice);
  + const auto Jp = Read(J, J.Capacity(), useDevice);
  + const auto Ap = Read(A, J.Capacity(), useDevice);
  +
  + MFEM_FORALL_SWITCH(useDevice, i, height,
  {
  - double resi = b(i), norm = 0.0;
  - for (int j = I[i]; j < I[i+1]; j++)
  + double resi = bp[i], norm = 0.0;
  + for (int j = Ip[i]; j < Ip[i+1]; j++)
  {
  - resi -= A[j] * x0(J[j]);
  - norm += fabs(A[j]);
  + resi -= Ap[j] * x0p[Jp[j]];
  + if (useFabs)
  + {
  + norm += fabs(Ap[j]);
  + }
  + else
  + {
  + norm += Ap[j];
  + }
  }
  if (norm > 0.0)
  {
  - x1(i) = x0(i) + sc * resi / norm;
  + x1p[i] = x0p[i] + sc * resi / norm;
  }
  else
  {
  - MFEM_ABORT("L1 norm of row " << i << " is zero.");
  - }
  + MFEM_ABORT_KERNEL("L1 norm of row is zero.");
  }
  + });
  }
  -void SparseMatrix::Jacobi3(const Vector &b, const Vector &x0, Vector &x1,
  +void SparseMatrix::Jacobi2(const Vector &b, const Vector &x0, Vector &x1,
  double sc) const
  {
  MFEM_VERIFY(Finalized(), "Matrix must be finalized.");
  -
  - for (int i = 0; i < height; i++)
  - {
  - double resi = b(i), sum = 0.0;
  - for (int j = I[i]; j < I[i+1]; j++)
  - {
  - resi -= A[j] * x0(J[j]);
  - sum += A[j];
  - }
  - if (sum > 0.0)
  - {
  - x1(i) = x0(i) + sc * resi / sum;
  + JacobiDispatch<true>(b,x0,x1,I,J,A,height,sc);
  }
  - else
  +
  +void SparseMatrix::Jacobi3(const Vector &b, const Vector &x0, Vector &x1,
  + double sc) const
  {
  - MFEM_ABORT("sum of row " << i << " is zero.");
  - }
  - }
  + MFEM_VERIFY(Finalized(), "Matrix must be finalized.");
  + JacobiDispatch<false>(b,x0,x1,I,J,A,height,sc);
  }
  void SparseMatrix::AddSubMatrix(const Array<int> &rows, const Array<int> &cols,
  