Unify `kernelmatrix`, `kernelkronmat`, `kronecker_kernelmatrix`, and `kernelpdmat`?

It's late and I haven't thought much about it - but maybe we could unify the different kernel matrix functions by using

# maybe the default dispatch  should depend on `k`, `x`, and `y`?
kernelmatrix(k, x, y) = kernelmatrix(Matrix, k, x, y)
kernelmatrix(k, x) = kernelmatrix(Matrix, k, x)

function kernelmatrix(::Type{T}, k, x, y)::T where {T<:Matrix}
    ...
end
function kernelmatrix(::Type{T}, k, x)::T where {T<:Matrix}
    ...
end

# instead of kronecker_kernelmatrix
function kernelmatrix(::Type{T}, k::MOKernel, x::IsotopicMOInputsUnion, y::IsotopicMOInputsUnion)::T where {T<:KroneckerProduct}
    ...
end
...

# instead of kernelpdmat
function kernelmatrix(::Type{T}, k, x)::T where {T<:PDMat}
    ...
end

I'm sure some parts of the sketch could be improved but the main idea would be to unify the different functions by using an additional argument with the desired (possibly abstract) container type for dispatching.

I assume this would also make it a bit more convenient to eg extend kernelmatrix to CuArray or static arrays.

One could even include kernelkronmat (I assume) if we define a special input type that we can dispatch on:

function kernelmatrix(::Type{T}, k, xgrid::Gridtype)::T where {T<:KroneckerProduct}
...
end

Start by reading the existing implementations and call sites for kernelmatrix, kernelkronmat, kronecker_kernelmatrix, and kernelpdmat. Compare their dispatch and container-type behavior, then define the scope and compatibility requirements for a unified API; done means the maintainers agree on the design and the affected functionality remains covered.