Eigenfactorization of BandedMatrix with underlying FillArrays storage fails due to strides

[jagot]

I know this is not a problem exclusive to BandedMatrices.jl, but maybe we can figure out where/how to fix it.

(tmp) pkg> st
Status `/private/tmp/Project.toml`
  [aae01518] BandedMatrices v0.15.21
  [1a297f60] FillArrays v0.9.6

julia> using LinearAlgebra, BandedMatrices, FillArrays
[ Info: Precompiling BandedMatrices [aae01518-5342-5314-be14-df237901396f]

julia> A = BandedMatrices._BandedMatrix(Ones{Float64}(3,15), Base.OneTo(15), 1, 1)
15×15 BandedMatrix{Float64,Ones{Float64,2,Tuple{Base.OneTo{Int64},Base.OneTo{Int64}}},Base.OneTo{Int64}}:
 1.0  1.0   ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅
 1.0  1.0  1.0   ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅
  ⋅   1.0  1.0  1.0   ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅
  ⋅    ⋅   1.0  1.0  1.0   ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅
  ⋅    ⋅    ⋅   1.0  1.0  1.0   ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅
  ⋅    ⋅    ⋅    ⋅   1.0  1.0  1.0   ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅
  ⋅    ⋅    ⋅    ⋅    ⋅   1.0  1.0  1.0   ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅
  ⋅    ⋅    ⋅    ⋅    ⋅    ⋅   1.0  1.0  1.0   ⋅    ⋅    ⋅    ⋅    ⋅    ⋅
  ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅   1.0  1.0  1.0   ⋅    ⋅    ⋅    ⋅    ⋅
  ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅   1.0  1.0  1.0   ⋅    ⋅    ⋅    ⋅
  ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅   1.0  1.0  1.0   ⋅    ⋅    ⋅
  ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅   1.0  1.0  1.0   ⋅    ⋅
  ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅   1.0  1.0  1.0   ⋅
  ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅   1.0  1.0  1.0
  ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅   1.0  1.0

julia> eigen(Symmetric(A))
ERROR: MethodError: no method matching strides(::Ones{Float64,2,Tuple{Base.OneTo{Int64},Base.OneTo{Int64}}})
Closest candidates are:
  strides(::SubArray) at subarray.jl:331
  strides(::Union{Base.ReinterpretArray{T,N,S,A} where S where A<:Union{SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray} where N where T, Base.ReshapedArray{T,N,A,MI} where MI<:Tuple{Vararg{Base.MultiplicativeInverses.SignedMultiplicativeInverse{Int64},N} where N} where A<:Union{Base.ReinterpretArray{T,N,S,A} where S where A<:Union{SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray} where N where T, SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray} where N where T, DenseArray}) at reinterpretarray.jl:76
  strides(::PermutedDimsArray{T,N,perm,iperm,AA} where AA<:AbstractArray where iperm) where {T, N, perm} at permuteddimsarray.jl:63
  ...
Stacktrace:
 [1] strides(::SubArray{Float64,2,Ones{Float64,2,Tuple{Base.OneTo{Int64},Base.OneTo{Int64}}},Tuple{UnitRange{Int64},Base.Slice{Base.OneTo{Int64}}},false}) at ./subarray.jl:331
 [2] stride(::SubArray{Float64,2,Ones{Float64,2,Tuple{Base.OneTo{Int64},Base.OneTo{Int64}}},Tuple{UnitRange{Int64},Base.Slice{Base.OneTo{Int64}}},false}, ::Int64) at ./subarray.jl:339
 [3] stride1(::SubArray{Float64,2,Ones{Float64,2,Tuple{Base.OneTo{Int64},Base.OneTo{Int64}}},Tuple{UnitRange{Int64},Base.Slice{Base.OneTo{Int64}}},false}) at /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.5/LinearAlgebra/src/LinearAlgebra.jl:197
 [4] _chkstride1 at /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.5/LinearAlgebra/src/LinearAlgebra.jl:203 [inlined]
 [5] chkstride1 at /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.5/LinearAlgebra/src/LinearAlgebra.jl:201 [inlined]
 [6] sbtrd!(::Char, ::Char, ::Int64, ::Int64, ::SubArray{Float64,2,Ones{Float64,2,Tuple{Base.OneTo{Int64},Base.OneTo{Int64}}},Tuple{UnitRange{Int64},Base.Slice{Base.OneTo{Int64}}},false}, ::Array{Float64,1}, ::Array{Float64,1}, ::Array{LinearAlgebra.Givens{Float64},1}, ::Array{Float64,1}) at /Users/jagot/.julia/packages/BandedMatrices/i7LFa/src/symbanded/bandedeigen.jl:193
 [7] eigen!(::Symmetric{Float64,BandedMatrix{Float64,Ones{Float64,2,Tuple{Base.OneTo{Int64},Base.OneTo{Int64}}},Base.OneTo{Int64}}}) at /Users/jagot/.julia/packages/BandedMatrices/i7LFa/src/symbanded/bandedeigen.jl:37
 [8] eigen(::Symmetric{Float64,BandedMatrix{Float64,Ones{Float64,2,Tuple{Base.OneTo{Int64},Base.OneTo{Int64}}},Base.OneTo{Int64}}}) at /Users/jagot/.julia/packages/BandedMatrices/i7LFa/src/symbanded/bandedeigen.jl:26
 [9] top-level scope at REPL[13]:1

I guess one could add a dispatch to eigen that automatically converts A.data to a dense array before proceeding?

[jagot]

Something like this:

function eigen(A::Symmetric{T,<:BandedMatrix{T,<:AbstractFill{T}}}) where T <: Real
    Ap = parent(A)
    B = BandedMatrices._BandedMatrix(Matrix(Ap.data), Ap.raxis, Ap.l, Ap.u)
    eigen(Symmetric(B))
end

[dlfivefifty]

I believe the "logic" for not including a fallback in Base is that it will crash Julia with a large sparse matrix, which isn't a worry we have here so it's fine.

Why not just make that the default for <:AbstractMatrix? <:AbstractFill seems too specialised.

Though.... it seems the wrong access point. eigen(...) surely calls eigen!(...) and so weakening the type restructions there might be better

[jagot]

Why not just make that the default for <:AbstractMatrix? <:AbstractFill seems too specialised.

You mean change the default implementation to

# Default implementation
function eigen!(A::Symmetric{T,<:BandedMatrix{T,<:StridedMatrix{T}}) where T <: Real
...
end

?

Though.... it seems the wrong access point. eigen(...) surely calls eigen!(...) and so weakening the type restructions there might be better

Well, I'd assume the mutating version eigen! to not allocate, and materializing e.g. a FillArray requires the allocation of a dense matrix. I'm fine with eigen allocating though.

[dlfivefifty]

No change wherever eigen(...) is overloaded to work for all <:AbstractMatrix, but still call eigen!