search

JuliaApproximation / SemiclassicalOrthogonalPolynomials.jl

A Julia repository for semiclassical orthogonal polynomials
MIT License
7 stars 3 forks source link

Support computing connection matrices when decrementing parameters #109

Closed DanielVandH closed 2 months ago

DanielVandH commented 2 months ago

This adds some support for decrementing parameters. The implementation has something wrong with it apparently.. the tests pass, but occasionally I reach a segfault. Any ideas? The block

    @testset "Simultaneous" begin 
        A = SemiclassicalJacobi(t, a, b, c)
        BB = SemiclassicalJacobi.(t, a .+ (-2:2), b .+ (-2:2), c .+ (-2:2))
        for B in BB 
            @test (A \ B)[1:100, 1:100] ≈ ApplyArray(inv, B \ A)[1:100, 1:100]
        end
    end

usually always segfaults, but it's not consistent where it does. I've seen one segfault in the Decrementing a testset but I can't seem to reproduce any of this reliably. I tried changing the calls to use the _direct method but the same thing occurs.

dlfivefifty commented 2 months ago

Try running with julia --check-bounds=yes

DanielVandH commented 2 months ago

With that option I now get no segfault but I also get no BoundsErrors... it seems to be related to trying to index into the inverse:

julia> A = SemiclassicalJacobi(2.0, 3.0, 4.0, 5.0);

julia> B = SemiclassicalJacobi(2.0, 4.0, 5.0, 6.0);

julia> R = A \ B; # fine

julia> R[1, 1] # not fine
ERROR: InterruptException:
Stacktrace:
  [1] getindex(A::LazyBandedMatrices.SymTridiagonal{Float64, SubArray{…}, SubArray{…}}, i::Int64, j::Int64)
    @ LazyBandedMatrices C:\Users\User\.julia\packages\LazyBandedMatrices\gqWs9\src\tridiag.jl:301
  [2] getindex
    @ .\subarray.jl:290 [inlined]
  [3] (BandedMatrix{…})(A::SubArray{…}, bnds::Tuple{…})
    @ BandedMatrices C:\Users\User\.julia\packages\BandedMatrices\d4g9c\src\banded\BandedMatrix.jl:214
  [4] BandedMatrix
    @ C:\Users\User\.julia\packages\BandedMatrices\d4g9c\src\banded\BandedMatrix.jl:207 [inlined]
  [5] BandedMatrix
    @ C:\Users\User\.julia\packages\BandedMatrices\d4g9c\src\banded\BandedMatrix.jl:219 [inlined]
  [6] _BandedMatrix
    @ C:\Users\User\.julia\packages\BandedMatrices\d4g9c\src\banded\BandedMatrix.jl:245 [inlined]
  [7] BandedMatrix
    @ C:\Users\User\.julia\packages\BandedMatrices\d4g9c\src\banded\BandedMatrix.jl:246 [inlined]
  [8] sub_materialize
    @ C:\Users\User\.julia\packages\BandedMatrices\d4g9c\src\generic\indexing.jl:28 [inlined]
  [9] sub_materialize
    @ C:\Users\User\.julia\packages\ArrayLayouts\3byqH\src\ArrayLayouts.jl:131 [inlined]
 [10] sub_materialize
    @ C:\Users\User\.julia\packages\ArrayLayouts\3byqH\src\ArrayLayouts.jl:132 [inlined]
 [11] copy
    @ C:\Users\User\.julia\packages\ArrayLayouts\3byqH\src\ArrayLayouts.jl:135 [inlined]
 [12] getindex(C::LazyBandedMatrices.SymTridiagonal{…}, kr::UnitRange{…}, jr::UnitRange{…})
    @ ClassicalOrthogonalPolynomials C:\Users\User\.julia\packages\ClassicalOrthogonalPolynomials\JEEfo\src\choleskyQR.jl:90
 [13] _fillqrbanddata!(J::ClassicalOrthogonalPolynomials.QRJacobiData{:Q, Float64}, inds::UnitRange{Int64})
    @ ClassicalOrthogonalPolynomials C:\Users\User\.julia\packages\ClassicalOrthogonalPolynomials\JEEfo\src\choleskyQR.jl:238
 [14] resizedata!(K::ClassicalOrthogonalPolynomials.QRJacobiData{:Q, Float64}, n::Int64, m::Int64)
    @ ClassicalOrthogonalPolynomials C:\Users\User\.julia\packages\ClassicalOrthogonalPolynomials\JEEfo\src\choleskyQR.jl:229
 [15] getindex(C::LazyBandedMatrices.SymTridiagonal{…}, kr::UnitRange{…}, jr::UnitRange{…})
    @ ClassicalOrthogonalPolynomials C:\Users\User\.julia\packages\ClassicalOrthogonalPolynomials\JEEfo\src\choleskyQR.jl:217
 [16] _fillqrbanddata!(J::ClassicalOrthogonalPolynomials.QRJacobiData{:Q, Float64}, inds::UnitRange{Int64})
    @ ClassicalOrthogonalPolynomials C:\Users\User\.julia\packages\ClassicalOrthogonalPolynomials\JEEfo\src\choleskyQR.jl:238
 [17] resizedata!(K::ClassicalOrthogonalPolynomials.QRJacobiData{:Q, Float64}, n::Int64, m::Int64)
    @ ClassicalOrthogonalPolynomials C:\Users\User\.julia\packages\ClassicalOrthogonalPolynomials\JEEfo\src\choleskyQR.jl:229
 [18] getindex(K::ClassicalOrthogonalPolynomials.QRJacobiData{:Q, Float64}, n::Int64, m::Int64)
    @ ClassicalOrthogonalPolynomials C:\Users\User\.julia\packages\ClassicalOrthogonalPolynomials\JEEfo\src\choleskyQR.jl:210
 [19] getindex
    @ .\subarray.jl:290 [inlined]
 [20] _getindex
    @ .\abstractarray.jl:1341 [inlined]
 [21] getindex
    @ .\abstractarray.jl:1291 [inlined]
 [22] _broadcast_getindex
    @ .\broadcast.jl:662 [inlined]
 [23] _getindex
    @ .\broadcast.jl:706 [inlined]
 [24] _getindex
    @ .\broadcast.jl:705 [inlined]
 [25] _broadcast_getindex
    @ .\broadcast.jl:681 [inlined]
 [26] getindex
    @ .\broadcast.jl:636 [inlined]
 [27] getindex
    @ C:\Users\User\.julia\packages\LazyArrays\MLFsy\src\lazybroadcasting.jl:93 [inlined]
 [28] getindex(A::LazyBandedMatrices.SymTridiagonal{Float64, BroadcastVector{…}, BroadcastVector{…}}, i::Int64, j::Int64)
    @ LazyBandedMatrices C:\Users\User\.julia\packages\LazyBandedMatrices\gqWs9\src\tridiag.jl:306
 [29] getindex
    @ .\subarray.jl:290 [inlined]
 [30] (BandedMatrix{…})(A::SubArray{…}, bnds::Tuple{…})
    @ BandedMatrices C:\Users\User\.julia\packages\BandedMatrices\d4g9c\src\banded\BandedMatrix.jl:214
 [31] BandedMatrix
    @ C:\Users\User\.julia\packages\BandedMatrices\d4g9c\src\banded\BandedMatrix.jl:207 [inlined]
 [32] BandedMatrix
    @ C:\Users\User\.julia\packages\BandedMatrices\d4g9c\src\banded\BandedMatrix.jl:219 [inlined]
 [33] _BandedMatrix
    @ C:\Users\User\.julia\packages\BandedMatrices\d4g9c\src\banded\BandedMatrix.jl:245 [inlined]
 [34] BandedMatrix
    @ C:\Users\User\.julia\packages\BandedMatrices\d4g9c\src\banded\BandedMatrix.jl:246 [inlined]
 [35] sub_materialize
    @ C:\Users\User\.julia\packages\BandedMatrices\d4g9c\src\generic\indexing.jl:28 [inlined]
 [36] sub_materialize
    @ C:\Users\User\.julia\packages\ArrayLayouts\3byqH\src\ArrayLayouts.jl:131 [inlined]
 [37] sub_materialize
    @ C:\Users\User\.julia\packages\ArrayLayouts\3byqH\src\ArrayLayouts.jl:132 [inlined]
 [38] layout_getindex
    @ C:\Users\User\.julia\packages\ArrayLayouts\3byqH\src\ArrayLayouts.jl:138 [inlined]
 [39] getindex
    @ C:\Users\User\.julia\packages\ArrayLayouts\3byqH\src\ArrayLayouts.jl:153 [inlined]
 [40] resizedata!(::BandedMatrices.BandedColumns{…}, ::SymTridiagonalLayout{…}, B::LazyArrays.CachedArray{…}, n::Int64, m::Int64)
    @ LazyArraysBandedMatricesExt C:\Users\User\.julia\packages\LazyArrays\MLFsy\ext\LazyArraysBandedMatricesExt.jl:492
 [41] resizedata!
    @ C:\Users\User\.julia\packages\LazyArrays\MLFsy\src\cache.jl:218 [inlined]
 [42] partialcholesky!(F::InfiniteLinearAlgebra.AdaptiveCholeskyFactors{…}, n::Int64)
    @ InfiniteLinearAlgebra C:\Users\User\.julia\packages\InfiniteLinearAlgebra\cI65l\src\infcholesky.jl:28
 [43] ldiv!(R::UpperTriangular{…}, B::LazyArrays.CachedArray{…})
    @ InfiniteLinearAlgebra C:\Users\User\.julia\packages\InfiniteLinearAlgebra\cI65l\src\infcholesky.jl:115
 [44] copyto!
    @ C:\Users\User\.julia\packages\ArrayLayouts\3byqH\src\triangular.jl:210 [inlined]
 [45] copy
    @ C:\Users\User\.julia\packages\ArrayLayouts\3byqH\src\ldiv.jl:21 [inlined]
 [46] materialize
    @ C:\Users\User\.julia\packages\ArrayLayouts\3byqH\src\ldiv.jl:22 [inlined]
 [47] ldiv
    @ C:\Users\User\.julia\packages\ArrayLayouts\3byqH\src\ldiv.jl:98 [inlined]
 [48] \
    @ C:\Users\User\.julia\packages\ArrayLayouts\3byqH\src\ldiv.jl:199 [inlined]
 [49] _copy_ldiv_ldiv
    @ C:\Users\User\.julia\packages\LazyArrays\MLFsy\src\linalg\inv.jl:120 [inlined]
 [50] _copy_ldiv_ldiv
    @ C:\Users\User\.julia\packages\LazyArrays\MLFsy\src\linalg\inv.jl:121 [inlined]
 [51] copy
    @ C:\Users\User\.julia\packages\LazyArrays\MLFsy\src\linalg\inv.jl:122 [inlined]
 [52] materialize(M::Ldiv{…}; kwds::@Kwargs{})
    @ ArrayLayouts C:\Users\User\.julia\packages\ArrayLayouts\3byqH\src\ldiv.jl:22
 [53] ldiv(A::ApplyArray{…}, B::LazyArrays.CachedArray{…}; kwds::@Kwargs{})
    @ ArrayLayouts C:\Users\User\.julia\packages\ArrayLayouts\3byqH\src\ldiv.jl:98
 [54] ldiv
    @ C:\Users\User\.julia\packages\ArrayLayouts\3byqH\src\ldiv.jl:98 [inlined]
 [55] \
    @ C:\Users\User\.julia\packages\ArrayLayouts\3byqH\src\ldiv.jl:199 [inlined]
 [56] sub_materialize
    @ C:\Users\User\.julia\packages\LazyArrays\MLFsy\src\linalg\inv.jl:173 [inlined]
 [57] sub_materialize
    @ C:\Users\User\.julia\packages\InfiniteArrays\7v9AI\src\infarrays.jl:399 [inlined]
 [58] sub_materialize
    @ C:\Users\User\.julia\packages\ArrayLayouts\3byqH\src\ArrayLayouts.jl:131 [inlined]
 [59] sub_materialize
    @ C:\Users\User\.julia\packages\ArrayLayouts\3byqH\src\ArrayLayouts.jl:132 [inlined]
 [60] layout_getindex
    @ C:\Users\User\.julia\packages\ArrayLayouts\3byqH\src\ArrayLayouts.jl:138 [inlined]
 [61] getindex
    @ C:\Users\User\.julia\packages\ArrayLayouts\3byqH\src\ArrayLayouts.jl:156 [inlined]
 [62] getindex
    @ C:\Users\User\.julia\packages\LazyArrays\MLFsy\src\linalg\inv.jl:177 [inlined]
 [63] getindex
    @ C:\Users\User\.julia\juliaup\julia-1.10.4+0.x64.w64.mingw32\share\julia\stdlib\v1.10\LinearAlgebra\src\adjtrans.jl:329 [inlined]
 [64] getindex
    @ .\subarray.jl:290 [inlined]
 [65] _default_blasmul!(::IndexCartesian, α::Bool, A::SubArray{…}, B::Vector{…}, β::Bool, C::LazyArrays.CachedArray{…})
    @ ArrayLayouts C:\Users\User\.julia\packages\ArrayLayouts\3byqH\src\muladd.jl:252
 [66] default_blasmul!
    @ C:\Users\User\.julia\packages\ArrayLayouts\3byqH\src\muladd.jl:259 [inlined]
 [67] materialize!(M::MulAdd{…})
    @ ArrayLayouts C:\Users\User\.julia\packages\ArrayLayouts\3byqH\src\muladd.jl:289
 [68] muladd!
    @ C:\Users\User\.julia\packages\ArrayLayouts\3byqH\src\muladd.jl:75 [inlined]
 [69] mul!
    @ C:\Users\User\.julia\packages\ArrayLayouts\3byqH\src\mul.jl:144 [inlined]
 [70] mul!
    @ C:\Users\User\.julia\packages\ArrayLayouts\3byqH\src\mul.jl:386 [inlined]
 [71] mul!
    @ C:\Users\User\.julia\juliaup\julia-1.10.4+0.x64.w64.mingw32\share\julia\stdlib\v1.10\LinearAlgebra\src\matmul.jl:237 [inlined]
 [72] *(A::SubArray{Float64, 2, Transpose{…}, Tuple{…}, false}, x::Vector{Float64})
    @ LinearAlgebra C:\Users\User\.julia\juliaup\julia-1.10.4+0.x64.w64.mingw32\share\julia\stdlib\v1.10\LinearAlgebra\src\matmul.jl:57
 [73] sub_materialize
    @ C:\Users\User\.julia\packages\LazyArrays\MLFsy\src\linalg\mul.jl:261 [inlined]
 [74] sub_materialize
    @ C:\Users\User\.julia\packages\ArrayLayouts\3byqH\src\ArrayLayouts.jl:132 [inlined]
 [75] map(f::typeof(ArrayLayouts.sub_materialize), t::Tuple{SubArray{…}, SubArray{…}})
    @ Base .\tuple.jl:292
 [76] sub_materialize
    @ C:\Users\User\.julia\packages\LazyArrays\MLFsy\src\linalg\mul.jl:261 [inlined]
 [77] sub_materialize
    @ C:\Users\User\.julia\packages\ArrayLayouts\3byqH\src\ArrayLayouts.jl:132 [inlined]
 [78] layout_getindex
    @ C:\Users\User\.julia\packages\ArrayLayouts\3byqH\src\ArrayLayouts.jl:138 [inlined]
 [79] getindex
    @ C:\Users\User\.julia\packages\ArrayLayouts\3byqH\src\ArrayLayouts.jl:159 [inlined]
 [80] map
    @ .\tuple.jl:342 [inlined]
 [81] _mul_getindex(args::Tuple{ApplyArray{…}, ApplyArray{…}}, k::Int64, j::Int64)
    @ LazyArrays C:\Users\User\.julia\packages\LazyArrays\MLFsy\src\linalg\mul.jl:218
 [82] getindex(::ApplyArray{Float64, 2, typeof(*), Tuple{ApplyArray{…}, ApplyArray{…}}}, ::Int64, ::Int64)
    @ LazyArrays C:\Users\User\.julia\packages\LazyArrays\MLFsy\src\linalg\mul.jl:184
Some type information was truncated. Use `show(err)` to see complete types.

julia> MemoryLayout(R)
LazyArrays.ApplyLayout{typeof(*)}()

Maybe related to #108's layout issues

DanielVandH commented 2 months ago

Not 100% sure but it seems the segfault goes away when I remove threads here

https://github.com/JuliaApproximation/ClassicalOrthogonalPolynomials.jl/blob/52564bfd4fc33667ebca6d8f2cb35c21f5e38426/src/choleskyQR.jl#L264-L276

TSGut commented 2 months ago

Hm, the resizing prior to the loop should make that safe if the sizes are correctly chosen, is this with --check-bounds=yes?

DanielVandH commented 2 months ago

Without --check-bounds=yes and with threads, I saw

  copyto! with padded/zeros only supported with equal axes
  Stacktrace:
    [1] error(s::String)
      @ Base .\error.jl:35
    [2] _copyto!(::LazyArrays.PaddedColumns{ArrayLayouts.DenseColumnMajor}, ::ArrayLayouts.ZerosLayout, dest::SubArray{Float64, 1, LazyArrays.CachedArray{Float64, 1, Vector{Float64}, Zeros{Float64, 1, Tuple{InfiniteArrays.OneToInf{Int64}}}}, Tuple{UnitRange{Int64}}, false}, src::Zeros{Float64, 1, Tuple{Base.OneTo{Int64}}})
      @ LazyArrays C:\Users\User\.julia\packages\LazyArrays\MLFsy\src\padded.jl:154
    [3] _copyto!
      @ C:\Users\User\.julia\packages\ArrayLayouts\3byqH\src\ArrayLayouts.jl:263 [inlined]
    [4] copyto!(dest::SubArray{Float64, 1, LazyArrays.CachedArray{Float64, 1, Vector{Float64}, Zeros{Float64, 1, Tuple{InfiniteArrays.OneToInf{Int64}}}}, Tuple{UnitRange{Int64}}, false}, src::Zeros{Float64, 1, Tuple{Base.OneTo{Int64}}})
      @ ArrayLayouts C:\Users\User\.julia\packages\ArrayLayouts\3byqH\src\ArrayLayouts.jl:271
    [5] __cat(::LazyArrays.CachedArray{Float64, 1, Vector{Float64}, Zeros{Float64, 1, Tuple{InfiniteArrays.OneToInf{Int64}}}}, ::Tuple{Infinities.InfiniteCardinal{0}}, ::Tuple{Bool}, ::Zeros{Float64, 1, Tuple{Base.OneTo{Int64}}}, ::Vararg{Any})
      @ InfiniteArrays C:\Users\User\.julia\packages\InfiniteArrays\7v9AI\src\InfiniteArrays.jl:126
    [6] _cat_t
      @ .\abstractarray.jl:1805 [inlined]
    [7] typed_vcat
      @ .\abstractarray.jl:1946 [inlined]
    [8] vcat
      @ C:\Users\User\.julia\juliaup\julia-1.10.4+0.x64.w64.mingw32\share\julia\stdlib\v1.10\SparseArrays\src\sparsevector.jl:1239 [inlined]
    [9] sub_materialize
      @ C:\Users\User\.julia\packages\LazyArrays\MLFsy\src\linalg\inv.jl:173 [inlined]
   [10] sub_materialize
      @ C:\Users\User\.julia\packages\InfiniteArrays\7v9AI\src\infarrays.jl:399 [inlined]
   [11] sub_materialize
      @ C:\Users\User\.julia\packages\ArrayLayouts\3byqH\src\ArrayLayouts.jl:131 [inlined]
   [12] sub_materialize
      @ C:\Users\User\.julia\packages\ArrayLayouts\3byqH\src\ArrayLayouts.jl:132 [inlined]
   [13] layout_getindex
      @ C:\Users\User\.julia\packages\ArrayLayouts\3byqH\src\ArrayLayouts.jl:138 [inlined]
   [14] getindex
      @ C:\Users\User\.julia\packages\ArrayLayouts\3byqH\src\ArrayLayouts.jl:156 [inlined]
   [15] getindex
      @ C:\Users\User\.julia\packages\LazyArrays\MLFsy\src\linalg\inv.jl:177 [inlined]
   [16] getindex
      @ .\subarray.jl:290 [inlined]
   [17] _getindex
      @ .\abstractarray.jl:1341 [inlined]
   [18] getindex
      @ .\abstractarray.jl:1291 [inlined]
   [19] iterate(A::SubArray{Float64, 2, ApplyArray{Float64, 2, typeof(inv), Tuple{ApplyArray{Float64, 2, typeof(*), Tuple{Diagonal{Float64, BroadcastVector{Float64, typeof(*), Tuple{BroadcastVector{Float64, typeof(sign), Tuple{SubArray{Float64, 1, UpperTriangular{Float64, InfiniteLinearAlgebra.AdaptiveQRFactors{Float64, BandedMatrix{Float64, Matrix{Float64}, Base.OneTo{Int64}}, LazyBandedMatrices.SymTridiagonal{Float64, SubArray{Float64, 1, ClassicalOrthogonalPolynomials.QRJacobiData{:Q, Float64}, Tuple{Base.Slice{InfiniteArrays.OneToInf{Int64}}, Int64}, false}, SubArray{Float64, 1, ClassicalOrthogonalPolynomials.QRJacobiData{:Q, Float64}, Tuple{Base.Slice{InfiniteArrays.OneToInf{Int64}}, Int64}, false}}}}, Tuple{InfiniteLinearAlgebra.InfBandCartesianIndices}, false}}}, Fill{Float64, 1, Tuple{InfiniteArrays.OneToInf{Int64}}}}}}, UpperTriangular{Float64, InfiniteLinearAlgebra.AdaptiveQRFactors{Float64, BandedMatrix{Float64, Matrix{Float64}, Base.OneTo{Int64}}, LazyBandedMatrices.SymTridiagonal{Float64, SubArray{Float64, 1, ClassicalOrthogonalPolynomials.QRJacobiData{:Q, Float64}, Tuple{Base.Slice{InfiniteArrays.OneToInf{Int64}}, Int64}, false}, SubArray{Float64, 1, ClassicalOrthogonalPolynomials.QRJacobiData{:Q, Float64}, Tuple{Base.Slice{InfiniteArrays.OneToInf{Int64}}, Int64}, false}}}}}}}}, Tuple{UnitRange{Int64}, UnitRange{Int64}}, false}, state::Tuple{CartesianIndices{2, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}}, CartesianIndex{2}})
      @ Base .\abstractarray.jl:1217
   [20] copyto_unaliased!(deststyle::IndexLinear, dest::Matrix{Float64}, srcstyle::IndexCartesian, src::SubArray{Float64, 2, ApplyArray{Float64, 2, typeof(inv), Tuple{ApplyArray{Float64, 2, typeof(*), Tuple{Diagonal{Float64, BroadcastVector{Float64, typeof(*), Tuple{BroadcastVector{Float64, typeof(sign), Tuple{SubArray{Float64, 1, UpperTriangular{Float64, InfiniteLinearAlgebra.AdaptiveQRFactors{Float64, BandedMatrix{Float64, Matrix{Float64}, Base.OneTo{Int64}}, LazyBandedMatrices.SymTridiagonal{Float64, SubArray{Float64, 1, ClassicalOrthogonalPolynomials.QRJacobiData{:Q, Float64}, Tuple{Base.Slice{InfiniteArrays.OneToInf{Int64}}, Int64}, false}, SubArray{Float64, 1, ClassicalOrthogonalPolynomials.QRJacobiData{:Q, Float64}, Tuple{Base.Slice{InfiniteArrays.OneToInf{Int64}}, Int64}, false}}}}, Tuple{InfiniteLinearAlgebra.InfBandCartesianIndices}, false}}}, Fill{Float64, 1, Tuple{InfiniteArrays.OneToInf{Int64}}}}}}, UpperTriangular{Float64, InfiniteLinearAlgebra.AdaptiveQRFactors{Float64, BandedMatrix{Float64, Matrix{Float64}, Base.OneTo{Int64}}, LazyBandedMatrices.SymTridiagonal{Float64, SubArray{Float64, 1, ClassicalOrthogonalPolynomials.QRJacobiData{:Q, Float64}, Tuple{Base.Slice{InfiniteArrays.OneToInf{Int64}}, Int64}, false}, SubArray{Float64, 1, ClassicalOrthogonalPolynomials.QRJacobiData{:Q, Float64}, Tuple{Base.Slice{InfiniteArrays.OneToInf{Int64}}, Int64}, false}}}}}}}}, Tuple{UnitRange{Int64}, UnitRange{Int64}}, false})
      @ Base .\abstractarray.jl:1095
   [21] copyto!
      @ .\abstractarray.jl:1068 [inlined]
   [22] _copyto!
      @ C:\Users\User\.julia\packages\ArrayLayouts\3byqH\src\ArrayLayouts.jl:259 [inlined]
   [23] _copyto!
      @ C:\Users\User\.julia\packages\ArrayLayouts\3byqH\src\ArrayLayouts.jl:263 [inlined]
   [24] copyto!
      @ C:\Users\User\.julia\packages\ArrayLayouts\3byqH\src\ArrayLayouts.jl:272 [inlined]
   [25] copyto_axcheck!
      @ .\abstractarray.jl:1177 [inlined]
   [26] Matrix{Float64}(x::SubArray{Float64, 2, ApplyArray{Float64, 2, typeof(inv), Tuple{ApplyArray{Float64, 2, typeof(*), Tuple{Diagonal{Float64, BroadcastVector{Float64, typeof(*), Tuple{BroadcastVector{Float64, typeof(sign), Tuple{SubArray{Float64, 1, UpperTriangular{Float64, InfiniteLinearAlgebra.AdaptiveQRFactors{Float64, BandedMatrix{Float64, Matrix{Float64}, Base.OneTo{Int64}}, LazyBandedMatrices.SymTridiagonal{Float64, SubArray{Float64, 1, ClassicalOrthogonalPolynomials.QRJacobiData{:Q, Float64}, Tuple{Base.Slice{InfiniteArrays.OneToInf{Int64}}, Int64}, false}, SubArray{Float64, 1, ClassicalOrthogonalPolynomials.QRJacobiData{:Q, Float64}, Tuple{Base.Slice{InfiniteArrays.OneToInf{Int64}}, Int64}, false}}}}, Tuple{InfiniteLinearAlgebra.InfBandCartesianIndices}, false}}}, Fill{Float64, 1, Tuple{InfiniteArrays.OneToInf{Int64}}}}}}, UpperTriangular{Float64, InfiniteLinearAlgebra.AdaptiveQRFactors{Float64, BandedMatrix{Float64, Matrix{Float64}, Base.OneTo{Int64}}, LazyBandedMatrices.SymTridiagonal{Float64, SubArray{Float64, 1, ClassicalOrthogonalPolynomials.QRJacobiData{:Q, Float64}, Tuple{Base.Slice{InfiniteArrays.OneToInf{Int64}}, Int64}, false}, SubArray{Float64, 1, ClassicalOrthogonalPolynomials.QRJacobiData{:Q, Float64}, Tuple{Base.Slice{InfiniteArrays.OneToInf{Int64}}, Int64}, false}}}}}}}}, Tuple{UnitRange{Int64}, UnitRange{Int64}}, false})
      @ Base .\array.jl:673
   [27] Array
      @ .\boot.jl:500 [inlined]
   [28] sub_materialize_axes
      @ C:\Users\User\.julia\packages\ArrayLayouts\3byqH\src\ArrayLayouts.jl:129 [inlined]
   [29] sub_materialize
      @ C:\Users\User\.julia\packages\ArrayLayouts\3byqH\src\ArrayLayouts.jl:130 [inlined]
   [30] sub_materialize
      @ C:\Users\User\.julia\packages\ArrayLayouts\3byqH\src\ArrayLayouts.jl:131 [inlined]
   [31] sub_materialize
      @ C:\Users\User\.julia\packages\ArrayLayouts\3byqH\src\ArrayLayouts.jl:132 [inlined]
   [32] layout_getindex
      @ C:\Users\User\.julia\packages\ArrayLayouts\3byqH\src\ArrayLayouts.jl:138 [inlined]
   [33] getindex(A::ApplyArray{Float64, 2, typeof(inv), Tuple{ApplyArray{Float64, 2, typeof(*), Tuple{Diagonal{Float64, BroadcastVector{Float64, typeof(*), Tuple{BroadcastVector{Float64, typeof(sign), Tuple{SubArray{Float64, 1, UpperTriangular{Float64, InfiniteLinearAlgebra.AdaptiveQRFactors{Float64, BandedMatrix{Float64, Matrix{Float64}, Base.OneTo{Int64}}, LazyBandedMatrices.SymTridiagonal{Float64, SubArray{Float64, 1, ClassicalOrthogonalPolynomials.QRJacobiData{:Q, Float64}, Tuple{Base.Slice{InfiniteArrays.OneToInf{Int64}}, Int64}, false}, SubArray{Float64, 1, ClassicalOrthogonalPolynomials.QRJacobiData{:Q, Float64}, Tuple{Base.Slice{InfiniteArrays.OneToInf{Int64}}, Int64}, false}}}}, Tuple{InfiniteLinearAlgebra.InfBandCartesianIndices}, false}}}, Fill{Float64, 1, Tuple{InfiniteArrays.OneToInf{Int64}}}}}}, UpperTriangular{Float64, InfiniteLinearAlgebra.AdaptiveQRFactors{Float64, BandedMatrix{Float64, Matrix{Float64}, Base.OneTo{Int64}}, LazyBandedMatrices.SymTridiagonal{Float64, SubArray{Float64, 1, ClassicalOrthogonalPolynomials.QRJacobiData{:Q, Float64}, Tuple{Base.Slice{InfiniteArrays.OneToInf{Int64}}, Int64}, false}, SubArray{Float64, 1, ClassicalOrthogonalPolynomials.QRJacobiData{:Q, Float64}, Tuple{Base.Slice{InfiniteArrays.OneToInf{Int64}}, Int64}, false}}}}}}}}, kr::UnitRange{Int64}, jr::UnitRange{Int64})
      @ ArrayLayouts C:\Users\User\.julia\packages\ArrayLayouts\3byqH\src\ArrayLayouts.jl:153
   [34] macro expansion
      @ C:\Users\User\.julia\juliaup\julia-1.10.4+0.x64.w64.mingw32\share\julia\stdlib\v1.10\Test\src\Test.jl:669 [inlined]
   [35] macro expansion
      @ c:\Users\User\.julia\dev\SemiclassicalOrthogonalPolynomials.jl\test\runtests.jl:629 [inlined]
   [36] macro expansion
      @ C:\Users\User\.julia\juliaup\julia-1.10.4+0.x64.w64.mingw32\share\julia\stdlib\v1.10\Test\src\Test.jl:1577 [inlined]
   [37] top-level scope
      @ c:\Users\User\.julia\dev\SemiclassicalOrthogonalPolynomials.jl\test\runtests.jl:626

when running

   t, a, b, c = 2.0, 3.0, 4.0, 5.0
        for _ in 1:1000 # ???
            A = SemiclassicalJacobi(t, a, b, c)
            B = SemiclassicalJacobi(t, a + 2, b, c)
            @test (A \ B)[1:100, 1:100] ≈ ApplyArray(inv, B \ A)[1:100, 1:100]
            A = SemiclassicalJacobi(t, 2.15, 2.05, 3.0111) # check fractions (but integer difference)
            B = SemiclassicalJacobi(t, 3.15, 2.05, 3.0111)
            @test (A \ B)[1:100, 1:100] ≈ ApplyArray(inv, B \ A)[1:100, 1:100]
            A = SemiclassicalJacobi(t, 0.0, 5.0, 2.0) # simultaneous changes
            B = SemiclassicalJacobi(t, 1.0, 4.0, 2.0)
            @test (A \ B)[1:100, 1:100] ≈ ApplyArray(inv, B \ A)[1:100, 1:100]
        end

With --check-bounds=yes and keeping the threads.. nothing happens and the loop finishes. Bizarre

dlfivefifty commented 2 months ago

Why does this for loop use @threads ? I’d be shocked if it’s any faster. And probably it’s not thread safe

TSGut commented 2 months ago

It's faster enough that I kept it around when benchmarking but it's definitely not needed.

I don't see why it wouldn't be thread-safe though, each loop is a completely independent simple calculation and access to the arrays should be handled by locks. Maybe inbounds and threads interactions have changed recently and it somehow affects the lock?

In any case, it's fine to single thread that loop, though it's not a satisfying explanation.

TSGut commented 2 months ago

I am not convinced that's what's going wrong here... Removing @threads and @inbounds in the ClassicalOPs loop, the following still doesn't work from this PR branch

julia> A = SemiclassicalJacobi(2.0, 3.0, 4.0, 5.0);

julia> B = SemiclassicalJacobi(2.0, 4.0, 5.0, 6.0);

julia> R = A \ B; # fine

julia> R[1, 1] # still not fine

In fact I think the above code never even hits that loop, it gets stuck somewhere else?

DanielVandH commented 2 months ago

Hm.. maybe there's a couple things going on at once. I've seen errors a few times now where that function appears in the stack trace, but maybe it's a coincidence. I'll continue looking into it tomorrow then

On Sat, 15 June 2024, 10:35 pm Timon Salar Gutleb, @.***> wrote:

I am not convinced that's what's going wrong here... Removing @threads and @inbounds in the ClassicalOPs loop, the following still doesn't work from this PR branch

julia> A = SemiclassicalJacobi(2.0, 3.0, 4.0, 5.0);

julia> B = SemiclassicalJacobi(2.0, 4.0, 5.0, 6.0);

julia> R = A \ B; # fine

julia> R[1, 1] # still not fine

— Reply to this email directly, view it on GitHub https://github.com/JuliaApproximation/SemiclassicalOrthogonalPolynomials.jl/pull/109#issuecomment-2170878964, or unsubscribe https://github.com/notifications/unsubscribe-auth/AWZPH4FF6CRESHH3HLIGXXTZHSXRBAVCNFSM6AAAAABJLYIVWOVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDCNZQHA3TQOJWGQ . You are receiving this because you authored the thread.Message ID: <JuliaApproximation/SemiclassicalOrthogonalPolynomials. @.***>

TSGut commented 2 months ago

It's 100% a good idea to disable the multithreading as a first step because it makes debugging easier. Just keep in mind it may not be the real culprit and once we fix whatever else is happening with the indexing of the inverse we should double check whether we want to keep or discard the loop multithreading.

dlfivefifty commented 2 months ago

Array accessing doesn’t use locks. If it did it would extremely slow.

also to make that for loop fast there’s a lot of things you could do, eg not repeatively access the same entry in memory

TSGut commented 2 months ago

I meant writing when I said accessing (though you're right locks are manually set on that) but even without a lock, they each only ever write to independent entries and this should never have a race condition? Are you saying the following snippet isn't threadsafe in Julia? I would be quite surprised by that...

julia> function threadloop(N)
           a = zeros(N)
           b = view(ones(2*N),1:N)
           Threads.@threads for i = 1:N
               a[i] = b[i]
           end
           return a
       end
threadloop (generic function with 1 method)

julia> N = 10000000;

julia> threadloop(N) == ones(N)
true

That's exactly what that loop is, it's filling in two preallocated vectors entry by entry from views of existing arrays, completely independently.

BTW I really don't think it matters whether we multithread that loop (the Q version we actually use by default is not multithreaded, much more complex and runs perfectly fine too), the complexity of that loop is negligible if used in a serious application - I am just trying to understand whether this is actually an issue here.

DanielVandH commented 2 months ago

My concern with it was in case the getindex method is doing any computations in the background for the vectors. I haven't checked yet if that's what it actually does (it was just why I first thought about that @threads computation since I remembered it used threading somewhere). e.g. the following method can segfault easily when using threads

struct AdaptiveVector{T} <: AbstractVector{T}
    x::Vector{T}
end 
function Base.getindex(x::AdaptiveVector, i)
    if i > length(x.x)
        pl = length(x.x)
        resize!(x.x, i)
        x.x[pl+1:end] .= (pl+1):i
    end 
    return x.x[i]
end
function threadloop(N)
    a = AdaptiveVector(1:100|>collect)
    b = zeros(N)
    Threads.@threads for i in 1:N 
        b[i] = a[i]
    end
    return b
end
function loop(N)
    a = AdaptiveVector(1:100|>collect)
    b = zeros(N)
    for i in 1:N 
        b[i] = a[i]
    end
    return b
end
N = 10000
threadloop(N) == loop(N) # false, and sometimes threadloop(N) will segfault
dlfivefifty commented 2 months ago

If accessing U causes a cached array to resize then it’s not thread safe

TSGut commented 2 months ago

If accessing U causes a cached array to resize then it’s not thread safe

Definitely. That's why U is resized to the max size prior to the loop. There could however be an issue with that resize call.

DanielVandH commented 2 months ago

I am not convinced that's what's going wrong here... Removing @threads and @inbounds in the ClassicalOPs loop, the following still doesn't work from this PR branch

julia> A = SemiclassicalJacobi(2.0, 3.0, 4.0, 5.0);

julia> B = SemiclassicalJacobi(2.0, 4.0, 5.0, 6.0);

julia> R = A \ B; # fine

julia> R[1, 1] # still not fine

In fact I think the above code never even hits that loop, it gets stuck somewhere else?

This does seem to be part of the main issue. It gets stuck when trying to multiply two InvLayouts

using SemiclassicalOrthogonalPolynomials, ArrayLayouts
using SemiclassicalOrthogonalPolynomials: semijacobi_ldiv
using LazyArrays
A = SemiclassicalJacobi(2.0, 3.0, 4.0, 5.0);
B = SemiclassicalJacobi(2.0, 4.0, 5.0, 6.0);
Q, P = A, B 
Qt, Qa, Qb, Qc = Q.t, Q.a, Q.b, Q.c 
Δa, Δb, Δc = Qa - P.a, Qb - P.b, Qc - P.c
R = ApplyArray(inv, semijacobi_ldiv(P.t, P.a, P.b, P.c, SemiclassicalJacobi(P.t, Qa, P.b, P.c)))
ApplyArray(*, R, R) # can't print (i.e. can't index)

Looking into it. The InvLayouts isn't the issue exactly, since e.g.

X = [1.0 2.0; 3.0 4.0]
Y = ApplyArray(inv, ApplyArray(*, X, [1.0 0.0; 0.0 1.0]))
Y * Y

works but MemoryLayout(Y) == MemoryLayout(R).

dlfivefifty commented 2 months ago

I tried to get the tests working but realised it was very slow because of some decisions in the design

DanielVandH commented 2 months ago

Is that because of me not using the constructor SemiClassicalJacobi(t, a, b, c, P), or is my idea for decrementing parameters one at a time slow?