Multiplication involving transpose of a `Matrix{<:SMatrix}`

[maltezfaria]

Hi,

I ran into the following (wrong?) behaviour when trying to multiply the following matrices of SMatrix:

using StaticArrays

A = @SMatrix rand(3, 3)
B = @SMatrix rand(3, 3)

blockA = fill(A, 3, 2)
blockB = fill(B, 3, 2)

val1 = blockA * transpose(blockB)
val2 = blockA * collect(transpose(blockB))

val1 ≈ val2 # false, by a lot

On my machine, the difference goes away if either A or B is changed to an MMatrix, which I found a bit puzzling. Here is some info on my version:

Julia Version 1.6.2
Commit 1b93d53fc4 (2021-07-14 15:36 UTC)
Platform Info:
  OS: macOS (x86_64-apple-darwin18.7.0)
  CPU: Intel(R) Core(TM) i9-9880H CPU @ 2.30GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-11.0.1 (ORCJIT, skylake)
Environment:
  JULIA_EDITOR = code
  JULIA_NUM_THREADS = 8

Thanks!

[mateuszbaran]

I think this may be an issue with Julia matmul code here: https://github.com/JuliaLang/julia/blob/1b93d53fc4bb59350ada898038ed4de2994cce33/stdlib/LinearAlgebra/src/matmul.jl#L831 . SMatrix is an isbits type so it tries to do that tiled thing and fails.

Ref JuliaLang/LinearAlgebra.jl#500 .

[thchr]

I guess this could be closed cf. the discussion in JuliaLang/LinearAlgebra.jl#500.

[maltezfaria]

Thanks for tracking this. Should the problem be fixed with JuliaLang/julia#42309? I tried running the example above on the nightly build, and I still get the same behaviour.

[maltezfaria]

I ran into another strange behaviour as well when trying to rely on the generic matmul for some operations on matrices of static matrices which I think should make sense. The example below illustrates the issue:

using StaticArrays
using LinearAlgebra

smat = @SMatrix [1 2; 3 4]

# mul! with matrices works fine
A = fill(smat, 1, 1)
B = fill(smat, 1, 1)
C = fill(smat, 1, 1)

mul!(copy(C),A,B,1,0) ≈ A*B # ok
mul!(copy(C),A,B,1,1) ≈ A*B + C # ok

# mul! with `matvecmul` throws an error
A = fill(smat, 1, 1)
B = fill(smat, 1)
C = fill(smat, 1)

mul!(copy(C),A,B,1,0) ≈ A*B # ok
mul!(copy(C),A,B,1,1) ≈ A*B + C # error with MulAddAmul

Thanks for the help!

[thchr]

Thanks for tracking this. Should the problem be fixed with JuliaLang/julia#42309? I tried running the example above on the nightly build, and I still get the same behaviour.

Huh, yeah, I confirm that this is indeed not fixed after all. Testing on 1.8.0-DEV.763, the previous test case still fails.
I suppose this should be reopened then :(

[maltezfaria]

Fixed in JuliaLang/LinearAlgebra.jl#1033