I just posted my first blog post, on using BandedMatrices.jl to solve ordinary differential equations via finite differences / spectral methods: Banded matrices and ordinary differential equations
This is first of a series outlining the results of the NumFocus small development grant to add support for GPUs and other backends to BlockBandedMatrices.jl
I accidentally posted an old working draft. Iāve uploaded the correct version which has some minor changes.
This is a very nice article and introduction to BandedMatrices.jl, thanks! It is disheartening to hear about the OpenBlas issues with banded solvers since installing MKL with Julia still requires some extra workā¦
A few very minor comments/typos in case you are updating it at some point:
.+ 1E-40 in a comment but actually uses 1E-100.btime the lu(L) and lu(LĢ) I assume the comment for the latter meant pivoted LU for sparse matrices?Me too on the boxes and scrolling.
It would be nice if you could make the width of your text and code boxes larger. Right now many of the code boxes require horizontal scrolling to read all the text. (At least on Chrome on my Mac.)
Yes, the formatting still could use improvement. Part of the issue is the Weave.jl generates raw HTML, and only Blogger supports this from what I can tell. Blogger is a bit daunting and I havenāt figured out the ārightā way to adjust the default layout. In particular, I donāt know why the column of text is so narrow.
Your corrections are correct, thanks.
Iāve increased the widths, does this improve things for you?
Very interesting! Do you have an idea of the trade-offs (in terms of bandwith size perhaps) between this approach and that of using iterative solvers with FFTs to do the matvecs?
New widths look much better to me!
Under practically all circumstances iterative methods are not competitive for ODEs, and this is certainly the case for singularly perturbed examples where one would need a very good preconditioner to overcome the ill-conditioning. I donāt know of such a preconditioner.
PDEs are a different story as the complexity of applying an operator is better than inverting an operator. Iāll talk about PDEs in the next post.
Great work. I like the package and your notes. I never thought about obtaining the fourth derivative matrix as a multiplication of two central derivative matrices ā very nice!. Looking forward to the PDE notes.
Iām not sure. On my screen, thereās still more than twice as much horizontal empty space as there is scroll box. I guess if your tools donāt make a āresponsiveā page, thereās not much you can doā¦ I wouldnāt want the text to fill all the space anyway.
Nice post and blog !
I wonder about the opportunity to add fixed size bandwidth (tri/penta/hepta/N-diag) specializations. We implement such N-diag matrices julia types for one of ourās PDE solver and I should compare the performance with this package. We usually use in-place LLT or LDT factorization to solve corresponding symmetric systems.
This sounds like StaticArray v regular array. I donāt think thereās much benefit though in knowing the bandwidths at compile time, as one is usually dealing with large dimensional problems where the actual computation swamps things.
Yes, In our cases which maybe not generally relevant, the inner loop operates through the bandwidth. When it is known statically, the loop is unrolled and this can lead to significant speed-up. This is the case for a cholesky factorization for instance. We define our own macros for this purpose but I guess that it maybe implemented efficiently with StaticArrays.
I forgot to mention that we deal with larges systems but small bandwidths (3-9).
If you are sure knowing that at compile time is useful then itās a good argument. But it would need a new SBandedMatix type that knows the bandwidths in the template variables.
That is exactly what we have implemented.
Please make a PR if you think itās of general use
Hi,
I have tried your package on a 2d problem and it is quite slower than using sparse. My matrix line numbers is Nx * Ny, same for columns. It has 5 bands at 0, Ā±1, Ā±Nx. It seems that you store a big matrix to hold this banded structure and that may be the reason of the slowdown I am experiencing. Do you have any further advice please?
Thank you
I wouldnāt expect a banded matrix structure to do well there since you want something like ādense bandednessā. But @dlfivefifty have you considered this kind of āsparse bandednessā?