I agree that this seems more like a scipy feature than a numpy feature. Users with structured matrices often use a sparse matrix format, though the API for using them in solvers could use some work. (I have a work-in-progress PR along those lines here: https://github.com/scipy/scipy/pull/6331) Perhaps this polyalgorithm approach could be used to dispatch sparse matrices to the appropriate solver, while optionally checking dense matrices for structure before dispatching them as well. Usage might look like: # if A is sparse, use scipy.sparse.linalg.solve, otherwise use scipy.linalg.solve scipy.linalg.generic_solve(A, b) # converts A to banded representation and calls scipy.linalg.solveh_banded, regardless if A is sparse or dense scipy.linalg.generic_solve(A, b, symmetric=True, banded=(-5, 5)) # runs possibly-expensive checks, then dispatches to the appropriate solver scipy.linalg.generic_solve(A, b, detect_structure=True) (I'm not advocating for "generic_solve" as the final name, I just needed a placeholder.) On Tue, Jan 10, 2017 at 4:58 AM, Ilhan Polat <ilhanpolat@gmail.com> wrote:
I've done some benchmarking and it seems that the packed storage comes with a runtime penalty which agrees with a few links I've found online
https://blog.debroglie.net/2013/09/01/lapack-and-packed-storage/ http://stackoverflow.com/questions/8941678/lapack-are- operations-on-packed-storage-matrices-faster
The access of individual elements in packed stored matrices is expected to be more costly than in full storage, because of the more complicated indexing necessary. Hence, I am not sure if this justifies the absence just by having a dedicated solver for a prescribed structure.
Existence of these polyalgorithms in matlab and not having in lapack should not imply FORTRAN users always know the structure in their matrices. I will also ask in LAPACK message board about this for some context.
But thanks tough. As usual there is more to it than meets the eye probably, ilhan
On Tue, Jan 10, 2017 at 4:16 AM, Robert Kern <robert.kern@gmail.com> wrote:
On Mon, Jan 9, 2017 at 7:10 PM, Ilhan Polat <ilhanpolat@gmail.com> wrote:
Yes, that's precisely the case but when we know the structure we can
just choose the appropriate solver anyhow with a little bit of overhead. What I mean is that, to my knowledge, FORTRAN routines for checking for triangularness etc. are absent.
I'm responding to that. The reason that they don't have those FORTRAN routines for testing for structure inside of a generic dense matrix is that in FORTRAN it's more natural (and efficient) to just use the explicit packed structure and associated routines instead. You would only use a generic dense matrix if you know that there isn't structure in the matrix. So there are no routines for detecting that structure in generic dense matrices.
-- Robert Kern
_______________________________________________ NumPy-Discussion mailing list NumPy-Discussion@scipy.org https://mail.scipy.org/mailman/listinfo/numpy-discussion
_______________________________________________ NumPy-Discussion mailing list NumPy-Discussion@scipy.org https://mail.scipy.org/mailman/listinfo/numpy-discussion