Support different matrix sizes in cauchy_dot_product

Right now product_by_order is secretly only valid for 2x2 matrices, and will silently produce nonsense with other array sizes. We should either check that the matrices have proper sizes or better support arbitrary sizes.