I noticed that for large matrices, MATMUL is crashing with a stack overflow. I can fix this with /heap-arrays0. The program does not crash when calling dgemm from MKL. I ran some tests, and the results from MATMUL and dgemm are identical. However MATMUL needs a large stack, and dgemm doesn't. Is this the correct behavior, or is there a bug somewhere?
Most of the time, the MATMUL intrinsic is implemented by in-line code. That is, the actual instructions to multiply each element of the matrices is generated by the compiler (including the loops to go from element to element, etc). That's different than a call to dgemm, where it is a simple routine call.
That said, because using /heap-arrays resolves the stack overflow, I'm going to guess that the compiler cannot detect that there won't be overlap between the result variable and the two operands, and for safety reasons it creates a temporary array. By default, this is on the stack.
To answer your question, I would characterize this as "expected behavior", and not likely a bug.
I ran a test, and the performance of dgemm is slightly better than matmul. Based on your reply, this might be because there is an extra step where values are copied from the temporary memory to the result variable.