lu(A) is not banded

[Veenty]

I have a matrix with band with (2,2), but when I do lu(A), the L matrix is not banded. It seems its because the subdiagonals are too big so LU prefers to do some scrabling around. Is there a way to avoid this problem (should I even care?)

[dlfivefifty]

Can you give code?

[Veenty]

For context, I'm doing a 1D boundary value problem solver using an integral formulation. It might be super case specific, so I'm giving exactly how I'm generating the matrix.

So L is a standard matrix and U is banded. I have tried doing lu(G, NoPivot()), but that generates 2 standard matrices

using LinearAlgebra
using BandedMatrices

function MultiplicationMatrix(a, n)


    #a are the chebyshev coefficients of the function to multiply by
    #it has already been choped of

    b = length(a)

    if n < b
        error("n must be greater than the number of coefficients")
    end

    M = BandedMatrix{Float64}(zeros(n, n), (b-1, b-1))

    #Toeplitz part
    for i = 1:n

        M[i,i] = a[1]
        for j = 1:b-1
            if i-j > 0 
                M[i-j,i] += a[j+1]/2
            end
        end

        for j = 1:b-1
            if i+j ≤ n
                M[i+j,i] += a[j+1]/2
            end
        end

    end

    #Almost Hankel part

    for j=1:b
        for i = 2:b
            if i+j -1 ≤ b 
                M[i,j] += a[ i+j-1   ]/2
            end
            
        end
    end

    return M


end


function IntegralMatrices(n, scale) 

    if n < 2
        error("n must be greater than 1")
    end

    ∫ = BandedMatrix{Float64}(zeros(n, n), (1, 1))

    ∫[2,1] = 1

    ∫[3,2] = 1/4

    for m =2:n-1

        ∫[m,m+1] = -1/(2*(m-1))
        if m < n-1
            ∫[m+2,m+1] = +1/(2*(m+1))
        end
        # ∫[m+2,m+1] = +1/(2*(m+1))

    end

    ∫∫ = BandedMatrix{Float64}(zeros(n, n), (2, 2))

    ∫∫[3,1] = 1/4

    ∫∫[2,2] = -1/8
    ∫∫[4,2] = 1/(8*3)

    ∫∫[3,3] = - 1/4*(1/6 + 1/2)
    ∫∫[5,3] = 1/(3*4*4)


    for m =3:n-1

        k = m+1
        ∫∫[k-2,k] = 1/(4*(m-1)*(m-2))
        ∫∫[k,k] = -1/(m^2 -1)/2
        if k +2 ≤ n
            ∫∫[k+2,k] = 1/(4*(m+1)*(m+2))
        end
        
    end

    return (scale/2)*∫, (scale/2)^2*∫∫


end

#parameters

k = 5
δ = 0.1
κ = 1.
v = 2π/κ
λ² = 0.
Nr = 20;

ψ_coef = SA[v*(1 - κ*δ/2) , v*(κ*δ/2) ] 
ψ²_coef = SA[ψ_coef[1]^2 + ψ_coef[2]^2 /2, 2*ψ_coef[1]*ψ_coef[2], ψ_coef[2]^2/2]


∫, ∫∫ = IntegralMatrices(Nr+2, δ)

Ma = MultiplicationMatrix(ψ²_coef, Nr+2)
Mb = MultiplicationMatrix(v*κ*ψ_coef, Nr+2)

G = Ma + Mb*∫ - (λ² + (v*κ*k)^2).*∫∫

lu(G)

[dlfivefifty]

Just post the output

[Veenty]

Sure, I changed Nr = 7, so have a smaller output. I think its just not using the lu from BandedMatrices?

BandedMatrices.BandedLU{Float64, BandedMatrix{Float64, Matrix{Float64}, Base.OneTo{Int64}}}
L factor:
7×7 Matrix{Float64}:
  1.0           0.0           0.0          0.0          0.0         0.0    0.0
  0.157459      1.0           0.0          0.0          0.0         0.0    0.0
 -0.0145028     0.0663674     1.0          0.0          0.0         0.0    0.0
 -5.37639e-20  -0.00184158    0.0608003    1.0          0.0         0.0    0.0
  2.30627e-20  -1.48511e-20  -0.000513882  0.0590486    1.0         0.0    0.0
  2.7122e-22    2.7288e-21   -4.88637e-20  1.35344e-19  0.0578368   1.0    0.0
 -4.75262e-23   8.87586e-23  -1.20219e-21  3.43052e-20  0.00025367  0.057  1.0
U factor:
7×7 BandedMatrix{Float64} with bandwidths (2, 4):
 35.728   1.87522  -3.46945e-18   0.0        0.0         ⋅           ⋅ 
  0.0    35.7288    0.937612     -0.296088   0.0        0.0          ⋅ 
  0.0     0.0      36.0112        1.42607   -0.086359   0.0         0.0
   ⋅      0.0       0.0          35.7394     1.56794   -0.0328987   0.0
   ⋅       ⋅        0.0           0.0       35.6649     1.64276    -0.0111033
   ⋅       ⋅         ⋅            0.0        0.0       35.633       1.68834
   ⋅       ⋅         ⋅             ⋅         0.0        0.0        35.6127

[Veenty]

So I went through the source code, and it seems that the LU decomposition is fine, but it seems that is not on how the L factor is computed for display reasons?

[dlfivefifty]

that sounds right a PR would be appreciated

[RobertGregg]

I ran into a problem using ldiv! which I think is related to this issue. Here is a link to discourse thread. Basically, ldiv! is about twice as slow as a sparse linear solve, however, using the backslash operator is faster for BandedMatrices.