A good way to recycle Gingko data structures at high speeds.

[QuasarJSPark]

Hello, I would like to apply an Iterative method using Ginkgo to the Matrix solving routine of the Newton Solver in the numerical simulator.

The simulator manages variables for matrix solving with a custom data type that is a combination of unordered map and vector to handle Jacobian matrices, and std::vector to handle residue vectors.

From the example files I've seen, it seems that gingko basically holds the data of a matrix or vector in gko::matrix_data<> and then reads it into a Matrix created in shared_ptr format. In the example code, the code to successively generate elements in the Diagonal term I created and put them into the Csr Matrix is as follows.

 using SpMtxType = gko::matrix::Csr<double, int>;

 const int nozeros = 10000;
 const auto ompexec = gko::OmpExecutor::create();
 auto SpMatrix = gko::share(SpMtxType::create(ompexec));
 gko::matrix_data<> mtx_data{gko::dim<2>(nonzeros, nonzeros)};

    for (auto i = 0; i < nonzeros; ++i) {
        mtx_data.nonzeros.emplace_back(i, i, i+2);
    }

    SpMatrix->read(mtx_data);

However, it seems like it would be overkill to recreate the datatype and form a Matrix for every Newton Iteration. In my old Trilinos implementation, I was able to create a CrsGraph and reuse the datatype specific to the solver via replace{global, local}value function.

So my question is this.

1. how can I implement the fastest way to convert from STL type of datastructure{multiple vectors for COO, vector) to Ginkgo's Matrix type?

2. Is there a way to reuse the Matrix type of databases once used by just changing the values? I haven't done a performance test yet, but I think it would be more efficient than allocating new memory each time.

[upsj]

gko::matrix_data is mostly suitable for IO and applications where the matrix assembly is not a costly operation. That doesn't seem to be the case in your application, so I want to suggest either updating the values in-place inside the generated matrix (the ordering of values stays intact when moving from a row-major sorted matrix_data object to a Csr matrix). We have internal functionality to find the nonzero location at a specific column, but we don't expose it yet.
For the initial matrix assembly, you can use device_matrix_data, which is stored in structure-of-arrays (SoA) storage format, which is more suitable for high-performance applications than the array-of-structs (AoS) used in matrix_data

There is some additional complexity involved if you want to use a preconditioner with your problem, because the preconditioner also needs to be recomputed with updated values, but we are also working on a way to provide a clean interface for that.

Finally, if your code is really performance-critical, I want to suggest replacing std::unordered_map (which has poor performance due to design constraints) by a better implementation like the "Swiss tables" from https://github.com/abseil/abseil-cpp

[QuasarJSPark]

Additionally, thanks for the suggestion about Google's Abseil. I'll look into that implementation and the
boost::unordered_flat_map (1.83), the latter performs slightly better, so I'm adopting that implementation.

[upsj]

device_matrix_data is meant for an already fully assembled data, potentially with duplicate elements and zeros (which is what the sum_duplicates and remove_zeros functions are for). Does the sparsity pattern of your matrix stay fixed, or does it change over subsequent nonlinear solves?

[QuasarJSPark]

It is fixed while analyzing the same linear coupled system.

[upsj]

In that case, I would suggest assembling the matrix yourself (e.g. with matrix_assembly_data, which uses unordered_map under the hood), and updating the values in-place in the matrix for subsequent outer iterations.

[QuasarJSPark]

For gko::matrix_assembly_data, set_value() allows you to change the value of an element in a given graph structure, so I think that would be a good idea. I will try that method, thanks for your kind reply.