Repository: stfc/PSyclone

Issue:
matmul2code transformation produces incorrect code
Embarrassingly, my extension of the matmul2code trans does not produce the correct code for matrix-matrix multiplication.
Currently it gives:

          do j = 1, 3, 1
            do i = 1, 3, 1
              jac(i,j) = 0.0_r_def
              do ii = 1, 3, 1
                jac(i,j) = jac(i,j) + jac_inv(ii,i) * jac_inv_transpose(j,ii)
              enddo
            enddo
          enddo

when it should be:

      jac(i,j) = jac(i,j) + jac_inv(i,ii) * jac_inv_transpose(ii,j)

Discussion:

Patch:
diff --git a/changelog b/changelog
index 90ff1c785..ac45fbeb7 100644
--- a/changelog
+++ b/changelog
@@ -68,6 +68,8 @@
 
 	23) PR #1725. Allows parsing an empty Fortran file.
 
+	24) PR #1716 for #1713. Fix bug in matmul transformation. 
+
 release 2.2.0 17th March 2022
 
 	1) PR #1439 for #1074. Adds a GOcean example for the use of the NaN-
diff --git a/src/psyclone/psyir/transformations/intrinsics/matmul2code_trans.py b/src/psyclone/psyir/transformations/intrinsics/matmul2code_trans.py
index 89ef0e59b..2225d1a6b 100644
--- a/src/psyclone/psyir/transformations/intrinsics/matmul2code_trans.py
+++ b/src/psyclone/psyir/transformations/intrinsics/matmul2code_trans.py
@@ -138,11 +138,11 @@ class Matmul2CodeTrans(Operator2CodeTrans):
 
         do i=1,N
             R(i) = 0.0
-                do j=1,M
-                    R(i) = R(i) + A(i,j) * B(j)
+            do j=1,M
+                R(i) = R(i) + A(i,j) * B(j)
 
     For a matrix-matrix multiplication, if the dimensions of ``R``, ``A``,
-    and ``B`` are ``R(P,M)``, ``A(N,M)``, ``B(P,N)``, the MATMUL is replaced
+    and ``B`` are ``R(P,M)``, ``A(P,N)``, ``B(N,M)``, the MATMUL is replaced
     with the following code:
 
     .. code-block:: fortran
@@ -151,7 +151,7 @@ class Matmul2CodeTrans(Operator2CodeTrans):
             do i=1,P
                 R(i,j) = 0.0
                 do ii=1,N
-                    R(i,j) = R(i,j) + A(ii,i) * B(j,ii)
+                    R(i,j) = R(i,j) + A(i,ii) * B(ii,j)
 
     Note that this transformation does *not* support the case where ``A`` is
     a rank-1 array.
@@ -427,22 +427,22 @@ class Matmul2CodeTrans(Operator2CodeTrans):
         result_ref = _create_matrix_ref(result.symbol,
                                         [i_loop_sym, j_loop_sym],
                                         result.children[2:])
-        # Create "matrix2(j,ii)"
+        # Create "matrix2(ii,j)"
         m2_array_reference = _create_matrix_ref(matrix2.symbol,
-                                                [j_loop_sym, ii_loop_sym],
+                                                [ii_loop_sym, j_loop_sym],
                                                 matrix2.children[2:])
-        # Create "matrix1(ii,i)"
+        # Create "matrix1(i,ii)"
         m1_array_reference = _create_matrix_ref(matrix1.symbol,
-                                                [ii_loop_sym, i_loop_sym],
+                                                [i_loop_sym, ii_loop_sym],
                                                 matrix1.children[2:])
-        # Create "matrix1(ii,i) * matrix2(j,ii)"
+        # Create "matrix1(i,ii) * matrix2(ii,j)"
         multiply = BinaryOperation.create(
             BinaryOperation.Operator.MUL, m1_array_reference,
             m2_array_reference)
-        # Create "result(i,j) + matrix1(ii,i) * matrix2(j,ii)"
+        # Create "result(i,j) + matrix1(i,ii) * matrix2(ii,j)"
         rhs = BinaryOperation.create(
             BinaryOperation.Operator.ADD, result_ref, multiply)
-        # Create "result(i,j) = result(i,j) + matrix1(ii,i) * matrix2(j,ii)"
+        # Create "result(i,j) = result(i,j) + matrix1(i,ii) * matrix2(ii,j)"
         assign = Assignment.create(result_ref.copy(), rhs)
         # Create ii loop and add the above code as a child
         # Work out the bounds


Test patch:
diff --git a/src/psyclone/tests/psyir/transformations/intrinsics/matmul2code_trans_test.py b/src/psyclone/tests/psyir/transformations/intrinsics/matmul2code_trans_test.py
index 5f2095ec0..bf5401385 100644
--- a/src/psyclone/tests/psyir/transformations/intrinsics/matmul2code_trans_test.py
+++ b/src/psyclone/tests/psyir/transformations/intrinsics/matmul2code_trans_test.py
@@ -734,7 +734,7 @@ def test_apply_matmat_no_indices(tmpdir, fortran_reader, fortran_writer):
         "    do i = 1, 5, 1\n"
         "      result(i,j) = 0.0\n"
         "      do ii = 1, 3, 1\n"
-        "        result(i,j) = result(i,j) + jac(ii,i) * jac_inv(j,ii)\n"
+        "        result(i,j) = result(i,j) + jac(i,ii) * jac_inv(ii,j)\n"
         "      enddo\n"
         "    enddo\n"
         "  enddo\n"
@@ -773,7 +773,7 @@ def test_apply_matmat_extra_indices(tmpdir, fortran_reader, fortran_writer):
         "      result(i,j,2) = 0.0\n"
         "      do ii = 1, 3, 1\n"
         "        result(i,j,2) = result(i,j,2) + "
-        "jac(ii,i,1) * jac_inv(j,ii,2)\n"
+        "jac(i,ii,1) * jac_inv(ii,j,2)\n"
         "      enddo\n"
         "    enddo\n"
         "  enddo\n" in out)
@@ -811,7 +811,7 @@ def test_apply_matmat_name_clashes(tmpdir, fortran_reader, fortran_writer):
         "      result(i_1,j_1,2) = 0.0\n"
         "      do ii_1 = 1, 3, 1\n"
         "        result(i_1,j_1,2) = result(i_1,j_1,2) + "
-        "jac(ii_1,i_1,1) * jac_inv(j_1,ii_1,2)\n"
+        "jac(i_1,ii_1,1) * jac_inv(ii_1,j_1,2)\n"
         "      enddo\n"
         "    enddo\n"
         "  enddo\n" in out)