search

GeoStat-Framework / pentapy

A Python toolbox for pentadiagonal linear systems
https://pentapy.readthedocs.io/
MIT License
14 stars 4 forks source link

Flattened vs Banded #13

Closed joglekara closed 4 years ago

joglekara commented 4 years ago

Very nice work, this is super useful!

I have a quick question, more of a clarification. The documentation says

# create a flattened pentadiagonal matrix
M_flat = (np.random.random((5, size)) - 0.5) * 1e-5

Here, I interpret M_flat to represent each of the 5 diagonals of a pentadiagonal matrix. I believe this is similar to a banded structure.

I see that pentapy supports the matrix being input as Full and Flat. The Full portion makes sense to me, this numpy array has shape (N_col, M_row). I'm a bit confused about Flat because flat typically refers to an (N, N) -> (N^2, 1) transformation (see https://numpy.org/doc/stable/reference/generated/numpy.ndarray.flatten.html). However, the above leads me to believe that Flat is actually just the 5 diagonals, i.e. a numpy array of shape (5, N).

So, this is a long winded way of asking, is it safe to assume that pentapy supports solving pentadiagonal systems represented in a sparse format using just the 5 diagonals to represent the operator (and the right hand side)?

MuellerSeb commented 4 years ago

Hey there! You are totally right. The term flat is ill-chosen at this point. What is meant, is the packed version of the matrix containing only the 5 diagonals defining the given penta-diagonal matrix. So this packed matrix has the shape (5, N). This is described in detail in the documentation, along with the answer of your question: https://geostat-framework.readthedocs.io/projects/pentapy/en/stable/core.html#pentapy.core.solve

But the issue with the term flat would have been a good one during the review for the paper ;-)