Fix & test * with matrix-of-adjoint

src/matrix_multiply_add.jl

@@ -248,6 +248,38 @@
     return _mul!(TSize(dest), mul_parent(dest), Size(A), Size(B), A, B, NoMulAdd{TMul, TDest}())
 end
 
+@static if Base.VERSION < v"1.11"
+    let
+       function mulfix!(dest, A, B)
+            axdest, axA, axB = axes(dest), axes(A), axes(B)
+            axA[2] == axB[1] || throw(DimensionMismatch("Tried to multiply arrays with axes $axA and $axB"))
+            axdest == (axA[1], axB[2]) || throw(DimensionMismatch("Tried to multiply arrays with axes $axA and $axB and assign to array with axes $axdest"))
+            fill!(dest, zero(eltype(dest)))
+            # This is not maximally efficient, but it produces the right answer on older Julia versions.
+            for ij in CartesianIndices(dest)
+                for k in axA[2]
+                    dest[ij] += A[ij[1], k] * B[k, ij[2]]
+                end
+            end
+            return dest
+       end
+        function LinearAlgebra.mul!(
+            dest::AbstractMatrix{<:StaticMatrix},
+            A::LinearAlgebra.AbstractTriangular{<:StaticVector},
+            B::Adjoint{Adjoint{Float64, V}, <:AbstractMatrix{V}}
+        ) where V<:StaticVector
+            mulfix!(dest, A, B)
+        end
+        function LinearAlgebra.mul!(
+            dest::AbstractMatrix{<:StaticMatrix},
+            A::AbstractMatrix{<:StaticVector},
+            B::Adjoint{Adjoint{Float64, V}, <:AbstractMatrix{V}}
+        ) where V<:StaticVector
+            mulfix!(dest, A, B)
+        end
+    end
+end
+
 """
     multiplied_dimension(A, B)
 

test/matrix_multiply_add.jl

@@ -1,6 +1,7 @@
 using StaticArrays
 using LinearAlgebra
 using BenchmarkTools
+using Random
 using Test
 
 macro test_noalloc(ex)
@@ -226,7 +227,7 @@ function test_wrappers_for_size(N, test_block)
         A_block = rand(SMatrix{N,N,SMatrix{2,2,Int,4}})
         B_block = rand(SMatrix{N,N,SMatrix{2,2,Int,4}})
         bv_block = rand(SVector{N,SMatrix{2,2,Int,4}})
-        
+
         # matrix-vector
         for wrapper in mul_add_wrappers
             # LinearAlgebra can't handle these
@@ -252,3 +253,12 @@ end
     test_wrappers_for_size(8, false)
     test_wrappers_for_size(16, false)
 end
+
+@testset "Adjoints and covariances" begin
+    X = [randn(SVector{3,Float64}) for _ in CartesianIndices((1:2, 1:3))]
+    μ = sum(X; dims=2)/size(X, 2)
+    ΔX = X .- μ
+    @test (ΔX*ΔX')[1, 1] ≈ sum(dx * dx' for dx in ΔX[1, :])
+    B = [randn(SVector{3,Float64}) for _ in CartesianIndices((1:2, 1:2))]
+    @test_throws DimensionMismatch ΔX * B'
+end