public class MM
{
  // Assumes A, B, C have the same domain
  void MM(double array<2> C, double array<2> A, double array<2> B)
  {
    RectDomain<2> d = C.domain;
    Point<2> upperleft = d.min(), bottomright = d.max(), stride = d.stride();

    foreach (c within d)
      {
	double sum = 0;

	// Do the dot product of a row of A with a column of B
	foreach (a within [ [ c[0], upperleft[1] ] : [ c[0], bottomright[1] ] : stride ],
		 b within [ [ upperleft[0], c[1] ] : [ bottomright[0], c[1] ] : stride ])
	  sum += A[a] * B[b];
	C[c] = sum;
      }
  }

  // Using border, this only requires that A, B, C have the same size and strides
  void MMborder(double array<2> C, double array<2> A, double array<2> B)
  {
    RectDomain<2> d = C.domain;
    Point<2> upperleft = d.min();
    RectDomain<2> Arow0 = A.domain.border(1, -1, 0);
    RectDomain<2> Bcolumn0 = B.domain.border(1, -2, 0);

    foreach (c within d)
      {
	double sum = 0;

	// Do the dot product of a row of A with a column of B
	foreach (a within Arow0 + Point<2>.direction(0, c[0] - upperleft[0]),
		 b within Bcolumn0 + Point<2>.direction(1, c[1] - upperleft[1]))
	  sum += A[a] * B[b];
	C[c] = sum;
      }
  }
}