Adds multi-dimensionality to ak.matmul to match numpy

[drculhane]

Closes #4961.

Adds multi-dim behavior to ak.matmul to match numpy.

As noted in comments in the matmul function, there are multiple cases:

If either of both arguments is/are 1-D, then matmul matches dot, so the args are sent there.

If both arrays are 2-D, the computation is handled as a conventional matrix-matrix multiplication. The shapes must be compatible. (Addendum: I had to fix a circular calling issue here -- dot also allows 2-D matrix-matrix multiplication, and was sending its arguments to matmul. Since that had each calling the other, I pulled the 2-D computation into a separate helper function named _matmul2D, that both dot and matmul call.)

If one array is > 2-D and the other is >= 2-D, the numpy rules for broadcasting shapes are followed, e.g.

• given argument shapes of [m,n,k,p] * [1,n,p,q], the second argument is first broadcast to shape [m,n,p,q]
• then [m,n] sub-results of shape [k,q] are calculated
• giving a final result of shape [m,n,k,q].

Python-side, the matmul function has been significantly rewritten.

There are new chapel functions to handle the multi-dim case (multidimmatmul and matmulShape) in LinalgMsg.chpl.

There are unit tests that check the multi-dim case up to 3-D.

Note about dot:

In the course of sending things to dot, I discovered that it had a case I hadn't handled when I wrote ak.dot: where the left argument is 1-D and the right isn't. numpy.dot does this, although the online documentation doesn't explicitly describe it. So I added that.

[ajpotts]

Logic looks good! I had one question about performance.

[1RyanK]

Looks good!