ClassicalOrthogonalPolynomials.jl icon indicating copy to clipboard operation
ClassicalOrthogonalPolynomials.jl copied to clipboard

Truncating the identity operator loses sparsity

Open putianyi889 opened this issue 1 year ago • 1 comments

To avoid type piracy, Jacobi(0,1)\Jacobi(0,1) should return a (0,0)-banded matrix.

julia> Jacobi(0,1)\Jacobi(0,1)
ℵ₀×ℵ₀ Diagonal{Float64, FillArrays.Ones{Float64, 1, Tuple{InfiniteArrays.OneToInf{Int64}}}} with indices OneToInf()×OneToInf()

julia> (Jacobi(0,1)\Jacobi(0,1))[1:10,1:10]
10×10 Matrix{Float64}:
 1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  1.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  1.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  1.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  1.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  1.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  1.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  1.0

julia> (Jacobi(0,1)\Jacobi(0,0))[1:10,1:10]
10×10 BandedMatrix{Float64} with bandwidths (0, 1):
 1.0  0.333333   ⋅    ⋅         ⋅         ⋅         ⋅         ⋅         ⋅         ⋅ 
  ⋅   0.666667  0.4   ⋅         ⋅         ⋅         ⋅         ⋅         ⋅         ⋅
  ⋅    ⋅        0.6  0.428571   ⋅         ⋅         ⋅         ⋅         ⋅         ⋅
  ⋅    ⋅         ⋅   0.571429  0.444444   ⋅         ⋅         ⋅         ⋅         ⋅ 
  ⋅    ⋅         ⋅    ⋅        0.555556  0.454545   ⋅         ⋅         ⋅         ⋅
  ⋅    ⋅         ⋅    ⋅         ⋅        0.545455  0.461538   ⋅         ⋅         ⋅
  ⋅    ⋅         ⋅    ⋅         ⋅         ⋅        0.538462  0.466667   ⋅         ⋅
  ⋅    ⋅         ⋅    ⋅         ⋅         ⋅         ⋅        0.533333  0.470588   ⋅
  ⋅    ⋅         ⋅    ⋅         ⋅         ⋅         ⋅         ⋅        0.529412  0.473684
  ⋅    ⋅         ⋅    ⋅         ⋅         ⋅         ⋅         ⋅         ⋅        0.526316

putianyi889 avatar Aug 13 '24 06:08 putianyi889

Just use layout_getindex for now to avoid this

dlfivefifty avatar Oct 21 '24 08:10 dlfivefifty