Strang matrix factorization - MethodError

[braamvandyk]

Minimal example:

julia> using LinearAlgebra

julia> N = 5
5

julia> A = Tridiagonal(-1 .* ones(N - 1), 2 .* ones(N), -1 .* ones(N-1))
5×5 Tridiagonal{Float64, Vector{Float64}}:
  2.0  -1.0    ⋅     ⋅     ⋅
 -1.0   2.0  -1.0    ⋅     ⋅
   ⋅   -1.0   2.0  -1.0    ⋅
   ⋅     ⋅   -1.0   2.0  -1.0
   ⋅     ⋅     ⋅   -1.0   2.0

julia> factorize(A)
LU{Float64, Tridiagonal{Float64, Vector{Float64}}}
L factor:
5×5 Matrix{Float64}:
  1.0   0.0        0.0    0.0  0.0
 -0.5   1.0        0.0    0.0  0.0
  0.0  -0.666667   1.0    0.0  0.0
  0.0   0.0       -0.75   1.0  0.0
  0.0   0.0        0.0   -0.8  1.0
U factor:
5×5 Matrix{Float64}:
 2.0  -1.0   0.0       0.0    0.0
 0.0   1.5  -1.0       0.0    0.0
 0.0   0.0   1.33333  -1.0    0.0
 0.0   0.0   0.0       1.25  -1.0
 0.0   0.0   0.0       0.0    1.2

julia> using SpecialMatrices

julia> B = Strang(5)
5×5 Strang{Float64}:
  2  -1   0   0   0
 -1   2  -1   0   0
  0  -1   2  -1   0
  0   0  -1   2  -1
  0   0   0  -1   2

julia> factorize(B)
ERROR: MethodError: no method matching factorize(::Strang{Float64})
Closest candidates are:
  factorize(::StridedMatrix{T}) where T at C:\Users\Braam\.julia\juliaup\julia-1.7.2+0~x64\share\julia\stdlib\v1.7\LinearAlgebra\src\dense.jl:1302
  factorize(::Adjoint) at C:\Users\Braam\.julia\juliaup\julia-1.7.2+0~x64\share\julia\stdlib\v1.7\LinearAlgebra\src\dense.jl:1376
  factorize(::Transpose) at C:\Users\Braam\.julia\juliaup\julia-1.7.2+0~x64\share\julia\stdlib\v1.7\LinearAlgebra\src\dense.jl:1377
  ...
Stacktrace:
 [1] top-level scope
   @ REPL[17]:1

[JeffFessler]

Thanks for the report. Probably Strang should just use the SymTridiagonal type from LinearAlgebra instead of re-inventing the wheel. There is a WIP PR #51 that will address this issue. Comments on that PR are welcome! I will leave this open until we merge that PR.

[JeffFessler]

@braamvandyk PR #78 should address this issue as long as you are using Julia >= v1.8. If you are using an earlier version feel free to make a PR to handle that version. (I don't know why it fails with 1.6 and don't have time to figure it out.) Thanks for reporting this issue - I hadn't realized that Strang has such a simple LDL' factorization!

You can see a verification of it in the docs here: https://julialinearalgebra.github.io/SpecialMatrices.jl/dev/generated/examples/1-overview

So I close the issue but please reopen if I've overlooked something.