Extend the intrinsic `matmul` to handle more than just two matrices.

[loiseaujc]

Motivation

The Fortran intrinisc matmul can handle only two matrices, i.e. result = matmul(A, B). While working on problems involving multiple low-rank matrices, I've been faced recently with expressions like result = A * B * C, or result = A * B * C * D, where A, B, C and D tend to be tall or wide matrices. Depending on how you evaluate these matrix-matrix multiplications, you can either have a pretty fast code or a terribly slow one.

Except for myself, would it be of interest to anyone to extend the intrinsic matmul to handle expressions involving multiple matrices (say between 2 and 5, I don't think there are situations where you need much more than that). The corresponding interfaces would simply be

result = matmul(A, B) ! Intrinsic matmul
result = matmul(A, B, C)
result = matmul(A, B, C, D)
result = matmul(A, B, C, D, E)

Given the dimensions of each array, the optimal order in which to perform the operations can be determined by solving the corresponding optimization problem using dynamic programming. Some details can be found here.

Prior Art

• numpy exposes the numpy.linalg.multi_dot(arrays, *, out=None) function which chains calls to numpy.dot using the optimal parenthesization of the matrices.

• I'm pretty sure julia has a similar feature but a quick search didn't bring actual references to it

As far as I know, Fortran does not allow for arbitrary long lists of dummy arguments (except for intrisic max and min and a couple of others I suppose). In practice, I don't think the extended matmul interface would need to handle more than say 5-ish arrays though.

Additional Information

Pinging @perazz, @jvdp1 and @jalvesz.