Wrong evaluation in sparse/higherorderfns.jl fix#23857

[KlausC]

Division of a sparse matrix by row vector causes 0's to become NaN
#23857

[Sacha0]

Thanks @KlausC! Apologies in advance for not reviewing this pull request in the short term. (Time is short with feature freeze fast approaching, and I have been remiss in reviewing a number of higher priority items as it stands.) Best!

[KlausC]

@Sacha0 just a reminder - it's 30 days now

[Sacha0]

Regrettably, work falling under the heading bugfix will need to wait till after proper feature freeze (and likely then some) :).

[KlausC]

Maybe it is time now to have a lock on this bug fix. The change was rather simple, as far as I remember.

[KlausC]

@mbauman , @andreasnoack , wouldn't it make sense to merge this minor bug fix as well. I fear it is getting lost after similar changes have been done to the same file(s). #26474

[mbauman]

I'm so sorry that this has fallen through the cracks. The change here is indeed simpler than it looks — here's a better view: https://github.com/JuliaLang/julia/pull/24806/files?w=1

I'm hesitant here because a very similar structure is used in the other branches and functions here — do they need the same fix? I'd like to understand a bit more about why this is failing, but I'm not seeing it at the moment.

[KlausC]

I'd like to understand a bit more about why this is failing

It is nearly 4 month by now. I think, I understood it then, but have to recover now. I will try tomorrow, to explain the change. At least it did not break any tests and fixed the reported errors.

[Sacha0]

Having briefly refreshed my memory of these methods, I believe that the _iszero(fzAvB) check (and similar _iszero(fvAzB) check in the preceding branch where A rather than B expands vertically) should instead be something to the effect of fzAvB == fillvalue in the methods for the not-zero-preserving-f case. These branches should stay around, however, as otherwise performance may degrade appreciably. Regrettably a nanosoldier suite for these methods remains a "one day" project. Best!

[KlausC]

_broadcast_nonzeropres! is only called, if A or B have to be expanded vertically, that means each column of A or B has to be replaced by copies of the first element.

I think the _broadcast_nozeropres! function works as follows:

1. It expands destination array C to full size (C.m*C.n) and fills all value cells with value f(zero(), zero()). (lines 644ff)
2. According to the shapes of A and B it uses different "branches", which have some structure in common
3. The outer loop proceeds column by column - index j.
4. assuming A is to be expanded, Ax is the value of A[1,j], which has to be joined with all B[i,j] for all row indices i.
5. An optimization is tried as follows: if f(Ax, zero(()) == fillvalue only the elements B[i,j] different from zero can contribute new values to the result. So only nonzero part of B[:,j] has to be looped, thus exploiting the sparsity of B . The other slots of C already are populated with fillvalue.
6. In the other case, f(Ax, B[i,j]) has to be stored for all i. In this case, the pre-calculated f(Ax, zero()) can be re-used for all i with B[i,j] == zero().

The bug is, that instead of f(Ax, zero(()) == fillvalue the program checks for f(Ax, zero()) == zero(). That is presumably a copy-error form _broadcst_zeropres!.

[mbauman]

Thanks for the explanation and change. Makes sense, and this now matches the vararg _broadcast_notzeropres! method below.

[Sacha0]

Thanks for persisting with this @KlausC! :)

[Sacha0]

(Absent objections, requests for time, or someone beating me to it, I plan to merge these changes tomorrow morning PST or later. Best!)