Cholesky broken for diagonal matrices

[dlfivefifty]

julia> cholesky(Symmetric(BandedMatrix(0 => 1:∞)))
Cholesky{Float64, InfiniteLinearAlgebra.AdaptiveCholeskyFactors{Float64, BandedMatrix{Float64, Matrix{Float64}, Base.OneTo{Int64}}, BandedMatrix{Float64, ApplyArray{Float64, 2, typeof(*), Tuple{Matrix{Int64}, ApplyArray{Float64, 2, typeof(vcat), Tuple{ApplyArray{Float64, 2, typeof(hcat), Tuple{Zeros{Float64, 2, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}}, InfiniteArrays.ReshapedArray{Int64, 2, InfiniteArrays.InfUnitRange{Int64}, Tuple{Int64, Infinities.InfiniteCardinal{0}}, Tuple{}}}}}}}}, InfiniteArrays.OneToInf{Int64}}}}
U factor:
ℵ₀×ℵ₀ UpperTriangular{Float64, InfiniteLinearAlgebra.AdaptiveCholeskyFactors{Float64, BandedMatrix{Float64, Matrix{Float64}, Base.OneTo{Int64}}, BandedMatrix{Float64, ApplyArray{Float64, 2, typeof(*), Tuple{Matrix{Int64}, ApplyArray{Float64, 2, typeof(vcat), Tuple{ApplyArray{Float64, 2, typeof(hcat), Tuple{Zeros{Float64, 2, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}}, InfiniteArrays.ReshapedArray{Int64, 2, InfiniteArrays.InfUnitRange{Int64}, Tuple{Int64, Infinities.InfiniteCardinal{0}}, Tuple{}}}}}}}}, InfiniteArrays.OneToInf{Int64}}}} with indices OneToInf()×OneToInf():
Error showing value of type Cholesky{Float64, InfiniteLinearAlgebra.AdaptiveCholeskyFactors{Float64, BandedMatrix{Float64, Matrix{Float64}, Base.OneTo{Int64}}, BandedMatrix{Float64, ApplyArray{Float64, 2, typeof(*), Tuple{Matrix{Int64}, ApplyArray{Float64, 2, typeof(vcat), Tuple{ApplyArray{Float64, 2, typeof(hcat), Tuple{Zeros{Float64, 2, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}}, InfiniteArrays.ReshapedArray{Int64, 2, InfiniteArrays.InfUnitRange{Int64}, Tuple{Int64, Infinities.InfiniteCardinal{0}}, Tuple{}}}}}}}}, InfiniteArrays.OneToInf{Int64}}}}:
ERROR: BoundsError: attempt to access 0-element UnitRange{Int64} at index [1]
Stacktrace:
  [1] throw_boundserror(A::UnitRange{Int64}, I::Int64)
    @ Base ./abstractarray.jl:734
  [2] getindex
    @ ./range.jl:930 [inlined]
  [3] _shift
    @ ~/.julia/packages/BandedMatrices/AzQRm/src/banded/BandedMatrix.jl:884 [inlined]
  [4] similar
    @ ~/.julia/packages/BandedMatrices/AzQRm/src/banded/BandedMatrix.jl:891 [inlined]
  [5] _convert_common_container
    @ ~/.julia/packages/BandedMatrices/AzQRm/src/banded/BandedMatrix.jl:128 [inlined]
  [6] convert(::Type{BandedMatrix{Float64, Matrix{…}, Base.OneTo{…}}}, M::SubArray{Float64, 2, BandedMatrix{Float64, Matrix{…}, Base.OneTo{…}}, Tuple{UnitRange{…}, UnitRange{…}}, false})
    @ BandedMatrices ~/.julia/packages/BandedMatrices/AzQRm/src/banded/BandedMatrix.jl:137
  [7] convert(::Type{BandedMatrix{Float64, Matrix{…}, Base.OneTo{…}}}, M::SubArray{Float64, 2, BandedMatrix{Float64, Matrix{…}, Base.OneTo{…}}, Tuple{UnitRange{…}, UnitRange{…}}, false})
    @ BandedMatrices ~/.julia/packages/BandedMatrices/AzQRm/src/banded/BandedMatrix.jl:147 [inlined]
  [8] materialize!(M::ArrayLayouts.MulAdd{BandedMatrices.BandedRows{…}, BandedMatrices.BandedColumns{…}, BandedMatrices.BandedColumns{…}, Float64, Adjoint{…}, SubArray{…}, SubArray{…}})
    @ BandedMatrices ~/.julia/packages/BandedMatrices/AzQRm/src/generic/matmul.jl:182
  [9] muladd!
    @ ~/.julia/packages/ArrayLayouts/FrPmc/src/muladd.jl:70 [inlined]
 [10] partialcholesky!(F::InfiniteLinearAlgebra.AdaptiveCholeskyFactors{Float64, BandedMatrix{…}, BandedMatrix{…}}, n::Int64)
    @ InfiniteLinearAlgebra ~/.julia/packages/InfiniteLinearAlgebra/TQWtv/src/infcholesky.jl:35
 [11] getindex
    @ InfiniteLinearAlgebra ~/.julia/packages/InfiniteLinearAlgebra/TQWtv/src/infcholesky.jl:42 [inlined]
 [12] getindex(A::UpperTriangular{Float64, InfiniteLinearAlgebra.AdaptiveCholeskyFactors{Float64, BandedMatrix{…}, BandedMatrix{…}}}, i::Int64, j::Int64)
    @ LinearAlgebra ~/Projects/julia-1.10/usr/share/julia/stdlib/v1.10/LinearAlgebra/src/triangular.jl:241
 [13] alignment(io::IOContext{…}, X::AbstractVecOrMat, rows::Vector{…}, cols::Vector{…}, cols_if_complete::Int64, cols_otherwise::Int64, sep::Int64, ncols::Infinities.InfiniteCardinal{…})
    @ Base ./arrayshow.jl:69
 [14] _print_matrix(io::IOContext{…}, X::AbstractVecOrMat, pre::String, sep::String, post::String, hdots::String, vdots::String, ddots::String, hmod::Int64, vmod::Int64, rowsA::InfiniteArrays.InfUnitRange{…}, colsA::InfiniteArrays.InfUnitRange{…})
    @ Base ./arrayshow.jl:207
 [15] print_matrix(io::IOContext{…}, X::UpperTriangular{…}, pre::String, sep::String, post::String, hdots::String, vdots::String, ddots::String, hmod::Int64, vmod::Int64)
    @ Base ./arrayshow.jl:171
 [16] 
    @ Base ./arrayshow.jl:171 [inlined]
 [17] print_array
    @ ./arrayshow.jl:358 [inlined]
 [18] show(io::IOContext{Base.TTY}, ::MIME{Symbol("text/plain")}, X::UpperTriangular{Float64, InfiniteLinearAlgebra.AdaptiveCholeskyFactors{Float64, BandedMatrix{…}, BandedMatrix{…}}})
    @ Base ./arrayshow.jl:399
 [19] show(io::IOContext{Base.TTY}, mime::MIME{Symbol("text/plain")}, C::Cholesky{Float64, InfiniteLinearAlgebra.AdaptiveCholeskyFactors{Float64, BandedMatrix{…}, BandedMatrix{…}}})
    @ LinearAlgebra ~/Projects/julia-1.10/usr/share/julia/stdlib/v1.10/LinearAlgebra/src/cholesky.jl:563
 [20] (::REPL.var"#55#56"{REPL.REPLDisplay{REPL.LineEditREPL}, MIME{Symbol("text/plain")}, Base.RefValue{Any}})(io::Any)
    @ REPL ~/Projects/julia-1.10/usr/share/julia/stdlib/v1.10/REPL/src/REPL.jl:273
 [21] with_repl_linfo(f::Any, repl::REPL.LineEditREPL)
    @ REPL ~/Projects/julia-1.10/usr/share/julia/stdlib/v1.10/REPL/src/REPL.jl:569
 [22] display(d::REPL.REPLDisplay, mime::MIME{Symbol("text/plain")}, x::Any)
    @ REPL ~/Projects/julia-1.10/usr/share/julia/stdlib/v1.10/REPL/src/REPL.jl:259
 [23] display(d::REPL.REPLDisplay, x::Any)
    @ REPL ~/Projects/julia-1.10/usr/share/julia/stdlib/v1.10/REPL/src/REPL.jl:278
 [24] display(x::Any)
    @ Base.Multimedia ./multimedia.jl:340
 [25] #invokelatest#2
    @ ./essentials.jl:887 [inlined]
 [26] invokelatest
    @ ./essentials.jl:884 [inlined]
 [27] print_response(errio::IO, response::Any, show_value::Bool, have_color::Bool, specialdisplay::Union{Nothing, AbstractDisplay})
    @ REPL ~/Projects/julia-1.10/usr/share/julia/stdlib/v1.10/REPL/src/REPL.jl:315
 [28] (::REPL.var"#57#58"{REPL.LineEditREPL, Pair{Any, Bool}, Bool, Bool})(io::Any)
    @ REPL ~/Projects/julia-1.10/usr/share/julia/stdlib/v1.10/REPL/src/REPL.jl:284
 [29] with_repl_linfo(f::Any, repl::REPL.LineEditREPL)
    @ REPL ~/Projects/julia-1.10/usr/share/julia/stdlib/v1.10/REPL/src/REPL.jl:569
 [30] print_response(repl::REPL.AbstractREPL, response::Any, show_value::Bool, have_color::Bool)
    @ REPL ~/Projects/julia-1.10/usr/share/julia/stdlib/v1.10/REPL/src/REPL.jl:282
 [31] (::REPL.var"#do_respond#80"{Bool, Bool, REPL.var"#93#103"{REPL.LineEditREPL, REPL.REPLHistoryProvider}, REPL.LineEditREPL, REPL.LineEdit.Prompt})(s::REPL.LineEdit.MIState, buf::Any, ok::Bool)
    @ REPL ~/Projects/julia-1.10/usr/share/julia/stdlib/v1.10/REPL/src/REPL.jl:911
 [32] (::VSCodeServer.var"#101#104"{REPL.var"#do_respond#80"{Bool, Bool, REPL.var"#93#103"{…}, REPL.LineEditREPL, REPL.LineEdit.Prompt}})(mi::REPL.LineEdit.MIState, buf::IOBuffer, ok::Bool)
    @ VSCodeServer ~/.vscode/extensions/julialang.language-julia-1.66.2/scripts/packages/VSCodeServer/src/repl.jl:122
 [33] #invokelatest#2
    @ Base ./essentials.jl:887 [inlined]
 [34] invokelatest
    @ Base ./essentials.jl:884 [inlined]
 [35] run_interface(terminal::REPL.Terminals.TextTerminal, m::REPL.LineEdit.ModalInterface, s::REPL.LineEdit.MIState)
    @ REPL.LineEdit ~/Projects/julia-1.10/usr/share/julia/stdlib/v1.10/REPL/src/LineEdit.jl:2656
 [36] run_frontend(repl::REPL.LineEditREPL, backend::REPL.REPLBackendRef)
    @ REPL ~/Projects/julia-1.10/usr/share/julia/stdlib/v1.10/REPL/src/REPL.jl:1312
 [37] (::REPL.var"#62#68"{REPL.LineEditREPL, REPL.REPLBackendRef})()
    @ REPL ~/Projects/julia-1.10/usr/share/julia/stdlib/v1.10/REPL/src/REPL.jl:386
Some type information was truncated. Use `show(err)` to see complete types.

[TSGut]

We can fix this by changing partialcholesky! to check whether it hits the diagonal special case:

function partialcholesky!(F::AdaptiveCholeskyFactors{T,<:BandedMatrix}, n::Int) where T
    if n > F.ncols
        _,u = bandwidths(F.data.array)
        resizedata!(F.data,n+u,n+u);
        ncols = F.ncols
        kr = ncols+1:n
        factors = view(F.data.data,kr,kr)
        banded_chol!(factors, UpperTriangular)
        # multiply remaining columns
        kr2 = max(n-u+1,kr[1]):n
        U1 = UpperTriangular(view(F.data.data,kr2,kr2))
        B = view(F.data.data,kr2,n+1:n+u)
        ldiv!(U1',B)
        if u > 0
            muladd!(-one(T),B',B,one(T),view(F.data.data,n+1:n+u,n+1:n+u))
        else # Diagonal special case
            muladd!(-one(T),zeros(T,0,0),zeros(T,0,0),one(T),view(F.data.data,n+1:n+u,n+1:n+u))
        end
        F.ncols = n
    end
    F
end

The only change is the else case of the if statement at the bottom. A more proper fix I suppose would be to fix this upstream in ArrayLayouts. The issue is basically that it's handing an empty banded matrix with bandwidths (0,-1) which breaks MulAdd somewhere. (the placement of the if statement should happen sooner for efficiency and we can avoid muladd altogether in the diagonal case if we special case it but it's more clear what breaks for this comment when written this way)

Basically two options:

1. Is this what you want MulAdd to do when given a degenerate banded matrix in which case I can clean the above up and turn it into a PR here or
2. Is this something to instead fix upstream in ArrayLayouts.jl?