ginkgo-project / ginkgo

Numerical linear algebra software package
https://ginkgo-project.github.io/
BSD 3-Clause "New" or "Revised" License
401 stars 88 forks source link

Issue with batch_ell matrix when solving with non zero per line > 1 #1552

Closed Geoflow closed 7 months ago

Geoflow commented 7 months ago

Hello,

I am actually using Ginkgo in order to solve a batch of linear systems. I have some issues on the execution of the solver (batch BiCgStab ).

In fact the building of the whole structure is Ok. If i take a diagonal matrix(ie nonzero per raw =1 ) the results are good. If i take number of non zeros components per line > 1, the obtained result is false..

My question is more about the storage format. What i understood from test files for batch ell, is that it loops over lines and stores columns index. Is this assumption ok ? In the other hand, i have read in the code that it uses a column major storage (which is supposed to store line indexes )

Thank you for help

upsj commented 7 months ago

I think you might be misunderstanding the ELL format slightly. All sparse matrix formats in Ginkgo are row-based (i.e. they store all column indices for each sparse row), but the ELL format itself internally arranges these entries in a column-major order for better vectorization. To be fair, row-major and column-major are less useful concepts in sparse matrices than they are in dense matrices, which is where some of the confusion might be coming from. If you want to assign the column index or value to the nz th nonzero of row row, you need to use the storage location nz * stride + row where stride usually defaults to num_rows.

Geoflow commented 7 months ago

Thank you for the answer, For the stride i took the number of non zero per line instead of the matrix size eg: ell_index= row*non_zero_per_row+ j (it is usual when you want to "linearize" 2d array )

where j refers to the j-th non zero component of the line.

upsj commented 7 months ago

That would be the incorrect row-major storage layout, yes, the stride needs to be multiplied with the nonzero index inside the row. Only with non_zero_per_row == 1 are the two storage layouts identical

Geoflow commented 7 months ago

Thanks a lot for your advices, it is more clear and it works for me !

upsj commented 7 months ago

That's great to hear, good luck with your experiments ;)