Add 'requires' clauses to PETSc solvers.

[bangerth]

@stefanozampini FYI

Relates to #14840.

[stefanozampini]

Wouldn't be better to rename this to can_wrap_petsc_vector_type or something?

[stefanozampini]

this is to reflect the fact that the supported classes must have the ability to construct new instances by wrapping a Vec and can return the Vec via the petsc_vector() method.

[stefanozampini]

same for the matrices

[stefanozampini]

A naive question: can this be done using traits? or with enable_if checking for explicit type::type(Vec) and Vec type::petsc_vector()?

[bangerth]

Yes, this can be done. The most explicit approach would be to do something like this:

template <typename VectorType>
concept can_wrap_petsc_vec =
  requires {
    VectorType vec_object (std::declval<Vec>());  // checks that an object can be constructed from a Vec object
  };

and then to use this concept in the same place as where I'm currently already using the other concept.

But it turns out that others have found such a concept useful in the past, and so it is also std::constructible_from<VectorType,Vec> :-)

[stefanozampini]

thanks. So the whole point of std::concepts is to make the error message from the compilation stage readable? (and not a million lines long?)

[bangerth]

That's basically the idea. For example, if you have

  template <typename MeshType>
  int n_cells (const MeshType &m) {
    return m.size();
  }

then you will get awkward error messages if you call that function with an int argument -- something like error: 'int' variable 'm' does not have a member function 'size()'.

So you want to annotate the function with a "concept" like

  template <typename MeshType>
    requires is_mesh_type<MeshType>
  int n_cells (const MeshType &m) {
    return m.size();
  }

that encodes the kinds of properties the type MeshType has to have to make this function useful. If you call n_cells(1), the compiler will say something like error: can't call 'n_cells()' with argument 'MeshType=int' because 'int' does not satisfy the concept 'is_mesh_type'.

That is how I use concepts in this patch. There are other uses for concepts, but that's a different story. The idea for these other uses are of the kind where you want to use different algorithms based on properties of the involved types. Say something like this:

  template <typename MatrixType>
    requires (matrix_type_represents_symmetric_matrix<MatrixType> == false)
  std::pair<FullMatrix,FullMatrix>
  lu_decomposition (const MatrixType &m)
  {
    ...use LU algorithm...
  }

  template <typename MatrixType>
    requires (matrix_type_represents_symmetric_matrix<MatrixType> == true)
  std::pair<FullMatrix,FullMatrix>
  lu_decomposition (const MatrixType &m)
  {
    ...use Cholesky algorithm...
  }

Basically, concepts here serve the purpose of "dispatching" to overloaded versions of the same function name, based on whether or not the type satisfies certain properties or not. This was previously possible with things such as std::enable_if, but it was awkward at best.

[bangerth]

Take another look. I've chosen to keep the original constraint, but add an "or can be constructed from a Vec object" as an alternative. That second part is, strictly speaking, making the first part redundant -- but I think it is easier to read this way from a user perspective.

If you allow me to, I will take care of similar constraints for the matrix objects in a second step.

[stefanozampini]

bug in the formatter?

[bangerth]

Possibly -- who knows. clang-format is a thing unto its own...

[stefanozampini]

If you allow me to, I will take care of similar constraints for the matrix objects in a second step.

I allow :-)

[bangerth]

Ping?

[stefanozampini]

sorry, I taught I already approved

[peterrum]

Looks good. But is a bit annoying that you need to write the requires clauses to all the implementation. This artificially increases the number of lines of codes but without adding anything new...

[bangerth]

@peterrum Yes, I agree. I did not originally think that it was necessary to do that (and GCC 11.2 doesn't force you to do it either), but I ultimately realized that it is necessary. Take this code:

  template <int dim>
    requires (dim % 2 == 0)             // only for even dimensions
  class Triangulation {
    void create();
  };

  template <int dim>
    requires (dim % 2 == 1)             // only for odd dimensions
  class Triangulation {
    void create();
  };

  template <int dim>
  void Triangulation<dim>::create() {...};  // which of the two classes above does this relate to???

So I ended up understanding why the language requires duplicating the requires clauses, even though you're right that it seems superfluous and verbose.