add a missing transpose in LinearAlgebra.copy_transpose! to handle matmul of isbits array-elements with non-identity transpose

[KristofferC]

Fixes the issue reported in #26008 (comment)

julia> using StaticArrays

julia> A = @SMatrix [1 2; 3 4]
2×2 SMatrix{2, 2, Int64, 4} with indices SOneTo(2)×SOneTo(2):
 1  2
 3  4

julia> blockA = fill(A, 1, 1)
1×1 Matrix{SMatrix{2, 2, Int64, 4}}:
 [1 2; 3 4]

julia> blockA * transpose(blockA) - blockA * collect(transpose(blockA))
1×1 Matrix{SMatrix{2, 2, Int64, 4}}:
 [0 0; 0 0]

Should probably add a test but it is a bit it requires quite a bit of scaffolding to set up. Maybe someone has that already available? cc @thchr

[thchr]

Cool - thanks for fixing this!

It is a little tedious to construct a reproducible example without bringing in StaticArrays, indeed, but here's a reasonably minimal and self-contained example:

using Test
using LinearAlgebra: Transpose

struct SMatrix2x2 <: AbstractMatrix{Int}
    cols::NTuple{2, NTuple{2, Int}} # tuples of columns
end
Base.size(::SMatrix2x2) = (2,2)
function Base.getindex(A::SMatrix2x2, i::Int, j::Int)
    @boundscheck (1 ≤ i ≤ 2 && 1 ≤ j ≤ 2) || throw(BoundsError(A, (i,j)))
    return A.cols[j][i]
end
Base.IndexStyle(::Type{SMatrix2x2}) = IndexCartesian()

function Base.convert(::Type{SMatrix2x2}, A::Matrix{Int})
    size(A) == (2,2) || throw(DimensionMismatch("A must be a 2x2 matrix"))
    SMatrix2x2(((A[1,1], A[2,1]), (A[1,2], A[2,2])))
end
function Base.:*(A::SMatrix2x2, B::Union{SMatrix2x2, Transpose{Int, SMatrix2x2}})
    SMatrix2x2(((A[1,1]*B[1,1] + A[1,2]*B[2,1], A[2,1]*B[1,1] + A[2,2]*B[2,1] ),
                (A[1,1]*B[1,2] + A[1,2]*B[2,2], A[2,1]*B[1,2] + A[2,2]*B[2,2])))
end
function Base.:+(A::SMatrix2x2, B::Union{SMatrix2x2, Transpose{Int, SMatrix2x2}})
    SMatrix2x2(((A[1,1]+B[1,1], A[2,1]*B[2,1] ),
                (A[1,2]+B[1,2], A[2,2]*B[2,2])))
end

A = SMatrix2x2(((1,3),(2,4)))
bA = fill(A, 1, 1)
@test bA*transpose(bA) == bA*collect(bA')

The extensions of +, * are needed e.g. for a promote_op(matprod, ...) call at one point. The convert method is needed during mul!.

EDIT: This fails in the same way as SMatrix{2,2,Int} on current master but still fails on this PR; I'll see if something else is needed.

[thchr]

Ah, this change indeed makes blockA * transpose(blockA) == blockA * collect(transpose(blockA)) but neither of the variants are actually right now, cf:

julia> A, B = (@SMatrix [1 2; 3 4]), [1 2; 3 4]
julia> blockA, blockB = fill(A, 1, 1), fill(B, 1, 1)

julia> blockA * transpose(blockA)
1×1 Matrix{SMatrix{2, 2, Int64, 4}}:
 [7 15; 10 22]
julia> blockB * transpose(blockB)
1×1 Matrix{Matrix{Int64}}:
 [5 11; 11 25]

Without this PR, the blockA * collect(transpose(blockA)) variant is correct at least.

[KristofferC]

Hm, maybe some more is needed. Thanks for the example!

[N5N3]

Looks like _generic_matmatmul! call copy_transpose! to make A in A * B' dense in column, where the added transpose is not wanted.
Also notice than blockA * blockA' still throw errors even we handle inner transpose correctly, as conj! is not defined for SMatrix eltype.
So making non-number eltype fallback to the the naive algorithm should be the simplest fix.

Edit: looks like copy_transpose! is undocumented, so we can change it's api with no breakage, the following code should fix all above examples and blockA * blockA'

function copy_transpose!(B::AbstractVecOrMat, ir_dest::AbstractRange{Int}, jr_dest::AbstractRange{Int},
                         A::AbstractVecOrMat, ir_src::AbstractRange{Int}, jr_src::AbstractRange{Int}; op = transpose)
    if length(ir_dest) != length(jr_src)
        throw(ArgumentError(string("source and destination must have same size (got ",
                                   length(jr_src)," and ",length(ir_dest),")")))
    end
    if length(jr_dest) != length(ir_src)
        throw(ArgumentError(string("source and destination must have same size (got ",
                                   length(ir_src)," and ",length(jr_dest),")")))
    end
    @boundscheck checkbounds(B, ir_dest, jr_dest)
    @boundscheck checkbounds(A, ir_src, jr_src)
    idest = first(ir_dest)
    @inbounds for jsrc in jr_src
        jdest = first(jr_dest)
        for isrc in ir_src
            B[idest,jdest] = op(A[isrc,jsrc])
            jdest += step(jr_dest)
        end
        idest += step(ir_dest)
    end
    return B
end

function copyto!(B::AbstractVecOrMat, ir_dest::AbstractUnitRange{Int}, jr_dest::AbstractUnitRange{Int}, tM::AbstractChar, M::AbstractVecOrMat, ir_src::AbstractUnitRange{Int}, jr_src::AbstractUnitRange{Int})
    if tM == 'N'
        copyto!(B, ir_dest, jr_dest, M, ir_src, jr_src)
    elseif tM == 'T'
        LinearAlgebra.copy_transpose!(B, ir_dest, jr_dest, M, jr_src, ir_src)
    else
        LinearAlgebra.copy_transpose!(B, ir_dest, jr_dest, M, jr_src, ir_src; op = adjoint)
    end
    B
end

function copy_transpose!(B::AbstractMatrix, ir_dest::AbstractUnitRange{Int}, jr_dest::AbstractUnitRange{Int}, tM::AbstractChar, M::AbstractVecOrMat, ir_src::AbstractUnitRange{Int}, jr_src::AbstractUnitRange{Int})
    if tM == 'N'
        LinearAlgebra.copy_transpose!(B, ir_dest, jr_dest, M, ir_src, jr_src; op = identity)
    elseif tM == 'T'
        copyto!(B, ir_dest, jr_dest, M, jr_src, ir_src)
    else
        @views B[ir_dest, jr_dest] .= adjoint.(M[jr_src, ir_src])
    end
    B
end

[dkarrasch]

All previously failing test cases pass now, so I think this can be closed.