ENH: diagonal banded form in solve banded #11344
Conversation
| raise ValueError("Matrix must be square (has shape %s)" % (a.shape,)) | ||
| (nlower, nupper) = l_and_u | ||
| if nlower > n or nupper > n: | ||
| raise ValueError("Number of nonzero diagonals must be less than square dimension") |
There was a problem hiding this comment.
Neither this exception nor the one above are covered by unit tests are the moment.
There was a problem hiding this comment.
Unit tests updated :)
scipy/linalg/basic.py
Outdated
|
|
||
| for i in range(nlower): | ||
| for j in range(n - i - 1): | ||
| lower[i, j] = a[i + j + 1, j] |
There was a problem hiding this comment.
Maybe just a longer-term thought, but if the linalg experts like this addition, I would suggest there's some opportunity for improved efficiency given three separate sets of Python loops, two of them nested. Maybe Cythonization at some point.
scipy/linalg/basic.py
Outdated
| for j in range(n - i - 1): | ||
| lower[i, j] = a[i + j + 1, j] | ||
|
|
||
| return np.concatenate((upper, mid, lower)) |
There was a problem hiding this comment.
There's probably opportunity for efficiency improvement here too if you already know the shape of the final result, but that could be a future change since this is disabled by default in any case.
There was a problem hiding this comment.
I rewrote code, and there is only one matrix now. Concatenation has been removed.
scipy/linalg/basic.py
Outdated
| """ | ||
| Solve the equation a x = b for x, assuming a is banded matrix. | ||
|
|
||
| The matrix a is stored in `ab` using the matrix diagonal ordered form:: | ||
| Until 'diagonal_form' is disabled, the matrix a is stored in `ab` using the matrix diagonal ordered form:: |
There was a problem hiding this comment.
Maybe Until -> Unless could be slightly clearer? Was just thinking the former might imply mutability during the algorithm.
Reference issue
Closes #8362
What does this implement/fix?
It is helpful to have the possibility of passing a banded matrix in its original form (not only in diagonal ordered form as it is now).
Additional information
I needed this feature and I found the above issue, so I thought that it could be included in Scipy. There is an implementation given in #8362, but I've used my own implementation because it doesn't consume so much memory.
Example (matrix from solve_banded docs):