Template TpetraWrappers::SparseMatrix::vmult on the vector type. #16614

kinnewig · 2024-02-09T10:53:44Z

Rework from LinearAlgebra::TpetraWrappers::SparseMatrix<Number, MemorySpace>::vmult (and Tvmult, vmult_add, Tvmult_add) such that it takes arbritraty vector types.

Initially, this was part of #16588, but I moved it into a separate pull request because it contains a few points I would like to discuss.

Each function of vmult, Tvmult, vmult_add, and Tvmult_add contains a part that converts the raw data of the input vectors into a Kokkos::DualView. I was thinking about moving that conversion into a separate (internal) function (and marking that function as inline).
And this is the most convenient way to convert the raw data from a vector into a Kokkos::DualView I came up with. But if there is a more elegant way to do that, please let me know!

Furthermore, we need a more elegant way to initialize all templates.

Since some tests include a vmult between a TrilinosWrappers::SparseMatrix and a dealii::Vector, this should help with #16611 as well.

bangerth

I will leave this one to @masterleinad given that it does all of this view magic...

bangerth · 2024-02-09T23:10:24Z

include/deal.II/lac/trilinos_tpetra_sparse_matrix.templates.h

+      template <typename Number, typename MemorySpace>
+      Number *
+      begin(LinearAlgebra::TpetraWrappers::Vector<Number, MemorySpace> &V)
+      {
+        return V.trilinos_vector().getDataNonConst().get();
+      }
+
+      template <typename Number, typename MemorySpace>
+      const Number *
+      begin(const LinearAlgebra::TpetraWrappers::Vector<Number, MemorySpace> &V)
+      {
+        return V.trilinos_vector().getData().get();
+      }
+
+      template <typename Number, typename MemorySpace>
+      Number *
+      end(LinearAlgebra::TpetraWrappers::Vector<Number, MemorySpace> &V)
+      {
+        return V.trilinos_vector().getDataNonConst().get() +
+               V.trilinos_vector().getLocalLength();
+      }
+
+      template <typename Number, typename MemorySpace>
+      const Number *
+      end(const LinearAlgebra::TpetraWrappers::Vector<Number, MemorySpace> &V)
+      {
+        return V.trilinos_vector().getData().get() +
+               V.trilinos_vector().getLocalLength();
+      }


Instead of providing these overloads here, why not just add begin() and end() functions to the LinearAlgebra::TpetraWrappers::Vector class? Then you won't need any of these functions here.

Sure, I will change this!
But do all other vector classes In deal.II implement begin() and end() ?

I saw, at least dealii::vector<Number> has a begin() and end() method implemented.

Good question. I don't actually know whether all vectors have begin() and end(), but even if they did it's not clear what semantics that would follow: Only the locally owned? locally stored? all elements?

Perhaps what you have is better. In fact, perhaps it would be better to call these functions begin_locally_owned() and end_locally_owned() to make this clear.

masterleinad · 2024-02-10T01:46:01Z

include/deal.II/lac/trilinos_tpetra_sparse_matrix.templates.h

+      Kokkos::View<
+        typename VectorType::impl_scalar_type **,
+        Kokkos::LayoutLeft,
+        Kokkos::Device<typename MemorySpace::kokkos_space::execution_space,
+                       typename MemorySpace::kokkos_space>>
+        kokkos_view_src(raw_src, src_local_size);


Do you really need a two-dimensional View here?

Suggested change

Kokkos::View<

typename VectorType::impl_scalar_type **,

Kokkos::LayoutLeft,

Kokkos::Device<typename MemorySpace::kokkos_space::execution_space,

typename MemorySpace::kokkos_space>>

kokkos_view_src(raw_src, src_local_size);

Kokkos::View<

typename VectorType::impl_scalar_type *,

typename MemorySpace::kokkos_space>

kokkos_view_src(raw_src, src_local_size);

How do you know the memory space of InputVectorType?

How do you know the memory space of InputVectorType?

Good point!

I thought I had to provide the memory space of the target vector. But you are right; this is not true. We need to provide the memory space of the source vector. So we need to be way more careful here...

For most vector types, the memory space is always Host. The only exceptions I am aware of is the TpetraWrappers::Vector<Number, MemorySpace> and the distributed::Vector<Number, MemorySpace>.

Do you really need a two-dimensional View here?

How do you know the memory space of InputVectorType?

The Tpetra::Vector always expects a two-dimensional View, as far as I know.

OK. In this case, you still have to provide extents for both dimensions in the View constructor.

OK. In this case, you still have to provide extents for both dimensions in the View constructor.

I am not that experienced with Kokkos::View, what constructor should I use here?

It should be

Kokkos::View< typename VectorType::impl_scalar_type **, Kokkos::LayoutLeft, Kokkos::Device<typename MemorySpace::kokkos_space::execution_space, typename MemorySpace::kokkos_space>> kokkos_view_src(raw_src, src_local_size, 1);

so that you provide the pointer and an argument for the extent in every dimension. This should be strictly enforced in newer Kokkos versions. Also, there is not much of a point to provide a full Device template argument. In the end, we only need the MemorySpace. For managed Views you could provide an execution space instance as runtime argument to make sure it uses a particular stream (and fences less) or even decide not to initialize the allocation at all. None of that really matters for unmanaged Views like here, of course.

masterleinad · 2024-02-10T16:10:36Z

include/deal.II/lac/trilinos_tpetra_sparse_matrix.templates.h

+      // For the src vector:
+      // create a Kokkos::View from the src vector
+      typename VectorType::impl_scalar_type *raw_src =
+        reinterpret_cast<typename VectorType::impl_scalar_type *>(


Why do you need a reinterpret_cast here?

It seems like the Tpetra::Vector expects this type.
Since this is just a complicated alias for Number, I can give it a try and remove this cast.

I would be very suspicious of any reinterpret_cast we are doing in these classes since we don't really have a buffer that we would use with placement new. Whenever we need a different type we should copy values instead.

I agree with this. x.begin() does not always return a pointer. It could just be an iterator. reinterpret_cast will in such cases not result in anything useful, with no error message. I think there ought to be a better way.

What are the types you're casting here?

We do not need the reinterpret_cast here. But at another point, we need a const_cast.
I specialized the function so we know which vector type we are working on, so we can ensure that the cast is working correctly.

kinnewig · 2024-02-15T14:50:44Z

I reworked the vmult function.
It now calls an internal function such that one can define specializations more easily. When no specialization for the InputVectorType is provided, it throws an Assert. I guess this is the best approach since it ensures that the memory space is handled correctly

include/deal.II/lac/trilinos_tpetra_sparse_matrix.templates.h

masterleinad · 2024-02-22T14:47:11Z

include/deal.II/lac/trilinos_tpetra_sparse_matrix.templates.h

+        typename SparseMatrix<Number, MemorySpace>::VectorType tpetra_src(
+          M.trilinos_matrix().getDomainMap(), kokkos_dual_view_src);
+
+        M.trilinos_matrix().apply(tpetra_src, tpetra_dst, mode, alpha, beta);


Does this also sync the DualView?

I just assumed it does. I just verified it passed the test:

trilinos_tpetra/sparse_matrix_vector_08.cc

Since it passed that test, I assumed it was okay as it was.

If there is only one memory space, sync doesn't do anything anyway. I think it's reasonable to assume that apply makes sure the two Views are in sync, though.

masterleinad · 2024-02-26T13:28:01Z

include/deal.II/lac/trilinos_tpetra_sparse_matrix.templates.h

+        typename SparseMatrix<Number, MemorySpace>::VectorType tpetra_src(
+          M.trilinos_matrix().getDomainMap(), kokkos_dual_view_src);
+
+        M.trilinos_matrix().apply(tpetra_src, tpetra_dst, mode, alpha, beta);


If there is only one memory space, sync doesn't do anything anyway. I think it's reasonable to assume that apply makes sure the two Views are in sync, though.

include/deal.II/lac/trilinos_tpetra_sparse_matrix.templates.h

masterleinad · 2024-02-26T14:00:09Z

You'll have to fix the indentation.

diff --git a/include/deal.II/lac/trilinos_tpetra_sparse_matrix.templates.h b/include/deal.II/lac/trilinos_tpetra_sparse_matrix.templates.h
index 87a3723e0..6ba2d5613 100644
--- a/include/deal.II/lac/trilinos_tpetra_sparse_matrix.templates.h
+++ b/include/deal.II/lac/trilinos_tpetra_sparse_matrix.templates.h
@@ -83,9 +83,8 @@ namespace LinearAlgebra
           dst.begin(), dst_local_size, 1);
 
         // get a Kokkos::DualView
-        auto mirror_view_dst =
-          Kokkos::create_mirror_view_and_copy(typename MemorySpace::kokkos_space{},
-                                     kokkos_view_dst);
+        auto mirror_view_dst = Kokkos::create_mirror_view_and_copy(
+          typename MemorySpace::kokkos_space{}, kokkos_view_dst);
         typename SparseMatrix<Number, MemorySpace>::VectorType::dual_view_type
           kokkos_dual_view_dst(mirror_view_dst, kokkos_view_dst);
 
@@ -99,9 +98,8 @@ namespace LinearAlgebra
           const_cast<Number *>(src.begin()), src_local_size, 1);
 
         // get a Kokkos::DualView
-        auto mirror_view_src =
-          Kokkos::create_mirror_view_and_copy(typename MemorySpace::kokkos_space{},
-                                     kokkos_view_src);
+        auto mirror_view_src = Kokkos::create_mirror_view_and_copy(
+          typename MemorySpace::kokkos_space{}, kokkos_view_src);
         typename SparseMatrix<Number, MemorySpace>::VectorType::dual_view_type
           kokkos_dual_view_src(mirror_view_src, kokkos_view_src);

Otherwise, this looks good to me now.

kinnewig · 2024-02-26T15:32:01Z

Thank you for all the feedback!

I fixed the indentation and squashed everything into one commit.

masterleinad · 2024-02-26T19:54:19Z

/rebuild

kinnewig · 2024-02-28T08:32:29Z

I fixed the failing test.

bangerth reviewed Feb 9, 2024

View reviewed changes

jpthiele mentioned this pull request Feb 10, 2024

Add SolverDirect for Tpetra Wrappers #16602

Merged

4 tasks

masterleinad reviewed Feb 10, 2024

View reviewed changes

kinnewig requested a review from masterleinad February 20, 2024 07:41

masterleinad requested changes Feb 22, 2024

View reviewed changes

kinnewig force-pushed the tpetra_wrappers_vmult branch from 2aed16a to dda3577 Compare February 25, 2024 10:35

masterleinad reviewed Feb 26, 2024

View reviewed changes

masterleinad approved these changes Feb 26, 2024

View reviewed changes

kinnewig force-pushed the tpetra_wrappers_vmult branch from 11ca285 to 3874de7 Compare February 26, 2024 15:20

kinnewig force-pushed the tpetra_wrappers_vmult branch from 3874de7 to c205c21 Compare February 26, 2024 19:17

masterleinad added the ready to test label Feb 26, 2024

Template TpetraWrappers::SparseMatrix::vmult on the vector type.

9db41cd

kinnewig force-pushed the tpetra_wrappers_vmult branch from c205c21 to 9db41cd Compare February 28, 2024 08:31

masterleinad merged commit 1d56f85 into dealii:master Feb 28, 2024
16 checks passed

tamiko mentioned this pull request Feb 28, 2024

Regression tester regressed 1d56f8 #16701

Closed

kinnewig deleted the tpetra_wrappers_vmult branch March 18, 2024 18:18

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Template TpetraWrappers::SparseMatrix::vmult on the vector type. #16614

Template TpetraWrappers::SparseMatrix::vmult on the vector type. #16614

kinnewig commented Feb 9, 2024 •

edited

bangerth left a comment

bangerth Feb 9, 2024

kinnewig Feb 11, 2024

bangerth Feb 11, 2024

masterleinad Feb 10, 2024

kinnewig Feb 11, 2024 •

edited

kinnewig Feb 11, 2024

masterleinad Feb 11, 2024

kinnewig Feb 15, 2024

masterleinad Feb 15, 2024

masterleinad Feb 10, 2024

kinnewig Feb 11, 2024

masterleinad Feb 11, 2024

bangerth Feb 11, 2024

kinnewig Feb 15, 2024

kinnewig commented Feb 15, 2024 •

edited

masterleinad Feb 22, 2024

kinnewig Feb 25, 2024

masterleinad Feb 26, 2024

masterleinad Feb 26, 2024

masterleinad commented Feb 26, 2024

kinnewig commented Feb 26, 2024

masterleinad commented Feb 26, 2024

kinnewig commented Feb 28, 2024

Template TpetraWrappers::SparseMatrix::vmult on the vector type. #16614

Template TpetraWrappers::SparseMatrix::vmult on the vector type. #16614

Conversation

kinnewig commented Feb 9, 2024 • edited

bangerth left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

kinnewig Feb 11, 2024 • edited

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

kinnewig commented Feb 15, 2024 • edited

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

masterleinad commented Feb 26, 2024

kinnewig commented Feb 26, 2024

masterleinad commented Feb 26, 2024

kinnewig commented Feb 28, 2024

kinnewig commented Feb 9, 2024 •

edited

kinnewig Feb 11, 2024 •

edited

kinnewig commented Feb 15, 2024 •

edited