Refactor Linear Algebra's matmul to use a safer scratch buffer

[danielwe]

https://github.com/JuliaLang/julia/blob/54adf2193a6ea1b10cd1a184a7ba863aa985eede/stdlib/LinearAlgebra/src/matmul.jl#L536
https://github.com/JuliaLang/julia/blob/54adf2193a6ea1b10cd1a184a7ba863aa985eede/stdlib/LinearAlgebra/src/matmul.jl#L558

I stumbled upon these two comments while perusing the base linalg code. Related statement in JuliaLang/julia#23914, where they were introduced:

Also marked a few even more invalid places where managed pointer is passed to unsafe_wrap (the fix for those will be much more involved).

Not terribly reassuring seeing this in production code, so I figured it deserved an open issue.

[mbauman]

Discussed a few days ago on the #gripes slack channel. Marked in JuliaLang/julia#23914, invalid because "it’s pointing at an unaligned array though, and probably expecting to use simd instructions on it"

[mbauman]

Could you please explain the close, @yuyichao? I think this issue is useful to track the need to change this code.

[StefanKarpinski]

Reopening until we have an explanation or a fix.

[yuyichao]

The Base code is full of these. If you want to make a list of all the files to scan for any invalid use it's fine. Otherwise, this is just an issue with a uninformative title that can never be closed.

Also alignment isn't the issue the issue is explained in the original issue as

managed pointer is passed to unsafe_wrap

quoted above.

[mbauman]

In this specific case, though, the managed pointer is bound to a const global — it will always exist (edit: ah, unless it gets reshaped, I suppose). So (barring that) is this specific case really an issue if there's no alignment problem?

[StefanKarpinski]

Just because we don't currently have a complete set of issues for broken things with comments does not mean that we should not try to have issues for them. Comments tend not to get fixed or found as effectively as issues in an issue tracker. This is why we have issue trackers instead of just leaving comments all over the code. The issue tracker also has the benefit of supporting discussions like this one which help resolve and/or clarify the problem. In short, lack of a complete set of issues is not a justification for having no issue in this case.

[yuyichao]

Conceptually this is never allowed. Practically there are many compiler limitations and that's exactly why we have this issue in the first place.
Practically, I don't think this will be too bad with current compiler though future assumptions on aliasing could make this very bad. There is another issue with this code though that's unrelated with GC, it's not thread safe/reentrant.

lack of a complete set of issues is not a justification for having no issue in this case.

No. There has been many vague issues like this that are just left open until it's hard to tell if the issue is fixed or not. That's also why a uninformative title doesn't help since people will be opening duplicates and that go directly against keeping track of discussion. I don't want to rename it since I don't know what's the expected scope. It's a very broad issue but the issue for this specific case is also somewhat special as mentioned above. Now it's renamed and if people are actually interested in fixing it now, I have no problem keeping it open.

[danielwe]

Thanks for the rename, I lack the chops to fully understand the problem and write an appropriate title myself.

Based on a quick search there aren't that many fixme comments, and these two stand out as both particularly dramatic and opaque, and right in the middle of a core number crunching routine, so an issue seemed appropriate -- if only to rewrite them as less disconcerting and more informative. No intentions of spamming your issue tracker. https://github.com/JuliaLang/julia/search?p=2&q=fixme&type=&utf8=%E2%9C%93

[yuyichao]

there aren't that many fixme comments

Yes, I gave up after adding a few...

[thchr]

Figured it might be helpful to explicitly crosspost a likely consequence of this (noted by @mateuszbaran in JuliaArrays/StaticArrays.jl#966):

using StaticArrays
A = @SMatrix [1 2; 3 4]
blockA = fill(A, 1, 1)
blockA * transpose(blockA) - blockA * collect(transpose(blockA))

gives [2 -1; 4 -3] instead of [0 0; 0 0].

It doesn't seem like a StaticArrays issue (can reproduce for other simple isbits matrix type).

[KristofferC]

These should have been removed recently so I think this can be closed. Cc @tkf

[thchr]

Ah, cool! Should/will JuliaLang/julia#42309 be backported then?

[tkf]

I didn't know about this issue. Yeah, I think this can be closed.

Should/will JuliaLang/julia#42309 be backported then?

Yeah, I'd say JuliaLang/julia#42309 is a pretty safe patch to backport.

[thchr]

@maltezfaria pointed out that the test case I linked above actually still fails in exactly the same way on current master.
I tested this on 1.8.0-DEV.763 (commit cd19e97771) and can confirm. It's then not exactly clear if JuliaLang/julia#42309 actually fixed all issues.

[tkf]

Does blockA hit the tiling case? Isn't it too small?

[KristofferC]

2 and 3 are too small. 1 is big enough ;)

My guess is that the transpose is not applied recursively, and there is a bug somewhere here (maybe just missing a tranpose call): https://github.com/JuliaLang/julia/blob/df81a0d33d81e54c662932fd482cab9c67a09e39/stdlib/LinearAlgebra/src/transpose.jl#L197

julia> blockA
1×1 Matrix{SMatrix{2, 2, Int64, 4}}:
 [1 2; 3 4]

julia> blockB = similar(blockA);

julia> LinearAlgebra.copy_transpose!(blockB, 1:1, 1:1, blockA, 1:1, 1:1);

julia> blockB # not transposed
1×1 Matrix{SMatrix{2, 2, Int64, 4}}:
 [1 2; 3 4]

julia> BlockA = [1 3; 2 4];

julia> BlockB = similar(BlockA);

julia> LinearAlgebra.copy_transpose!(BlockB, 1:2, 1:2, BlockA, 1:2, 1:2);

julia> BlockB # transposed
2×2 Matrix{Int64}:
 1  2
 3  4