-
Notifications
You must be signed in to change notification settings - Fork 3
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
Conversation
Try running with |
With that option I now get no segfault but I also get no 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 |
Not 100% sure but it seems the segfault goes away when I remove |
Hm, the resizing prior to the loop should make that safe if the sizes are correctly chosen, is this with |
Without 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 |
Why does this for loop use |
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. |
I am not convinced that's what's going wrong here... Removing 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? |
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
<#109 (comment)>,
or unsubscribe
<https://github.com/notifications/unsubscribe-auth/AWZPH4FF6CRESHH3HLIGXXTZHSXRBAVCNFSM6AAAAABJLYIVWOVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDCNZQHA3TQOJWGQ>
.
You are receiving this because you authored the thread.Message ID:
<JuliaApproximation/SemiclassicalOrthogonalPolynomials.
***@***.***>
|
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. |
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 |
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. |
My concern with it was in case the 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 |
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. |
This does seem to be part of the main issue. It gets stuck when trying to multiply two 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 |
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))) |
There was a problem hiding this comment.
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.
There was a problem hiding this comment.
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)
?
There was a problem hiding this comment.
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.
There was a problem hiding this comment.
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)
There was a problem hiding this comment.
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)
.
There was a problem hiding this comment.
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.
There was a problem hiding this comment.
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.
There was a problem hiding this comment.
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?
There was a problem hiding this comment.
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
There was a problem hiding this comment.
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
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) |
There was a problem hiding this comment.
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 =
?
There was a problem hiding this comment.
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).
There was a problem hiding this comment.
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
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I'll remove it
I tried to get the tests working but realised it was very slow because of some decisions in the design |
Is that because of me not using the constructor |
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))) |
There was a problem hiding this comment.
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))) |
There was a problem hiding this comment.
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)
.
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
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.