Add an abstract base class ReadVector purely for vector access.

[drwells]

This is at the proof-of-concept stage and a viable direction for part of #15168.

In FEValues (and places like the Kelly error estimator and integrate_difference in which we only access vectors through FEValues) the only vector operation we need is "extract values from a global vector on a single cell". We could eliminate the use of VectorType completely in that class by making all of our vectors implement a single common interface for that task. I implemented this as new base class ReadVector<T> which contains virtual functions size() and extract_entries().

FEValues is one of the most expensive things to compile so reducing the number of instantiations in it by about an of magnitude (by switching from instantiation over all vector classes to instantiation over all scalars) would help a lot.

One disadvantage to this approach is that we need to either delete some checks on vector sizes or make size() virtual.

[drwells]

Another thing ReadVector could do is set up a Kokkos view, e.g., via

template <class MemorySpace>
PetscErrorCode VecGetKokkosView(Vec, Kokkos::View<const PetscScalar *, MemorySpace> *);

or by Tpetra::MultiVector::getLocalViewDevice().

[drwells]

Any thoughts? If people are OK with making

template <typename Number>
class ReadVector
{
public:
  using size_type = types::global_dof_index;

  virtual size_type
  size() const = 0;

  virtual void
  extract_elements(const ArrayView<const types::global_dof_index> &indices,
                   ArrayView<Number> &entries) const = 0;
};

a base class for all of our vector classes we could get rid of a lot of template instantiations - I'd like to get this done for 9.5. We can work on the Kokkos view functions and other advanced features at the workshop.

[Rombur]

Why are you using an ArrayView for the indices and not an IndexSet?

[drwells]

I used ArrayView since most of our vector classes already implement something like (e.g., for Vector)

  template <typename OtherNumber>
  void
  extract_subvector_to(const std::vector<size_type> &indices,
                       std::vector<OtherNumber> &    values) const;

I don't know if that's any better or worse than an IndexSet.

[Rombur]

The issue is that if you want all the locally owned elements, you first need to create large std::vector and then, you can create the ArrayView.

[bangerth]

What was your reasoning not to have a VectorView class? This would be a type-erasure class of the kind

class VectorView 
{
    class ViewBase {
       public:
          virtual unsigned int size() const = 0;  // to show only one of the member functions this class will have
    };

    template <typename VectorType>
    class View {
       public:
          View (VectorType &v) : v(v) {}
          virtual unsigned int size() const { return v.size(); }
       private:
          VectorType &v;
    };

    std::unique_ptr<ViewBase> view;

    public:
      // constructor; allows implicit initialization
      template <typename VectorType>
      VectorView (VectorType &v) : view(std::make_unique<View<VectorBase>>>(v) {}

      unsigned int size() { return view->size(); }
};

This has the advantage that you don't have to adjust the interface of the vector classes, including the fact that the size() function does not become virtual. I don't know that that particular issue makes a performance difference, but I like the fact that we wouldn't have to change anything about the existing vector classes.

[drwells]

The problem I couldn't figure out how to fix with VectorView is the one I described in #15168 (comment): If we have a function

template <typename Number>
void foo(const VectorView<Number> &v);

then AFAICT we cannot do

Vector<double> fe_vector;
foo(fe_vector);

since template resolution occurs before function overloading. Here's a complete example which doesn't use deal.II:

#include <iostream>
#include <vector>

template <typename Number>
struct Vector
{
  std::vector<Number> values;
};

template <typename Number>
struct View
{
  Number *a;

  View(Vector<Number> &b) : a(b.values.data())
  {}
};

void foo1(const View<double> &view)
{
    std::cout << "foo1()" << std::endl;
}

template <typename Number>
void foo2(const View<Number> &view)
{
    std::cout << "foo2()" << std::endl;
}

int main()
{
    Vector<double> v;
    foo1(v); // OK
    foo2(v); // error
}

fails to compile with the error

[drwells@archway overload]$ g++ overload.cc
overload.cc: In function 'int main()':
overload.cc:34:9: error: no matching function for call to 'foo2(Vector<double>&)'
   34 |     foo2(v); // error
      |     ~~~~^~~
overload.cc:25:6: note: candidate: 'template<class Number> void foo2(const View<Number>&)'
   25 | void foo2(const View<Number> &view)
      |      ^~~~
overload.cc:25:6: note:   template argument deduction/substitution failed:
overload.cc:34:9: note:   'Vector<double>' is not derived from 'const View<Number>'
   34 |     foo2(v); // error
      |     ~~~~^~~

Writing a vector view class in that way would be great but I don't see how to make it work - hence I went with the simpler version of just writing a base class.

[bangerth]

Ah, yes, that's a bummer. I guess that's why we have make_array_view() instead of relying on automatic conversions. But here we really don't want to use make_vector_view() in all of the places where we want to use a view :-( I tried to work around this using an explicit conversion operator in class Vector, but that also doesn't seem to work. I also tried a C++17-style deduction guide:

template <typename Number>
View(const Vector<Number> &) -> View<Number>;

But that, too, doesn't work unless you want to write

foo2(View(v));

which works but defeats the purpose :-(

[bangerth]

So what do others think? I don't think it would be too big of a problem to make size() virtual as well.

[drwells]

size() is already virtual for most of the parallel vector classes so I don't think it's that significant of a change.

The idea of a vector view is still the right way, IMO, to implement operator[] and local_element() with type erasure, so we should still do it in some other way (I especially like that we can do Kokkos views).

[drwells]

w.r.t. using IndexSet or std::vector - we already have the std::vector-based function for access to values corresponding to cell dofs in our vector classes. We can just add both to ReadVector to cover both use cases.

[drwells]

/rebuild

This patch lowers the total number of symbols in the library from 217k to 209k - another advantage is that it decouples FEValues and the other classes from all of our vector classes.

[bangerth]

What do we want to do with this patch? I love the direction, including that you've taken the time to also convert places such as KellyErrorEstimator that also have a large number of instantiations. At the same time, this is a large change, and I'd be more comfortable if we pushed that back to after the release.

[masterleinad]

I would prefer to discuss this after the release since this is not entirely uncontroversial.

[drwells]

Since its part of #745 etc. I think we should discuss it at the workshop.

This patch shows that we can support FEValues and things that only use FEValues with just two vector functions - knowing that can help inform our choices for other changes to the vector classes.

[Rombur]

What about instead of rolling our own class, we use std::mdspan(or std::span)? mdspan is used in the interface of std::linalg. Under the hood, the Kokkos::View will switch to mdspan soonish so it will be easy to get the mdspan associated with a Kokkos::View. std::mdspan is part of c++23 but here is an implementation that is compatible with earlier standard. It's a header-only library, so it's easy to include.

[drwells]

As an alternative to the proposed View class? We could, but I think we need RAII to handle something like PETSc where we need matching calls to VecGetArrayRead() / VecRestoreArrayRead(). mdspan can definitely be a part of that.

[Rombur]

We need to allow vectors to translate from global to local indices to resolve this.

We should be able to do that through the access policy of mdspan.

[drwells]

Do all of our dependencies use contiguous storage for vectors? PETSc does (we convert to a serial vector) but I don't know about the other vector classes.

[Rombur]

Do all of our dependencies use contiguous storage for vectors?

Trilinos and our own vector do use contiguous storage.

[bangerth]

But then we still have to convert manually from these vectors to std::span. That runs into the same problem we had with the VectorView. The only automatic conversion was to the base class.

[kronbichler]

I am looking forward to the workshop (albeit I only participate remotely and in some time slots). My view on this is that having a simple way to access these functions would be great. I understand the base class, so that would be good for me; maybe mdspan can also do it, but I do not get those details without seeing the code.

In either case, the performance impact should be small or negligible. We go into these functions with global indices, so they will be much slower than direct array access, and thus another indirection for size() and the extraction of elements through a virtual function will not matter. The only concern I have is that the vector should not be copied upon access (which is out of question anyway).

[luca-heltai]

@stefanozampini, I think this is one of the discussions @bangerth was referring to.

[tjhei]

I like the idea and the reduction in the number of instantiations is great.

As mentioned earlier, something like VectorBase seems to be a more fitting name. Who knows what other functions will end up in this base later?

[drwells]

It sounds like everyone is happy with the current state of this - in the meeting this morning we decided to keep the name ReadVector.

[kronbichler]

Are the changes in these integrators the only places that need to be touched? Do you happen to know why here? Is it because we have std::vector<std::vector<double>> that can somehow not be automatically converted? I think that might be something we need to describe more clearly in the incompatibilities section of a changelog?

[drwells]

It's because the template instantiation is different: since we instantiate on Number instead of VectorType the deduction process is slightly different and we need to explicitly pass an ArrayView (instead of relying on the implicit conversion constructor). I added a changelog entry summarizing the problem.

I reran all the tests which use FEValues or FEFaceValues locally and they all pass so I think I fixed all the relevant ones in this PR.

[bangerth]

There we go!

[bangerth]

This is the change that caused #15580.

[bangerth]

For reference, on my laptop:

Before:

-rwxrwxr-x 1 bangerth bangerth 1291636280 Jul  1 16:19 libdeal_II.g.so.9.6.0-pre
-rwxrwxr-x 1 bangerth bangerth  367992504 Jul  1 16:20 libdeal_II.so.9.6.0-pre

After:

-rwxrwxr-x 1 bangerth bangerth 1240677608 Jul  2 14:47 libdeal_II.g.so.9.6.0-pre
-rwxrwxr-x 1 bangerth bangerth  343967440 Jul  2 14:48 libdeal_II.so.9.6.0-pre

So about a 4% reduction for the debug library, and 6.5% for the release library. Nice!