<div dir="ltr">I agree that this seems more like a scipy feature than a numpy feature.<div><br></div><div>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: <a href="https://github.com/scipy/scipy/pull/6331">https://github.com/scipy/scipy/pull/6331</a>)</div><div><br></div><div>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:</div><div><br></div><div># if A is sparse, use scipy.sparse.linalg.solve, otherwise use scipy.linalg.solve</div>scipy.linalg.generic_solve(A, b)<div><br></div><div># converts A to banded representation and calls scipy.linalg.solveh_banded, regardless if A is sparse or dense<br><div>scipy.linalg.generic_solve(A, b, symmetric=True, banded=(-5, 5))</div><div><br></div><div># runs possibly-expensive checks, then dispatches to the appropriate solver</div>scipy.linalg.generic_solve(A, b, detect_structure=True)</div><div><br></div><div><br></div><div>(I'm not advocating for "generic_solve" as the final name, I just needed a placeholder.)</div></div><div class="gmail_extra"><br><div class="gmail_quote">On Tue, Jan 10, 2017 at 4:58 AM, Ilhan Polat <span dir="ltr"><<a href="mailto:ilhanpolat@gmail.com" target="_blank">ilhanpolat@gmail.com</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div dir="ltr"><div><div><div>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 <br><br><a href="https://blog.debroglie.net/2013/09/01/lapack-and-packed-storage/" target="_blank">https://blog.debroglie.net/<wbr>2013/09/01/lapack-and-packed-<wbr>storage/</a><br><a href="http://stackoverflow.com/questions/8941678/lapack-are-operations-on-packed-storage-matrices-faster" target="_blank">http://stackoverflow.com/<wbr>questions/8941678/lapack-are-<wbr>operations-on-packed-storage-<wbr>matrices-faster</a><br><br>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. <br><br></div>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.<br><br></div>But thanks tough. As usual there is more to it than meets the eye probably, <br></div>ilhan <br><div><div><div><br><br> <br></div></div></div></div><div class="gmail_extra"><br><div class="gmail_quote"><div><div class="h5">On Tue, Jan 10, 2017 at 4:16 AM, Robert Kern <span dir="ltr"><<a href="mailto:robert.kern@gmail.com" target="_blank">robert.kern@gmail.com</a>></span> wrote:<br></div></div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div><div class="h5"><div dir="ltr"><span>On Mon, Jan 9, 2017 at 7:10 PM, Ilhan Polat <<a href="mailto:ilhanpolat@gmail.com" target="_blank">ilhanpolat@gmail.com</a>> wrote:<br>><br>> 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. <br><br></span>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.<br> <br>--<br>Robert Kern</div>
<br></div></div><span class="">______________________________<wbr>_________________<br>
NumPy-Discussion mailing list<br>
<a href="mailto:NumPy-Discussion@scipy.org" target="_blank">NumPy-Discussion@scipy.org</a><br>
<a href="https://mail.scipy.org/mailman/listinfo/numpy-discussion" rel="noreferrer" target="_blank">https://mail.scipy.org/mailman<wbr>/listinfo/numpy-discussion</a><br>
<br></span></blockquote></div><br></div>
<br>______________________________<wbr>_________________<br>
NumPy-Discussion mailing list<br>
<a href="mailto:NumPy-Discussion@scipy.org">NumPy-Discussion@scipy.org</a><br>
<a href="https://mail.scipy.org/mailman/listinfo/numpy-discussion" rel="noreferrer" target="_blank">https://mail.scipy.org/<wbr>mailman/listinfo/numpy-<wbr>discussion</a><br>
<br></blockquote></div><br></div>