Constructor with pairs and size no longer works about bandedmatrices.jl HOT 6 CLOSED

I was using the following snippet:

using BandedMatrices, SparseArrays, LinearAlgebra, BenchmarkTools
n = 10000
h = 1/n
A = BandedMatrix{Float64}(undef, (n,n), (1,1))
A[band(0)] .= -2/h^2
A[band(1)] .= A[band(-1)] .= 1/h^2
B = sparse(A)
C = SymTridiagonal(A)
f = fill(1.0,n)

Then

julia> @btime $A\$f;
  643.493 μs (7 allocations: 469.03 KiB)

julia> @btime $B\$f;
  2.662 ms (64 allocations: 5.42 MiB)

julia> @btime $C\$f;
  173.811 μs (8 allocations: 234.66 KiB)

(I'm on Linux and have weened myself off MKL since they changed their free license to prohibit academic work -- not sure if bundling with JuliaPro changes that. But good to know! In any case, things become more interesting with 2D operators using kron and with BlockBandedMatrices for optimality systems for PDE-constrained optimization. Also, unlike BandedMatrix, SymTridiagonalMatrix has no lu (or rather, cholesky) method -- but then, neither does SymBandedMatrix, although it should, or isn't it exported?)

from bandedmatrices.jl.

Sorry, I’ll fix that now. The easiest is

A = BandedMatrix{Float64}(undef, (n,n), (1,1))
A[band(0)] .= -2/h^2
A[band(1)] .= A[band(-1)] .= 1/h^2

See also

https://github.com/JuliaMatrices/BandedMatrices.jl/blob/master/examples/clarrays.jl

from bandedmatrices.jl.

Ooh, that's clever -- thanks! (I didn't think to look into that example since it seemed too advanced. Maybe your snippet could be added to the documentation, since the discrete gradient or Laplacian are poster children for banded matrices and using a BandedMatrix rather than a SparseArray for the latter gives about a speedup factor of 4?)

from bandedmatrices.jl.

Yes I should add a finite difference example. If you make one I’d be happy to accept a PR.

Note that OpenBLAS (Julia’s default) is pretty slow for banded matrices. If you use MKL (JuliaPro) or Apples BLAS it’s much faster.

from bandedmatrices.jl.

Note you can also use SymTridiagonal in Base.

from bandedmatrices.jl.

The README constructor now works.

from bandedmatrices.jl.