New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Add a Tpetra-based version of Trilinos SparseMatrix and SparsityPattern. #16288
Conversation
Absolutely! |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I realize that I have a bunch of pending comments, likely from quite a long time ago already. Might as well flush them.
What is the status of this patch? Do you want us to look through it in its entirety?
#include <deal.II/lac/trilinos_tpetra_sparse_matrix.h> | ||
#include <deal.II/lac/trilinos_tpetra_vector.h> |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Since you make no other changes to this file, do you actually need these #include
s?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
That originates from some changes that are not included in this pull request.
I removed those includes.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Edit: These changes are actually needed. Otherwise, affine_constraints.inst
does not compile.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Could these lines be move to the .templates.h
file then? Or ideally to the .cc
file? This way they don't have to be included transitively by everyone who just wants the affine_constraints.h
file.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
They are already in the .templates.h
file
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
My mistake -- I now realize that I'm looking at the .templates.h file, not the .h file.
Thank you for your feedback so far! |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Here is batch #1. Nothing of great importance -- most of what you did was just follow the old interface, which is useful because there is little to review then. I'll get to the rest next.
using size_type = dealii::types::signed_global_dof_index; | ||
|
||
template <typename Node> | ||
using MapType = | ||
Tpetra::Map<int, dealii::types::signed_global_dof_index, Node>; | ||
|
||
template <typename Number, typename Node> | ||
using MatrixType = Tpetra:: | ||
CrsMatrix<Number, int, dealii::types::signed_global_dof_index, Node>; | ||
|
||
template <typename Node> | ||
using GraphType = | ||
Tpetra::CrsGraph<int, dealii::types::signed_global_dof_index, Node>; |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
MatrixType
depends on the scalar type, and these are all strictly speaking copies of what is used inside the SparseMatrix
class. Could we avoid making these using
declarations at the global level? I'll write a patch for the existing class as well.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I am not sure what I should change here.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I've made it down to operator/=
in the sparse_matrix.templates.h
file. I'll take care of the rest next year :-)
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
This is the rest of the sparse matrix stuff. This is looking good, I think I only found one real bug.
I will continue with the sparsity pattern things later.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Here is the sparsity pattern stuff. I think that most of it is right, but there are a few places where I'm unsure.
Teuchos::RCP<GraphType> | ||
trilinos_sparsity_pattern() const; |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
You can return a const Teuchos::RCP<...> &
instead of the copy. This avoids creating a copy and thereby playing with the use count, avoiding a couple of nanoseconds of work :-)
You can probably do the same in other places where you return the underlying RCP.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
This is against the idea of using Teuchos::RCP<...>
I don't think a few nanoseconds are worth it...
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I thought of it more as a question of what pattern we want to use in general. For all non-scalar types, getter functions that are const
should return a const
reference. It may not matter in this specific case, but I think as a general rule it is the right thing to do and it would be nice if we did the same here :-)
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
We could return a Teuchos::RCP<const ...>
in those places.
Teuchos::RCP<const MapType> | ||
domain_partitioner() const; |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Like here and in the following function.
Thanks for going through all of my comments. Once you think it's in good shape, squash into a small number of commits and let me know -- then we can merge! |
2af3047
to
f0eaee6
Compare
5d6ddf2
to
0b2b5f0
Compare
I answered most of the comments. And most of them are already marked as resolved. Maybe one remark: As discussed in the comments, the |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Let's see how this goes!
/rebuild |
Alright, let's get this merged and see what happens! |
This is the minimal set for a SparseMatrix implementation, where my goal is to solve the Laplace equation from step-40. This is a rather large pull request and hopefully a big step towards a working Tpetra interface.
But first, why a new class instead of reworking the existing SparseMatrix class?
I know the original plan was to rework the existing TrilinosWrappers::SparseMatrix class to use Tpetra instead of Epetra. But after I started to work on that (see PR: #16052), I quickly reached a point where it was not feasible to continue with that approach.
My three main reasons to write a new class instead are:
Two template Arguments
As mentioned above, part of the reason for introducing a new class instead of reworking the existing one is to expose the template parameter for the scalar type, so the TpetraWrappers::SparseMatrix class can also be used with different scalars (similar to other existing SparseMatrix classes in deal.II).
The second template parameter exposes the Kokkos Node. By choosing different Kokkos Nodes for the SparseMatrix class, the code can be easily adapted to be executed on other host devices (e.g., GPUs). More work must be done to use that feature to its full potential. But we have to start somewhere.
To make the class still as similar as possible to existing SparseMatrix classes in deal.II, a reasonable default value for that template parameter is set.
I would like your opinion, especially on the (Kokkos) Node template parameter.
Included functions in the SparseMatrix class
My first intermediate goal is to use the Tpetra wrapper to solve step-40. Therefore, I focused on the functions used in step-40. Anyhow, I tried to include the set of the most important functions. If I missed something that is frequently used, please let me know!
About the Instantiation:
I only implemented a rudimentary instantiation for the TpetraWrappers::SparseMatrix class class. @jpthiele implemented an extensive instantiation for this class (see tpetra_wrappers_inst). Once this PR is accepted, he will create a follow-up PR.
@masterleinad, especially the
TpetraWrappers::SparsityPattern
, is based on your work from here. Is it okay for you that I have included that in this pull request?So please take a close look at that class, and I look forward to your feedback!