Also, in a separate code path, for evalpoly(zeros(2,2), (1,2,3))

Seems like a Base issue where it is trying to add a number to a matrix? Perhaps evalpoly may add p[i]*J instead of p[i], where J is an appropriate identity element. I'm unsure which one to use in general: perhaps oneunit?

It's just a matter of how the sum is initialized — it should be initialized by one(x) * p[end] rather than by p[end].

(Note that it should be one, not oneunit, because the former corresponds to x^0.)

However, I experimented with this and it seems like we need a completely separate method when x is a matrix and the coefficients are scalars, because of #1162 — even if you fix the sum initialization, it gives completely wrong results because of the weird behavior of muladd.

(Also, a specialized method for matrix x could be a good idea anyway because it could work more in-place.)

So, for now, the fact that it throws an error is a good thing — better than giving the wrong answer.

Here is a generic implementation that we can use as a fallback, which should handle most matrix types AFAICT:

# non-inplace fallback
function _evalpoly(X::AbstractMatrix, p)
    Base.require_one_based_indexing(p)
    p0 = isempty(p) ? Base.reduce_empty_iter(+, p) : p[end]
    Xone = one(X)
    S = Base.promote_op(*, typeof(Xone), typeof(Xone))(Xone) * p0
    for i = length(p)-1:-1:1
        S = X * S + @inbounds(p[i] isa AbstractMatrix ? p[i] : p[i] * I)
    end
    return S
end

evalpoly(X::AbstractMatrix, p::Tuple) = _evalpoly(X, p)
evalpoly(X::AbstractMatrix, p::AbstractVector) = _evalpoly(X, p)

It would be nice to have an optimized implementation when p is an array or tuple of scalars that works in-place, since such matrix polynomial functions are pretty common, but it's a little tricky to decide what types of X it should handle. Just X::StridedMatrix{<:Number}, maybe? (Otherwise you can easily run into trouble for things like SMatrix or other exotic types.)

A generic in place version could work using similar. However, in the case of SMatrix it may sometimes be more performant to use an out of place version that returns an SMatrix than an inplace version that returns a heap allocated MMatrix. I think that your proposed generic fallback is worth merging.

A generic in place version could work using similar.

You can't just use similar(X), because for some matrix types (e.g. SymTridiagonal) the type of X * X is different from that of X or similar(X) — that's why I used Base.promote_op in the code above. I guess you could call the promote_op conversion on similar(X), though this might possibly make an extra copy. It also may be tricky to get in-place code that works well for sparse matrix types that don't support efficient random mutation.

I also agree that for SMatrix you would probably want the out-of-place version.

I was thinking of S = Base.promote_op(similar ∘ *, typeof(Xone), typeof(Xone))(Xone) * p0

You also run into trouble for matrices of matrices, but I guess those aren't well-supported in linear algebra anyway right now. e.g. one(X) throws a MethodError for X = [rand(2,2) for i=1:2, j=1:2], even though in principle this should be fine.