Fill device_matrix_data

[MarcelKoch]

This PR allows to fill the entries of a device_matrix_data with a specified value. The main use case is in combination with sum_duplicates and remove_zeros to simplify the assembly setup.

[upsj]

Is there any purpose in filling with something else than zeroes? I think I would prefer a more specific fill_zero that fills it with (0, 0, 0) entries.

[MarcelKoch]

This would also mean to change the constructor (which is more important to me). The alternative would be to introduce an enum for the two different initialization modes (i.e. fill_mode::uninitialized, fill_mode::zero), which would be fine with me. The only downside is that it would be longer, fill_mode::zero vs {0, 0, 0}.

[upsj]

Does it need to happen in a single function call though? Can't we ask the user to call fill_zero instead? fill_mode introduces a fundamentally new concept we haven't used so far, so I'm wondering whether that is the right thing to do

[MarcelKoch]

For me the constructor is more important than the fill method. Users usually know when they create the device matrix data if it's going to be used in an assembly setting or not. For the assembly, I assume that they almost always want to zero out the data, because otherwise sum_duplicates and remove_zeros doesn't work. So we can make it easier for users by allowing them to fill the data during the constructor.

[yhmtsai]

LGTM for the implementation.
another question: I assume users will allocate more than they need in practice, is it the case?

[yhmtsai]

Suggested change

	/**
	/**

[upsj]

Why does it need to be done in a single step? Can't it be two steps (construction + fill)? Can you describe a matrix assembly scenario where the number of values to be written isn't known beforehand? And if so, couldn't it be more efficient to compute the storage requirements beforehand and leave no gaps?

[MarcelKoch]

@upsj @yhmtsai it is pretty common in FEM applications to not use the total number of DOFs, but instead use num_elements * num_dofs_per_element. This is an overestimation, but usually not too excessive. This makes the assembly a lot easier, since you can then assemble each element in parallel, and then have only a post-processing step to create the correct connectivity pattern. AFAIK, this is also the idea behind the sum_duplicates functionality.
If you would use a structure with just size of num_total_dofs, then you would need atomic updates to handle DOFs that are shared by multiple elements, which are nearly all DOFs in FEM with continuous elements. From that perspective, the approach of overallocating and reducing at the end seems more efficient to me for a single matrix assembly. Of course, if you have transient problems, we might want to provide a better way to repeatedly assemble a matrix with the same sparsity structure. But IMO that is out of the scope for this PR.

[MarcelKoch]

@upsj doing it in one step is syntactic sugar. But so is setting the dimension. The one step is more convenient, so it might be easier for the user, both in terms of slightly less code, and discoverbility. Right now, do our users know that they will get uninitialized memory? I would not count on it.
Besides, I think how the memory is initialized is a natural property of the constructor, so it makes sense to specify it there.

[upsj]

I know the common approach of assembling elements independently (which is what I built sum_duplicates for), but what I'm not familiar with is the case where part of the preallocated storage is left unused. Is that related to boundary conditions? It should be straightforward for users to fix this in their own assembly routines (just set the values to 0), but I'm not opposed fundamentally to it, just want to make it an explicit operation.

On the syntactic sugar side, I slightly disagree: There are unchangeable properties of the device_matrix_data (size, nnz) and there are changeable properties (actual values), and I think it is unacceptable to leave unchangeable properties uninitialized, but the changeable ones are fine. We do the same thing across almost all data structures in Ginkgo (only structures with row_ptr semantics like Csr and Sellp get initialized partially, because otherwise an empty matrix would already crash many algorithms), which is why I at least want to think twice before introducing some initialization semantics. If we do want to add them, I would probably consider adding them to other matrix types as well in the long term.

[MarcelKoch]

Used storage can appear in FVM or FDM on the boundary. From the user perspective, it should not matter if they explicitly set those to zero or just skip them. In some cases, setting the values might seem a bit odd, because it would involve invalid non-zero coordinates.
If you think for example of a 3pt stencil, the idiomatic way to assemble it would be

if(i > 0) add(i, i - 1, -1);
add(i, i, 2);
if(i < n - 1) add(i, i + 1, -1);

where add(0, -1, 0) would seem a bit odd (and might even lead to invalid memory access). Similar things could happen on the FVM side with the boundary handling.

[upsj]

The question there is: Is the resulting row/column then zero itself? Because that would make the matrix singular, and maybe bring us back to #842 if we want to use it in solvers. Do we need some post-processing functionality for that use case as well?

[MarcelKoch]

Not in my example, since i only iterates over the interior points. But I think that gets a bit off-topic. Point is, users might not expect the memory to be uninitialized and just assume that it is zeroed out.

[upsj]

I see, that makes sense. But that would not work in a more generic grid, where the indexing isn't as simple. But yes, let's not get too much off track.

What I'm saying is: All of our matrix data structures do not initialize their arrays, so that is an expectation that would already be formed in other parts of the code. But again, I am not opposed to it, I just want to suggest that we do it everywhere, if we do it. In that case, fill_mode should be defined somewhere more central, and array, matrix types and device_matrix_data should all use them (of course not in this PR)

[MarcelKoch]

I agree that if we accept this PR, it should be added throughout Ginkgo. But IMO, this should be done in subsequent PRs, to keep the changes in each one short. I would also leave fill_mode where it is for now, since moving it to array.hpp for example but not adding it to gko::array seems quite incomplete to me.
I would also suggest changing the enum name to init_mode, since that might be a bit clearer.

[upsj]

The take-away from the discussion in today's meeting:

• The safety of pre-initialized data should outweigh the performance advantage of not initializing data, so the default case should be initializing it to zero, with the option to leave it uninitialized with a special constructor/create function.
• Move forward with some kind of initialized constructor, either using an enum to mark initialization status or with a separate function (create_)uninitialized with the same parameters as the constructor. My personal preference would be on different functions, among other things since the enum would imply a runtime choice rather than a static choice in the source code.
• Make all matrix types initialized to zero by default, with a create_uninitialized copy of the original create function
• Arrays should remain uninitialized, since they don't have a canonical default value.