Update

Author: odow

src/implementations/LinearAlgebra.jl

@@ -424,11 +424,14 @@ function operate(
     # concrete. Bad things can happen if S or T is abstract and we pick the
     # wrong type for C.
     if !(isconcretetype(S) && isconcretetype(T))
-        # We can't use A * B here, because this may StackOverflow via the chain:
-        #   *(A, B) -> mul(A, B) -> operate(*, A, B) -> *(A, B)
-        # See MutableArithmetics.jl#336
-        @assert axes(A, 2) == axes(B, 1)
-        return [sum(A[i, j] * B[j] for j in axes(B, 1)) for i in axes(A, 1)]
+        if S <: AbstractMutable || T <: AbstractMutable
+            # We can't use A * B here, because this will StackOverflow via:
+            #   *(A, B) -> mul(A, B) -> operate(*, A, B) -> *(A, B)
+            # See MutableArithmetics.jl#336
+            @assert axes(A, 2) == axes(B, 1)
+            return [sum(A[i, j] * B[j] for j in axes(B, 1)) for i in axes(A, 1)]
+        end
+        return A * B
     end
     C = undef_array(promote_array_mul(typeof(A), typeof(B)), axes(A, 1))
     return operate_to!(C, *, A, B)