Distributed vector

[MarcelKoch]

This PR adds a distributed (multi-)vector class. I've opted into using mostly multivector based descriptions instead of using dense for that. This might clash a bit with the rest, until we separate dense into multivector and dense matrix.

The vector is a distributed global vector, i.e. it has the size global_num_rows x num_cols. Only the rows are distributed. Each process stores only the rows that belong to itself in Dense matrix, according to the used partition.

Regarding the tests, I've put the tests that require MPI into /core/test/mpi/distributed, not sure if that is the best place.

Partially addresses #907

Note: should we use distributed-develop as target branch again? That way we could make sure that we have all basic things (vectors, matrices, solvers) together for the next release.

Main contributions are from @upsj and @pratikvn.

[MarcelKoch]

format-rebase!

[ginkgo-bot]

Error: Rebase failed, see the related Action for details

[MarcelKoch]

format!

[upsj]

A first look over the code. The class structure and implementation looks good to me, haven't checked the tests yet.

A fundamental discussion point I would like to raise is whether Vector should store Partition. In my original design, both Matrix and Vector exclusively used local indexing after they have been filled, and I believe this is a property we could use generally. If users need to translate global indices and local indices, they have the original partition available from their setup code. Using the vector and matrix as a convenient vehicle to carry them around seems like it might overload the functionality a bit.

[MarcelKoch]

Regarding the partition, my use case for that was repartitioning, which we should support in the long term. This could be achieved without storing the partition in the vector class, since the repartition would need to store it anyway. However, I don't see how we could make sure that the vector we are repartitioning is valid for that operation, if the vector does not store its partition. I think the additional safety we get from storing it is worth the downside. (I'm not entirely sure what you see as a downside, except that it is not necessary, perhaps you could elaborate that?)

Also, I think there is a strong connection between vector and partition. The vector is only valid in combination with this one partition, which is expressed in that way.

[pratikvn]

Mostly looks good. Some minor questions on function interfaces.

[pratikvn]

Regarding the base class to merge this into, distributed-develop might be a good choice, unless we decide that the interface more or less looks good, then we can try to merge this directly into develop

I think storing a partition in Vector also allows us to only Vector operations on data (with a different partition), without access to matrix or other partition storing classes. Therefore, I also think it would be useful to store the partition in the vector class.

[MarcelKoch]

I've worked in most of the reviews. There are two open questions so far (1 minor, 1 major):

1. Do agree on the current order of the constructor parameters and their default values?
2. Do we keep the partition in the vector?

For 2. the arguments against (as I understand it) are:

• There is no use case where the partition information is required, but won't be provided elsewhere. This especially holds for the SpMV

The arguments in favor are:

• Can ensure that we always operate on consistent data, e.g. SpMV. For example, we could have cases, where the local and global sizes all match up, but the partitions are not the same. (Note: this would require testing partitions for equality in constant time)
• Simple access to the partition, without relying on other classes that might store it.

[MarcelKoch]

I've added the DenseCache that is used in the distributed matrix also to this PR. If wanted, I could also extract that into its own PR, since the changes for that are quite small.

[pratikvn]

Mostly looks good. Some comments on the whether we want to make DenseCache public. Regarding the discussion points:

1. I think the constructor order is okay.
2. I think a distributed object should have always have partition member.
   • It should be something like a Concept require for distributed objects.
   • There are cases for vector operations, where the partition makes sense. For example dot products between two vectors.
   • Additionally, the read function already takes a partition. In my opinion, it might make more sense to not have the read function take a partition, but create the object with a partition and the read function use the object's partition.

[MarcelKoch]

@pratikvn you raise a very important point. The current implentations are only correct, if both vectors (think axpy, dot, ...) use the same partition. It is not enough that both local vectors have the same size. This is something that we should check for at some point.

[upsj]

I see, there seems to be a very fundamental difference between our approaches. In my original implementation, neither Matrix nor Vector store their partition. If you ignore the off-diagonal entries, you can imagine the matrix as a block-diagonal matrix with local entries in each row. I agree that having the partition available is useful in some cases, but in that case the current implementation has similar issues to what #753 is fixing: Copying the vector to another executor doesn't copy the Partition, which leads to issues when trying to read it from a kernel.

A very basic rule I try to follow in my designs is to try and keep the amount of state of each object as minimal as possible, and always prefer members with value semantics. This leads to the least weird empty states for objects. For communicator and executor it is strictly necessary to have them available to each object (also, we need to check that objects share the same commuicator before having them interact!), but all current functionality of Vector and Matrix (pointwise operations, reductions, SpMV) can be implemented entirely without knowing the partition. You can imagine the read_distributed function doing a global reordering that groups rows/columns from the same block together. The whole point of the interface is for the partition to become unnecessary once you have everything set up. I would assume that applications also tend to work on local indices rather than global indices? For more advanced usage of the halo, maybe we should consider adding a matrix type that takes care of overlaps similar to the P/P^T matrix in libCEED's L-Vector and T-Vector?

[pratikvn]

While I agree that minimizing state stored by a class is a good idea, I think it is important to have enough information to enable the class to be independent, without having to rely on other classes to provide valid information on how the data is being stored in this class. To me, a partition is crucial to a distributed object and provides information on the object's distribution. Regarding the issues you raise:

1. If copying across executors does not copy all the data of that class, then we need to either try to make sure we overload the copy_from to handle that, or disable copy_from and use the = operators with the correct behaviour.
2. Similar to the communicator, we technically do need to verify that the partition of two interacting objects are the same or atleast commute with each other. Example for dot products between two vectors, which can be stored in different ways on different ranks. Right now, we just assume they have the same local sizes. Without a partition object, you disable dot products between vectors even if they have the same global size.
3. The T-vector/L-vector is a good point, but that is also something the partition/index set provides, global to local and local to global maps.

[upsj]

I tend to disagree with the premise that a vector needs to contain information that allows you to understand what its row indices mean. Take Dense for example, without context (e.g. a FEM mesh), you don't know how to interpret the different entries. That context (mapping mesh entries to DOFs) is provided by the application and not stored in the vector. I see the distributed vector in the same way: The application provides a way to interpret what the local indices mean based on its FEM mesh and partition.

Also a side note we haven't talked about: For the most efficient way to store a partition (i.e. one contiguous range for each rank), the local size is entirely sufficient to check for consistency (assuming the default MPI configuration to quit if one rank aborts) Instead of checking the global property "matching partition" everywhere, you check the local property "matching size" everywhere, which is equivalent, because in this case, the read_distributed function doesn't do any (implicit) global reordering.

Finally, I may not have made the fine distinction between our current approach and L/T vectors clear enough:

• Our approach partitions the matrix by row only, storing the off-diagonal block separately and communicating only via a sparse gather before applying the off-diagonal matrix.
• The CEED approach symmetrically stores a slightly larger matrix (local row + overlap) in that it extends both the rows and columns by its set of off-diagonal columns (that of course only works for pattern-symmetric matrices) and communicates both via a sparse gather before the (in this case) single SpMV and via a sparse scatter-add afterwards (to collect contributions from neighboring elements on the boundary). That means more communication, but also makes the local handling much nicer, because you basically get a local domain made out of complete elements (I am not 100% sure about the details of it, but I guess with your experience in RAS and domain composition, this might sound understandable/familiar?)
  The main difference is that there is no longer a separation between diagonal and off-diagonal matrix, and you can imagine the P and P^T matrices as communication operators that do the gather/scatter operations.

EDIT: As a final thought, generally, it is always easier to add something later (store the partition after read_distributed) than it is to remove it later, since the latter involves an interface break. Because IMO the repartitioning is a really advanced use case, I would like to see that handled separately at least, because there may be smoother ways to express it?

EDIT2: I just remembered another issue: If we store the partition in Vector and Matrix, we also need to propagate it to all intermediate vectors in solvers, as well as distributed preconditioner wrappers. That functionality might impact a lot of classes with distributed-specific code.

EDIT3: After some discussions with Natalie, I believe the CEED case is pretty FEM-specific, since each nonzero involving the boundary of an element is computed from the sum of multiple submatrices involving the neighboring elements (see "E vector"), and such a decomposition into individual submatrices may not exist/be well-defined in the general case?

[ginkgo-bot]

Note: This PR changes the Ginkgo ABI:

Functions changes summary: 0 Removed, 0 Changed (40 filtered out), 48 Added functions
Variables changes summary: 0 Removed, 0 Changed, 0 Added variable

For details check the full ABI diff under Artifacts here

[sonarqubecloud]

Kudos, SonarCloud Quality Gate passed!   

0 Bugs
0 Vulnerabilities
0 Security Hotspots
50 Code Smells

87.9% Coverage
4.1% Duplication