Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Support computing connection matrices when decrementing parameters #109

Closed
wants to merge 3 commits into from

Conversation

DanielVandH
Copy link
Member

@DanielVandH DanielVandH commented Jun 15, 2024

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
Copy link
Member

Try running with julia --check-bounds=yes

@DanielVandH
Copy link
Member Author

DanielVandH commented Jun 15, 2024

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
Copy link
Member Author

@TSGut
Copy link
Member

TSGut commented Jun 15, 2024

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

@DanielVandH
Copy link
Member Author

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
Copy link
Member

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

@TSGut
Copy link
Member

TSGut commented Jun 15, 2024

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
Copy link
Member

TSGut commented Jun 15, 2024

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
Copy link
Member Author

DanielVandH commented Jun 15, 2024 via email

@TSGut
Copy link
Member

TSGut commented Jun 15, 2024

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
Copy link
Member

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
Copy link
Member

TSGut commented Jun 16, 2024

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
Copy link
Member Author

DanielVandH commented Jun 16, 2024

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
Copy link
Member

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

@TSGut
Copy link
Member

TSGut commented Jun 16, 2024

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
Copy link
Member Author

DanielVandH commented Jun 16, 2024

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).

This was referenced Jun 16, 2024
if isinteger(Δa) && isinteger(Δb) && isinteger(Δc) # (Δa,Δb,Δc) are integers -> use QR/Cholesky iteratively
if ((isone(Δa)||isone(Δa/2)) && iszero(Δb) && iszero(Δc)) || (iszero(Δa) && (isone(Δb)||isone(Δb/2)) && iszero(Δc)) || (iszero(Δa) && iszero(Δb) && (isone(Δc)||isone(Δc/2)))
if Δa < 0
R = ApplyArray(inv, semijacobi_ldiv(P.t, P.a, P.b, P.c, SemiclassicalJacobi(P.t, Qa, P.b, P.c)))
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The constructor SemiclassicalJacobi(P.t, Qa, P.b, P.c) should be avoided at all costs: we want to make sure we are building the new OP from the existing one as otherwise we have to redo an expensive Cholesky.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Should they all just be SemiclassicalJacobi(t, a, b, c, P)?

Copy link
Member

@TSGut TSGut Jun 21, 2024

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If Qa < P.a then SemiclassicalJacobi(P.t, Qa, P.b, P.c, P) won't work because lowering isn't implemented.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Then shouldn't this be doing it the other way around and not calling inv?

Btw the extra function semijacobi_ldiv(Qt, Qa, Qb, Qc, P::SemiclassicalJacobi) is confusing. I think there should just be the one semijacobi_ldiv(Q ::SemiclassicalJacobi, P ::SemiclassicalJacobi)

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Actually scratch that. I think it should be SemiclassicalJacobi(P.t, Qa, P.b, P.c, Q).

Copy link
Member

@TSGut TSGut Jun 21, 2024

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Just in case there is confusion as to the convention: Conceptually, semijacobi_ldiv_direct goes the whole distance in one decomposition step, so for stability reasons should only be used for changes of +1 in a parameter but in principle it can be used for any positive polynomial modification. semijacobi_ldiv iterates over semijacobi_ldiv_direct to modify the parameters up in +1 (/ -1 once implemented) increments and combines them, which is stable.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Would it make more sense then to just try and get lowering working first, then? I've not looked much into the details of it but maybe it makes things simpler.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

You should modify semijacobi_ldiv_direct with a case (-1,0,0) and so on, then call them in the same style it's currently done in semijacobi_ldiv, so that it iterates by calling itself recursively until all changes are accounted for.

To do the (-1,0,0) case it's totally fine to do the (1,0,0) case, i.e. the current code (note here the M are not optional)

M = cholesky(P.X).U
return ApplyArray(*, Diagonal(Fill(1/M[1],∞)), M)

and then instead return the inv of that. I think that is what you have been working towards?

Copy link
Member

@TSGut TSGut Jun 21, 2024

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If the inv plug and play there doesn't work (as an ApplyArray) then we can write a custom struct here to make it work - it would be nicer tho if inv of an UpperTriangular just worked there, perhaps we need to re-wrap with UpperTriangular or maybe just apply inv to M = cholesky(P.X).U (in this case and as appropriate in the others) and then undo the scaling by the diagonal manually by inverse rules.

But once that mod in semijacobi_ldiv_direct is done it should all just work automatically. I could do a working prototype pretty quickly if you guys really need lowering

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ah okay, yes I think that sounds like what I'm working towards. I didn't fully understand the semijacobi_ldiv_direct design - I should read through it and the details in the paper again. A working prototype would be pretty useful if you can, I'd appreciate that. We need lowering for my $b = -1$ work. I'd work at it now but I need to get my computer in for repairs sometime in the week before I can use Julia again.

it would be nicer tho if inv of an UpperTriangular just worked there, perhaps we need to re-wrap with UpperTriangular, I don't know if the structure wrapper is maintained.

I've been looking at inv of an UpperTriangular while trying to get all this working that would make a lot of the issues I've been having here simpler. Need to have something like JuliaArrays/LazyArrays.jl#320 resolved first (I think).

if ((isone(Δa)||isone(Δa/2)) && iszero(Δb) && iszero(Δc)) || (iszero(Δa) && (isone(Δb)||isone(Δb/2)) && iszero(Δc)) || (iszero(Δa) && iszero(Δb) && (isone(Δc)||isone(Δc/2)))
if Δa < 0
R = ApplyArray(inv, semijacobi_ldiv(P.t, P.a, P.b, P.c, SemiclassicalJacobi(P.t, Qa, P.b, P.c)))
M = ApplyArray(*, semijacobi_ldiv(Qt, Qa, Qb, Qc, SemiclassicalJacobi(P.t, Qa, P.b, P.c)), R)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why the need for M = ?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Not sure what you mean. Those M = were already there I think. I split up M and R into two lines to make it easier to read and make it easier to see that the ldiv is being applied in steps (decrementing one parameter, then onto the next).

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

But isn't the M = line just returned? So the variable is never used

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'll remove it

@dlfivefifty
Copy link
Member

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

@DanielVandH
Copy link
Member Author

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?

if isinteger(Δa) && isinteger(Δb) && isinteger(Δc) # (Δa,Δb,Δc) are integers -> use QR/Cholesky iteratively
if ((isone(Δa)||isone(Δa/2)) && iszero(Δb) && iszero(Δc)) || (iszero(Δa) && (isone(Δb)||isone(Δb/2)) && iszero(Δc)) || (iszero(Δa) && iszero(Δb) && (isone(Δc)||isone(Δc/2)))
if Δa < 0
R = ApplyArray(inv, semijacobi_ldiv(P.t, P.a, P.b, P.c, SemiclassicalJacobi(P.t, Qa, P.b, P.c)))
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Then shouldn't this be doing it the other way around and not calling inv?

Btw the extra function semijacobi_ldiv(Qt, Qa, Qb, Qc, P::SemiclassicalJacobi) is confusing. I think there should just be the one semijacobi_ldiv(Q ::SemiclassicalJacobi, P ::SemiclassicalJacobi)

if isinteger(Δa) && isinteger(Δb) && isinteger(Δc) # (Δa,Δb,Δc) are integers -> use QR/Cholesky iteratively
if ((isone(Δa)||isone(Δa/2)) && iszero(Δb) && iszero(Δc)) || (iszero(Δa) && (isone(Δb)||isone(Δb/2)) && iszero(Δc)) || (iszero(Δa) && iszero(Δb) && (isone(Δc)||isone(Δc/2)))
if Δa < 0
R = ApplyArray(inv, semijacobi_ldiv(P.t, P.a, P.b, P.c, SemiclassicalJacobi(P.t, Qa, P.b, P.c)))
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Actually scratch that. I think it should be SemiclassicalJacobi(P.t, Qa, P.b, P.c, Q).

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

3 participants