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')

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]

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