Infinite loop converting between ChebyshevU and Ultraspherical
Closed this issue · 1 comments
DanielVandH commented
Noticed accidentally that the code
julia> using ClassicalOrthogonalPolynomials
julia> U = ChebyshevU();
julia> U \ diff(U) never finishes. Ctrl+Cing this gives
julia> U \ diff(U)
ERROR: InterruptException:
Stacktrace:
[1] cache_filldata!(::LazyArrays.CachedArray{…}, ::UnitRange{…}, ::Base.OneTo{…})
@ LazyArrays C:\Users\User\.julia\packages\LazyArrays\JNdsG\src\cache.jl:222
[2] resizedata!(::ArrayLayouts.DenseColumnMajor, ::ArrayLayouts.ZerosLayout, ::LazyArrays.CachedArray{…}, ::Int64, ::Int64)
@ LazyArrays C:\Users\User\.julia\packages\LazyArrays\JNdsG\src\cache.jl:263
[3] resizedata!
@ C:\Users\User\.julia\packages\LazyArrays\JNdsG\src\cache.jl:218 [inlined]
[4] setindex!
@ C:\Users\User\.julia\packages\LazyArrays\JNdsG\src\cache.jl:82 [inlined]
[5] setindex!
@ .\subarray.jl:331 [inlined]
[6] _setindex!
@ .\abstractarray.jl:1426 [inlined]
[7] setindex!
@ .\abstractarray.jl:1396 [inlined]
[8] copyto_unaliased!(deststyle::IndexCartesian, dest::SubArray{…}, srcstyle::IndexCartesian, src::SubArray{…})
@ Base .\abstractarray.jl:1102
[9] copyto!
@ .\abstractarray.jl:1068 [inlined]
[10] _copyto!
@ C:\Users\User\.julia\packages\ArrayLayouts\ccyJu\src\ArrayLayouts.jl:255 [inlined]
[11] _copyto!
@ C:\Users\User\.julia\packages\ArrayLayouts\ccyJu\src\ArrayLayouts.jl:259 [inlined]
[12] copyto!
@ C:\Users\User\.julia\packages\ArrayLayouts\ccyJu\src\ArrayLayouts.jl:264 [inlined]
[13] ldiv!(Y::LazyArrays.CachedArray{…}, A::MatrixFactorizations.QR{…}, B::BandedMatrices.BandedMatrix{…})
@ LinearAlgebra C:\Users\User\.julia\juliaup\julia-1.10.3+0.x64.w64.mingw32\share\julia\stdlib\v1.10\LinearAlgebra\src\factorization.jl:179
[14] _ldiv!
@ C:\Users\User\.julia\packages\ArrayLayouts\ccyJu\src\ldiv.jl:84 [inlined]
[15] copyto!(dest::LazyArrays.CachedArray{…}, M::ArrayLayouts.Ldiv{…}; kwds::@Kwargs{})
@ ArrayLayouts C:\Users\User\.julia\packages\ArrayLayouts\ccyJu\src\ldiv.jl:113
[16] copyto!
@ C:\Users\User\.julia\packages\ArrayLayouts\ccyJu\src\ldiv.jl:113 [inlined]
[17] ldiv!
@ C:\Users\User\.julia\packages\ArrayLayouts\ccyJu\src\ldiv.jl:104 [inlined]
[18] _ldiv!
@ C:\Users\User\.julia\packages\ArrayLayouts\ccyJu\src\ldiv.jl:83 [inlined]
[19] copyto!
@ C:\Users\User\.julia\packages\ArrayLayouts\ccyJu\src\ldiv.jl:113 [inlined]
[20] copy
@ C:\Users\User\.julia\packages\ArrayLayouts\ccyJu\src\ldiv.jl:21 [inlined]
[21] materialize
@ C:\Users\User\.julia\packages\ArrayLayouts\ccyJu\src\ldiv.jl:22 [inlined]
[22] ldiv
@ C:\Users\User\.julia\packages\ArrayLayouts\ccyJu\src\ldiv.jl:98 [inlined]
[23] copy(M::ArrayLayouts.Mul{…})
@ LazyArrays C:\Users\User\.julia\packages\LazyArrays\JNdsG\src\linalg\inv.jl:90
[24] materialize
@ C:\Users\User\.julia\packages\ArrayLayouts\ccyJu\src\mul.jl:127 [inlined]
[25] mul
@ C:\Users\User\.julia\packages\ArrayLayouts\ccyJu\src\mul.jl:128 [inlined]
[26] *(A::LazyArrays.ApplyArray{…}, B::BandedMatrices.BandedMatrix{…})
@ ArrayLayouts C:\Users\User\.julia\packages\ArrayLayouts\ccyJu\src\mul.jl:213
[27] _copy_ldiv_mul
@ C:\Users\User\.julia\packages\LazyArrays\JNdsG\src\linalg\inv.jl:112 [inlined]
[28] copy
@ C:\Users\User\.julia\packages\LazyArrays\JNdsG\src\linalg\inv.jl:116 [inlined]
[29] basis_ldiv_size
@ C:\Users\User\.julia\packages\ContinuumArrays\0grIp\src\bases\bases.jl:327 [inlined]
[30] copy
@ C:\Users\User\.julia\packages\ContinuumArrays\0grIp\src\bases\bases.jl:323 [inlined]
[31] materialize
@ C:\Users\User\.julia\packages\ArrayLayouts\ccyJu\src\ldiv.jl:22 [inlined]
[32] ldiv
@ C:\Users\User\.julia\packages\ArrayLayouts\ccyJu\src\ldiv.jl:98 [inlined]
[33] \(A::ChebyshevU{Float64}, B::QuasiArrays.ApplyQuasiMatrix{Float64, typeof(*), Tuple{…}})
@ QuasiArrays C:\Users\User\.julia\packages\QuasiArrays\ZCn0F\src\matmul.jl:34
[34] top-level scope
@ REPL[9]:1
Some type information was truncated. Use `show(err)` to see complete types.I'm guessing that this should give a "Not implemented error" since the equivalent does
julia> Ultraspherical(1) \ diff(U)
ERROR: Not implemented
Stacktrace:
[1] error(s::String)
@ Base .\error.jl:35
[2] \(C2::Ultraspherical{Float64, Int64}, C1::Ultraspherical{Float64, Int64})
@ ClassicalOrthogonalPolynomials C:\Users\User\.julia\packages\ClassicalOrthogonalPolynomials\Hdhud\src\classical\ultraspherical.jl:211
[3] copy
@ C:\Users\User\.julia\packages\ContinuumArrays\0grIp\src\bases\bases.jl:0 [inlined]
[4] materialize
@ C:\Users\User\.julia\packages\ArrayLayouts\ccyJu\src\ldiv.jl:22 [inlined]
[5] ldiv
@ C:\Users\User\.julia\packages\ArrayLayouts\ccyJu\src\ldiv.jl:98 [inlined]
[6] \(A::Ultraspherical{Float64, Int64}, B::QuasiArrays.ApplyQuasiMatrix{Float64, typeof(*), Tuple{…}})
@ QuasiArrays C:\Users\User\.julia\packages\QuasiArrays\ZCn0F\src\matmul.jl:34
[7] top-level scope
@ REPL[11]:1
Some type information was truncated. Use `show(err)` to see complete types.dlfivefifty commented
The issue is the following: we are multiplying R and D that look like:
julia> R
inv(ℵ₀×ℵ₀ BandedMatrix{Float64} with bandwidths (0, 2) with data vcat(((-).((Float64) ./ (ℵ₀-element InfiniteArrays.InfStepRange{Float64, Float64} with indices OneToInf()) with indices OneToInf()) with indices OneToInf())' with indices Base.OneTo(1)×OneToInf(), 1×ℵ₀ FillArrays.Zeros{Float64, 2, Tuple{Base.OneTo{Int64}, InfiniteArrays.OneToInf{Int64}}} with indices Base.OneTo(1)×OneToInf(), ((Float64) ./ (ℵ₀-element InfiniteArrays.InfStepRange{Float64, Float64} with indices OneToInf()) with indices OneToInf())' with indices Base.OneTo(1)×OneToInf()) with indices Base.OneTo(3)×OneToInf() with indices OneToInf()×OneToInf()) with indices OneToInf()×OneToInf():
1.0 0.0 1.0 0.0 1.0 0.0 1.0 0.0 …
⋅ 2.0 0.0 2.0 0.0 2.0 0.0 2.0
⋅ ⋅ 3.0 0.0 3.0 0.0 3.0 0.0
⋅ ⋅ ⋅ 4.0 0.0 4.0 0.0 4.0
⋅ ⋅ ⋅ ⋅ 5.0 0.0 5.0 0.0
⋅ ⋅ ⋅ ⋅ ⋅ 6.0 0.0 6.0 …
⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 7.0 0.0
⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 8.0
⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅
⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅
⋮ ⋮ ⋱
julia> D
ℵ₀×ℵ₀ BandedMatrix{Float64} with bandwidths (-1, 1) with data 1×ℵ₀ FillArrays.Fill{Float64, 2, Tuple{Base.OneTo{Int64}, InfiniteArrays.OneToInf{Int64}}} with indices Base.OneTo(1)×OneToInf() with indices OneToInf()×OneToInf():
⋅ 2.0 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ …
⋅ ⋅ 2.0 ⋅ ⋅ ⋅ ⋅ ⋅
⋅ ⋅ ⋅ 2.0 ⋅ ⋅ ⋅ ⋅
⋅ ⋅ ⋅ ⋅ 2.0 ⋅ ⋅ ⋅
⋅ ⋅ ⋅ ⋅ ⋅ 2.0 ⋅ ⋅
⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 2.0 ⋅ …
⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 2.0
⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅
⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅
⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅
⋮ ⋮ ⋱ Their layouts are:
julia> MemoryLayout(R)
LazyArrays.InvLayout{BandedMatrices.BandedColumns{LazyArrays.LazyLayout}}()
julia> MemoryLayout(D)
BandedMatrices.BandedColumns{FillLayout}()We just need to add InvLayout{<:LazyBandedLayouts} to the overloads that leave multiplication as lazy. But Because this is all about to be moved from LazyBandedMatrices.jl to extensions in JuliaArrays/LazyArrays.jl#303 we can't fix it right now.