Faster dot product for symmetric matrices

[lkapelevich]

I want to get an inner product between two symmetric matrices, and dot seems to dispatch to a pretty slow method.
E.g.

using BenchmarkTools, LinearAlgebra
n = 5000;
A = Symmetric(randn(n, n));
B = Symmetric(randn(n, n));
function symdot(A::Symmetric{T, Matrix{T}}, B::Symmetric{T, Matrix{T}}) where {T}
    ret = zero(T)
    m = size(A, 2)
    @inbounds for j in 1:m
        for i in 1:(j - 1)
            ret += A[i, j] * B[i, j] * 2
        end
        ret += A[j, j] * B[j, j]
    end
    return ret
end

julia> @benchmark dot($A, $B)
benchmark symdot($A, $B)
BenchmarkTools.Trial:
  memory estimate:  0 bytes
  allocs estimate:  0
  --------------
  minimum time:     371.083 ms (0.00% GC)
  median time:      371.527 ms (0.00% GC)
  mean time:        371.907 ms (0.00% GC)
  maximum time:     377.313 ms (0.00% GC)
  --------------
  samples:          14
  evals/sample:     1

julia> @benchmark symdot($A, $B)
BenchmarkTools.Trial:
  memory estimate:  0 bytes
  allocs estimate:  0
  --------------
  minimum time:     39.275 ms (0.00% GC)
  median time:      39.359 ms (0.00% GC)
  mean time:        39.378 ms (0.00% GC)
  maximum time:     39.814 ms (0.00% GC)
  --------------
  samples:          127
  evals/sample:     1

Any chance of getting a faster dot product in?

This is very specific, but related to #25565 and #16573.

[andreasnoack]

It should be fine to add a specialized method for symmetric matrices even though it's a quite special case. Please open a PR.

[tkoolen]

Might be significantly faster to linearly index into the parents if they are IndexLinear() (and possibly only if the uplos match). Still, that's a common case.

[goggle]

The suggested solution is only much faster if both matrices A and B are uplo='U'.
Here are some benchmarks:

using LinearAlgebra
using BenchmarkTools
n = 5000;
X = randn(n, n);
Y = randn(n, n);

A = Symmetric(X);
B = Symmetric(Y);

julia> @btime symdot(A, B)
  25.900 ms (1 allocation: 16 bytes)
-1182.2972818525195
julia> @btime dot(A, B)
  306.933 ms (1 allocation: 16 bytes)
-1182.2972818517771

B = Symmetric(Y, :L);

julia> @btime symdot(A, B)
  121.021 ms (1 allocation: 16 bytes)
194.59762244317895
julia> @btime dot(A, B)
  276.302 ms (1 allocation: 16 bytes)
194.59762244405476

A = Symmetric(X, :L);

julia> @btime symdot(A, B)
  266.479 ms (1 allocation: 16 bytes)
-6375.261986764326
julia> @btime dot(A, B)
  310.592 ms (1 allocation: 16 bytes)
-6375.261986764004

So I guess this needs some work to make the specialized method significantly faster in the other cases as well.