How can I update the values of an already assembled sparse matrix?

[edocento]

I cannot find a way to fill a larger matrix (Coo or Csr type) with smaller blocks at specific coordinates in such a way that values added in previous steps are summed up to new values at the same (row,col) coordinate. In other words, I would like to reproduce with your library a function like MatSetValuesLocal(...) in Petsc (with flag ADD_VALUES) .
(Here the docs of the function I am referring to https://petsc.org/main/docs/manualpages/Mat/MatSetValuesLocal.html )

[upsj]

We don't have a good high-performance interface for this yet, but you can start off by using the gko::matrix_assembly_data type and its add_data member function, which is meant for exactly this FEM-like matrix assembly use case. matrix_assembly_data can be passed directly to any matrix->read(...) function of the correct value and index type :) Note that this is meant for assembling a matrix once - if your matrix entries change repeatedly (e.g. nonlinear problem), I'd recommend modifying the entries in-place using the existing column index information in the matrix types.

[edocento]

Computing in advance the sparsity pattern was an option, but I was looking for other ways and "best practice" (maybe more efficient). Anyway, even if it is not high-performance, your answer is a good start for me. I will explore other solutions to my problem and I am interested in knowing how other people has accomplished this task for FEM-like matrices.
P.S. I would mark your answer, but I am new to GitHub and I am not sure if this allows or not new replies

[upsj]

Since you mention CUDA exclusively, a good starting point, in case you don't care about the overhead of storing duplicate entries, might be using Thrust, storing the individual contributions separately, sorting them in row-major order, extracting the sparsity pattern and combining them afterwards.
Knowing the sparsity pattern in advance is extremely useful, since it gives you both global memory requirements (how many nonzeros) as well as local index assignment (which indices belong to which row) in both Csr and Coo. Dealing with dynamic allocation on the GPU can be painful.

Just to break down the components you will probably need:

• device_pointer_cast to mark pointers as being on the device, so Thrust doesn't try to copy their contents first
• make_zip_iterator to combine multiple arrays into a "virtual" array of tuples. That might be necessary to combine row indices and column indices if you store them in a SoA layout.
• sort_by_key to group everything by nonzero location in row-major order
• unique_copy to extract the unique sparsity pattern locations
• reduce_by_key to combine multiple contributions to the same nonzero entry.

Now that I think of it, that might actually be a useful member function for gko::device_matrix_data, I'll think about adding one

[edocento]

That means that a matrix_assembly_data can only be written to sequentially from the host side, since it uses the non-thread-safe std::unordered_map internally. If you want to build everything on the GPU, you can initialize a matrix with a known number of nonzeros and fill it manually, or use device_matrix_data if you don't want to deal with matrix internals directly. Still in the end, you will need a way to check whether a certain nonzero is already present in the partially assembled matrix, and we don't have a nice interface for that yet, as it also involves the structure of the specific matrix formats (Honestly, we haven't had time to investigate how this should look for best performance and generality - probably a small lookup hashmap/compressed bitmap for each row if you are dealing with Csr or a Coo/Csr hybrid). Is it possible for you to compute the sparsity pattern first and fill it with values later? Maybe @nbeams has an idea how this works in MFEM?

Is this related to the fact that, calling write(std::cout, lend(stencil_matrix)) after stencil_matrix->read(data) in the heat_equation example (https://ginkgo-project.github.io/ginkgo/doc/develop/heat_equation.html) outputs a zero matrix? I used a similar filling technique in my project (with matrix_assembly_data in place of matrix_data) and it compiles, but I cannot be sure that the matrix is correctly filled

[upsj]

No, gko::write by default writes out the data in dense format (actually another good point we could probably address, thanks!), so you don't see the nonzeros among all the zeros. Setting the optional layout type parameter to gko::layout_type::coordinate changes this to sparse output.

[Slaedr]

A quick note with regard to Thrust: it should work with OpenMP and ROCm (AMD) as well, so it's actually a reasonable choice.