How to make custom LinOp convertible to CSR?

[lahwaacz]

Your ACM paper has an example showing how to define a custom LinOp named StencilMatrix. Later it says that it can be used in the example in Listing 1 instead of the CSR format, but the preconditioner has to be removed, because the StencilMatrix is not convertible to CSR.

Several questions:

• Other matrix formats in Ginkgo are convertible to CSR. Does it mean that if you have a matrix in a different format, let's say Ell, the preconditioner needs to always create a full copy to access the elements? That would be rather inefficient...
• How do you make a custom LinOp convertible to CSR?
• What if the custom LinOp is intended to wrap an external data structure which is essentially CSR already? Is it possible to avoid the conversion/copy?
• Do you have a formal specification of which solvers/preconditioners have which requirements on the LinOp?

[MarcelKoch]

Hi @lahwaacz, I will try to address the points you have raised.

1. Currently, our preconditioners expect that the matrix is in CSR format. If that is not the case, we try to convert the matrix into CSR format, which requires a copy. One exception is the Jacobi preconditioner with block_size=1 only reqiures that the diagonal can be extracted. We should improve our documentation to make this clear.
2. If you want to implement a new LinOp type that is convertible to CSR, you have to implement the ConvertibleTo<Csr<value_type, index_type>> interface. This requires implementing the two functions void convert_to(Csr<value_type, index_type>* result) const and void move_to(Csr<value_type, index_type>* result). Usually the move_to can be implemented using the convert_to, but how you implement the convert_to depends on your custom type. In the StencilMatrix example you've mentioned, this can be implemented with this approach:

matrix_data<...> md(this->get_size());
for(int row = 0; row < this->get_size()[0]; ++row){
  // need to handle boundary correctly!
  md.push_back({row, row - 1, coeffs[0]});
  md.push_back({row, row , coeffs[1]});
  md.push_back({row, row + 1, coeffs[2]});
}
result->read(md);

3. If your custom LinOp wraps your internal data, which is already in CSR format, you don't actually need that new format. If you internally already have the row pointer, column indices, and values available, you can just wrap each of them with our gko::array as a view, and use our CSR format directly. As an example:

gko::dim<2> size(...);
auto mtx = gko::matrix::Csr<>::create(exec, size, 
                                      gko::make_array_view(exec, num_rows + 1, your_row_ptrs),
                                      gko::make_array_view(exec, num_elements, your_col_idx),
                                      gko::make_array_view(exec, num_elements, your_values));

4. Sadly, we don't have a specification of the acceptable input formats for our solvers/preconditioners in our documentation. I've mentioned the preconditioners above, so I will focus now on the solvers. Generally, our solvers only require that the format implements the LinOp interface, i.e. the apply function must be available. The only exception I can think of is the Bicg solver, which additionaly requires that the format implements theTransposable interface.

I hope this could clarify your questions. We will work on making our documentation more expressive with regard to the points you brought up!

[lahwaacz]

Thanks, that clarifies a lot. Looking forward to your improved documentation.

[lahwaacz]

To make MyLinOp convertible to Csr, I also had to add using gko::EnableLinOp< MyLinOp >::convert_to; using gko::EnableLinOp< MyLinOp >::move_to; to avoid warnings due to overloaded virtual function in clang.

[upsj]

Hi Jakub, I'll try to answer the individual questions:

• Yes, since Csr is the only matrix format that can compactly represent arbitrary sparsity patterns with efficient row lookup, it's the format that every preconditioner is based on.
• You need to do two things to make a LinOp convertible to Csr:
  1. derive from ConvertibleTo<Csr<ValueType, IndexType>> and
  2. override/implement the convert_to and move_to functions.
• Yes, other matrix formats also sometimes profit from similarities in storage formats. Take the Csr <--> Coo conversion for example. Both use the same col_idxs and values arrays, just the rows are represented as offsets in Csr and a flat vector in Coo. move_to can thus avoid copies of these two arrays.
• No, we don't have such a specification at the moment, but mainly because all iterative solvers work with arbitrary types, and all preconditioners and direct solvers require Csr inputs.