Use new BandedMatrix \ for large bandwidth solves

[dlfivefifty]

I think the following timing suggests that we would perform significantly better using truncation and banded matrix :

#HelmholtzDirichlet.jl example
@time ∂u∂n = ⨍[G]\uiΓ   # 6.882111 seconds (20.04 M allocations: 4.365 GB, 10.21% gc time)
cfs=pad(coefficients(uiΓ,rangespace(⨍[G])),100)
@time ⨍[G][1:100,1:100]\cfs  # 0.786749 seconds (2.67 M allocations: 525.196 MB, 10.94% gc time)

It's based on LU, so there may be inaccuracies introduced, but in this case they match to 8.9e-12

[dlfivefifty]

For k=0.1 the two are tied! (Probably the cost is dominated by constructing the operator)

@time ∂u∂n = ⨍[G]\uiΓ    # 0.001622 seconds (9.53 k allocations: 871.625 KB)
@time u=⨍[G][1:10,1:10]\cfs  # 0.001138 seconds (4.63 k allocations: 306.406 KB)

[MikaelSlevinsky]

It could even use a Cholesky factorization for (formally) positive definite fundamental solutions.

[MikaelSlevinsky]

Cholesky or LDLT could work provided:

• LowRankFun(...;method=:Cholesky) returns a new type such as CholeskyFun or SymmetricLowRankFun
• getindex(::Operator,::CholeskyFun) returns a new type such as SPDBandedOperator or getindex(::Operator,::SymmetricLowRankFun) returns a new type such as SymmetricBandedOperator
• Populating a SPDBandedOperator results in a SPDBandedMatrix or populating a SymmetricBandedOperator results in a SymmetricBandedMatrix
• we add LL' for SPDBandedMatrix or LDL' for SymmetricBandedMatrix

[dlfivefifty]

I think the new qrfact mitigates the benefits of doing this automatically.