JuliaMath / JuliaMath/QuadGK.jl

infinite limits with units and segbuf broken

Open
#95 0 comments 0 reactions 0 assignees View on GitHub
Dominant language
Julia
Stars
314
Forks
41
PR merge metrics
No merged PRs in 30d

Description

The issue is that segbufs don't always work with infinite limits, such as when the limits have units.
This is a MWE
```julia
using QuadGK, Unitful

Tu = typeof(1.0u"m")
segbuf = QuadGK.alloc_segbuf(Tu, Tu, Tu)

quadgk(x -> 1/(1+oneunit(x)/x), 0.0u"m", 1.0u"m", segbuf=segbuf) # OK (0.3068528194400547 m, 1.9931833961095435e-11 m)

quadgk(x -> 1/(1+oneunit(x)/x), 0.0u"m", Inf*u"m", segbuf=segbuf) # broken
```
which gives the error
```julia
ERROR: DimensionError: m and 0.0 are not dimensionally compatible.
Stacktrace:
[1] convert(#unused#::Type{Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}}, x::Float64)
@ Unitful ~/.julia/packages/Unitful/jiUeO/src/conversion.jl:106
[2] QuadGK.Segment{Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}, Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}, Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}}(a::Float64, b::Float64, I::Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}, E::Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}})
@ QuadGK ~/.julia/dev/QuadGK/src/evalrule.jl:3
[3] convert(#unused#::Type{QuadGK.Segment{Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}, Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}, Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}}}, s::QuadGK.Segment{Float64, Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}, Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}})
@ QuadGK ~/.julia/dev/QuadGK/src/evalrule.jl:10
[4] setindex!(A::Vector{QuadGK.Segment{Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}, Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}, Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}}}, x::QuadGK.Segment{Float64, Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}, Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}}, i1::Int64)
@ Base ./array.jl:969
[5] macro expansion
@ ./broadcast.jl:974 [inlined]
[6] macro expansion
@ ./simdloop.jl:77 [inlined]
[7] copyto!
@ ./broadcast.jl:973 [inlined]
[8] copyto!
@ ./broadcast.jl:926 [inlined]
[9] materialize!
@ ./broadcast.jl:884 [inlined]
[10] materialize!
@ ./broadcast.jl:881 [inlined]
[11] do_quadgk(f::QuadGK.var"#18#27"{var"#11#12", Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}, Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}}, s::Tuple{Float64, Float64}, n::Int64, atol::Nothing, rtol::Nothing, maxevals::Int64, nrm::typeof(LinearAlgebra.norm), segbuf::Vector{QuadGK.Segment{Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}, Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}, Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}}})
@ QuadGK ~/.julia/dev/QuadGK/src/adapt.jl:39
[12] #50
@ ~/.julia/dev/QuadGK/src/adapt.jl:235 [inlined]
[13] handle_infinities(workfunc::QuadGK.var"#50#51"{Nothing, Nothing, Int64, Int64, typeof(LinearAlgebra.norm), Vector{QuadGK.Segment{Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}, Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}, Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}}}}, f::var"#11#12", s::Tuple{Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}, Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}})
@ QuadGK ~/.julia/dev/QuadGK/src/adapt.jl:126
[14] quadgk(::Function, ::Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}, ::Vararg{Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}}; atol::Nothing, rtol::Nothing, maxevals::Int64, order::Int64, norm::Function, segbuf::Vector{QuadGK.Segment{Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}, Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}, Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}}})
@ QuadGK ~/.julia/dev/QuadGK/src/adapt.jl:234
[15] top-level scope
@ REPL[11]:1
```

Contributor guide

No contributing guide indexed for this repository

Assessment

This issue has not been assessed yet.

Get new issues in your inbox

A short digest of beginner-friendly GitHub issues.