@@ -424,6 +424,15 @@ function operate(
424424 # concrete. Bad things can happen if S or T is abstract and we pick the
425425 # wrong type for C.
426426 if ! (isconcretetype (S) && isconcretetype (T))
427+ if S <: AbstractMutable || T <: AbstractMutable
428+ # We can't use A * B here, because this will StackOverflow via:
429+ # *(A, B) -> mul(A, B) -> operate(*, A, B) -> *(A, B)
430+ # See MutableArithmetics.jl#336
431+ if axes (A, 2 ) != axes (B, 1 )
432+ throw (DimensionMismatch ())
433+ end
434+ return [sum (A[i, j] * B[j] for j in axes (B, 1 )) for i in axes (A, 1 )]
435+ end
427436 return A * B
428437 end
429438 C = undef_array (promote_array_mul (typeof (A), typeof (B)), axes (A, 1 ))
@@ -439,6 +448,9 @@ function operate(
439448 # concrete. Bad things can happen if S or T is abstract and we pick the
440449 # wrong type for C.
441450 if ! (isconcretetype (S) && isconcretetype (T))
451+ # It's safe to use A * B here, because there is no fallback for
452+ # *(A, B) -> mul(A, B) -> operate(*, A, B)
453+ # See MutableArithmetics.jl#336
442454 return A * B
443455 end
444456 C = undef_array (
0 commit comments