Move reinterpret-based optimization for complex matrix * real vec/mat to lower level.

[N5N3]

Inspired by #44011.
This has 2 advantages:

1. reinterpret-based optimiztion is faster only for BLAS-compatible cases. This PR ensures that it does not cause regression.
2. reinterpret(T, A::StridedArray) might not be a StridedArray. But we can still call BLAS api as we have strides defined.

The current test should cover all cases. (see #29246)
Edit: Test size is enlarged to 5x5 to enable BLAS fallback.

[dkarrasch]

LGTM. The linux32 test failure is concerning, as it seems to hang in the early linalg tests, one of which is the matmul tests touched here. I'll restart CI.

[Moelf]

https://discourse.julialang.org/t/in-place-multiplication-is-too-much-slower-for-complexf64/75887/

while you're at it do you want to fix this too?

[dkarrasch]

https://discourse.julialang.org/t/in-place-multiplication-is-too-much-slower-for-complexf64/75887/

while you're at it do you want to fix this too?

I guess that's taken care of by #44011?

[N5N3]

In general we don't want extra allocation when use mul!(C,A,B). Thus the convert-based optimization is applied only to A*B.
Unfortunately, BLAS's gemm only support inputs with dense column. (Thus #44011 only handles "mat * vec" case)
MKL has a extended api cblas_?gemv_batch_strided. Looks like we can use it to do gemm with non-dense column.
But I never used it. Edit: have some try

using MKL_jll
const Lib = MKL_jll.libmkl_rt
import LinearAlgebra.BlasInt
import LinearAlgebra.BLAS: @blasfunc
for (gemv_batch_strided, elty) in ((:(:dgemv_batch_strided),:Float64))
    @eval function gemv_batch_strided!(transA::AbstractChar,
                alpha::Union{($elty), Bool},
                A::AbstractMatrix{$elty}, B::AbstractVecOrMat{$elty},
                beta::Union{($elty), Bool},
                C::AbstractVecOrMat{$elty})
        m, n = size(A)
        batch = size(C, 2)
        ccall(($gemv_batch_strided, Lib), Cvoid,
            (Ref{UInt8}, Ref{BlasInt}, Ref{BlasInt}, Ref{$elty}, Ptr{$elty}, Ref{BlasInt}, Ref{BlasInt},
            Ptr{$elty}, Ref{BlasInt}, Ref{BlasInt}, Ref{$elty}, Ptr{$elty}, Ref{BlasInt}, Ref{BlasInt},Ref{BlasInt}),
            transA, m, n, alpha, A, max(1,stride(A,2)), 0, 
            B, stride(B, 1), stride(B, 2), beta, C, stride(C, 1), stride(C, 2), batch)
        C
    end
end

The performance is not good though:

julia> a = randn(1000, 1000);
julia> b = view(randn(2000, 1000), 1:2:2000, :);
julia> c = zeros(1000,1000);

julia> gemv_batch_strided!('N', true, a, b, false, c) ≈ a * b
true

julia> @btime gemv_batch_strided!('N', true, $a, $b, false, $c);
  39.001 ms (0 allocations: 0 bytes) # MKL use 6 threads

julia> @btime mul!($c,$a,$b)
 687.478 ms (3 allocations: 30.75 KiB) # general matmul only use 1 thread