LinearAlgebra: make matprod_dest public (#55537)

Currently, in a matrix multiplication `A * B`, we use `B` to construct
the destination. However, this may not produce the optimal destination
type, and is essentially single-dispatch. Letting packages specialize
`matprod_dest` would help us obtain the optimal type by dispatching on
both the arguments. This may significantly improve performance in the
matrix multiplication. As an example:
```julia
julia> using LinearAlgebra, FillArrays, SparseArrays

julia> F = Fill(3, 10, 10);

julia> s = sprand(10, 10, 0.1);

julia> @Btime $F * $s;
  15.225 μs (10 allocations: 4.14 KiB)

julia> typeof(F * s)
SparseMatrixCSC{Float64, Int64}

julia> nnz(F * s)
80

julia> VERSION
v"1.12.0-DEV.1074"
```
In this case, the destination is a sparse matrix with 80% of its
elements filled and being set one-by-one, which is terrible for
performance. Instead, if we specialize `matprod_dest` to return a dense
destination, we may obtain
```julia
julia> LinearAlgebra.matprod_dest(F::FillArrays.AbstractFill, S::SparseMatrixCSC, ::Type{T}) where {T} = Matrix{T}(undef, size(F,1), size(S,2))

julia> @Btime $F * $s;
  754.632 ns (2 allocations: 944 bytes)

julia> typeof(F * s)
Matrix{Float64}
```
Potentially, this may be improved further by specializing `mul!`, but
this is a 20x improvement just by choosing the right destination.

Since this is being made public, we may want to bikeshed on an
appropriate name for the function.

Author: jishnub

stdlib/LinearAlgebra/src/LinearAlgebra.jl

@@ -174,7 +174,8 @@ public AbstractTriangular,
         isbanded,
         peakflops,
         symmetric,
-        symmetric_type
+        symmetric_type,
+        matprod_dest
 
 const BlasFloat = Union{Float64,Float32,ComplexF64,ComplexF32}
 const BlasReal = Union{Float64,Float32}

stdlib/LinearAlgebra/src/matmul.jl

@@ -123,7 +123,15 @@ function (*)(A::AbstractMatrix, B::AbstractMatrix)
     mul!(matprod_dest(A, B, TS), A, B)
 end
 
-matprod_dest(A, B, TS) = similar(B, TS, (size(A, 1), size(B, 2)))
+"""
+    matprod_dest(A, B, T)
+
+Return an appropriate `AbstractArray` with element type `T` that may be used to store the result of `A * B`.
+
+!!! compat
+    This function requires at least Julia 1.11
+"""
+matprod_dest(A, B, T) = similar(B, T, (size(A, 1), size(B, 2)))
 
 # optimization for dispatching to BLAS, e.g. *(::Matrix{Float32}, ::Matrix{Float64})
 # but avoiding the case *(::Matrix{<:BlasComplex}, ::Matrix{<:BlasReal})