sciml / muladdmacro.jl Goto Github PK

@JuliaRegistrator register()

This issue is used to trigger TagBot; feel free to unsubscribe.

If you haven't already, you should update your TagBot.yml to include issue comment triggers.
Please see this post on Discourse for instructions and more details.

If you'd like for me to do this for you, comment TagBot fix on this issue.
I'll open a PR within a few hours, please be patient!

I expected this to work, but it didn't. There is generally no problem with having minus signs inside muladd, because FMA instructions allow for that, so expressions like b-a*x should be converted to muladd(-a, x, b).

julia> macroexpand(:(MuladdMacro.@muladd b - a * x))
:(b - a * x)

julia> macroexpand(:(MuladdMacro.@muladd b + a * x))
:((muladd)(a, x, b))

┌ Warning: `a::Number + b::AbstractArray` is deprecated, use `a .+b` instead.
│   caller = muladd at math.jl:1014 [inlined]

is showing up in some tests. That's caused by

julia> muladd(1, 1, ones(1))
┌ Warning: `a::Number + b::AbstractArray` is deprecated, use `a .+ b` instead.
│   caller = muladd(::Int64, ::Int64, ::Array{Float64,1}) at math.jl:1014
└ @ Base.Math ./math.jl:1014

@devmotion does the newest version not catch broadcasting as muladd.(...)? This is showing up in OrdinaryDiffEq.jl and StochasticDiffEq.jl tests.

@muladd should not change the order in which operands are multiplied and added, since any reordering might introduce numerical instabilities. So the examples of the README should rather give

julia> macroexpand(:(@muladd k3 = f(t + c3*dt, @. uprev+dt*(a031*k1+a032*k2))))
:(k3 = f((muladd)(c3, dt, t), (muladd).(dt, (muladd).(a032, k2, *.(a031, k1)), uprev)))

julia> macroexpand(:(@muladd integrator.EEst = integrator.opts.internalnorm((update - dt*(bhat1*k1 + bhat4*k4 + bhat5*k5 + bhat6*k6 + bhat7*k7 + bhat10*k10))./ @. (integrator.opts.abstol+max(abs(uprev),abs(u))*integrator.opts.reltol))))
:(integrator.EEst = integrator.opts.internalnorm((update - dt * (muladd)(bhat10, k10, (muladd)(bhat7, k7, (muladd)(bhat6, k6, (muladd)(bhat5, k5, (muladd)(bhat4, k4, bhat1 * k1)))))) ./ (muladd).(max.(abs.(uprev), abs.(u)), integrator.opts.reltol, integrator.opts.abstol)))

This occured first in SciML/OrdinaryDiffEq.jl#93 (comment)

I would suggest changing @muladd in such a way that at least the order of evaluation of the summation-multiplication complexes is left unchanged (i.e. the implementation should be reversed) and summands with n factors are split in two groups of the first n-1 factors and the last factor (i.e. contrary to the current implementation; this should leave the evaluation order unchanged also inside of each summand). However, I'm not sure whether it would be best to evaluate the sum of all not-multiplication summands first (as currently done) or last.

MacroTools' pattern matching is very slow, mostly due to using exceptions to signal no-match. This causes JuliaLang/julia#28221. Ideally MacroTools can be fixed too, but in the meantime it seems easier to switch this package to use simple tests like isexpr(ex, :call) && ex.args[1] == :+ instead of MacroTools. Of course that can be abstracted into iscall(ex, :+), etc. --- feel free to go nuts with that.