`matmul!` with `ReshapedArray` of an `Adjoint` does not work

[nathanaelbosch]

Is there any way to get this to work? Example code:

  X = reshape(rand(2, 3)', 2, 3)
  N = rand(2, 2)
  V = rand(2, 3)
  typeof(X) # Base.ReshapedArray{Float64, 2, Adjoint{Float64, Matrix{Float64}}, Tuple{Base.MultiplicativeInverses.SignedMultiplicativeInverse{Int64}}}
  typeof(N) # Matrix{Float64} (alias for Array{Float64, 2})
  typeof(V) # Matrix{Float64} (alias for Array{Float64, 2})
  matmul!(X, N, V)

Errors with:

ERROR: MethodError: no method matching matmul!(::Base.ReshapedArray{Float64, 2, …}, ::Matrix{Float64}, ::Matrix{Float64}, ::StaticInt{1}, ::StaticInt{0}, ::Nothing, ::Nothing, ::Nothing)

Closest candidates are:
  matmul!(::AbstractVecOrMat, ::AbstractMatrix, ::AbstractVecOrMat, ::Any, ::Any, ::Any, ::Nothing, ::StaticInt{2})
   @ Octavian ~/.julia/packages/Octavian/iOEo8/src/matmul.jl:484
  matmul!(::AbstractVecOrMat, ::AbstractMatrix, ::AbstractVecOrMat, ::Any, ::Any, ::Any, ::Any, ::StaticInt)
   @ Octavian ~/.julia/packages/Octavian/iOEo8/src/matmul.jl:510
  matmul!(::AbstractVecOrMat, ::AbstractMatrix, ::AbstractVecOrMat, ::Any, ::Any, ::Any)
   @ Octavian ~/.julia/packages/Octavian/iOEo8/src/matmul.jl:474
  ...

Stacktrace:
 [1] matmul! @ Octavian ~/.julia/packages/Octavian/iOEo8/src/matmul.jl:437

[chriselrod]

julia> X = reshape(reshape(1:6,2,3)', 2, 3)
2×3 reshape(adjoint(reshape(::UnitRange{Int64}, 2, 3)), 2, 3) with eltype Int64:
 1  5  4
 3  2  6

That is an unusual memory layout. It'd be easy to make this work for A or B in matmul!(C, A, B), because Octavian already has temporary buffers for A and B that it copies data over into in order to control the memory layout. But it doesn't do something like that for C.

You're welcome to make a PR to add support. I can answer questions.

[nathanaelbosch]

That is an unusual memory layout.

It comes up in a mul!(X::Adjoint, K::KroneckerProduct, Y::AbstractMatrix context (roughly), and the Kronecker-trick needs the reshapes. I guess I might circumvent X being an adjoint and replace it with a standard matrix, but that makes things look quite a lot less natural in the code. So from my perspective it would be great to have this!

By the way, mul! is also quite inefficient for this, it seems to fall back to LinearAlgebra.generic_matmatmul! if I'm not mistaken.

You're welcome to make a PR to add support. I can answer questions.

I'm not so sure if I'm the right person for this, I am definitely not familiar with the codebase and I have no idea how involved it is to get this done. But depending on the feasibility and required effort I could try.