[patch, fortran] Inline MATMUL(A,TRANSPOSE(B)), PR 66094

[Thomas Koenig <tkoenig at netcologne dot de>]

Hello world,

after the recent discussion about MATMUL(A,TRANSPOSE(B)) I have prepared
a patch which implements this.

The patch is a rather straightforward implementation using the machinery
that is already in place for matmul inlining, so the risk of introducing
a regression should be quite low.

Because of this, and beacuse this is a significant use case, I think
this could still be OK for trunk.  Jerry already indicated in the PR
that this could be the case.

Regression-tested.

So... comments?  Toon, would this help you?  Could yo maybe give this
a spin?

Regards

	Thomas


2016-01-17  Thomas Koenig  <tkoenig@gcc.gnu.org>

        PR fortran/66094
        * frontend-passes.c (enum matrix_case):  Add case A2B2T for
        MATMUL(A,TRANSPoSE(B)) where A and B are rank 2.
        (inline_limit_check):  Also add A2B2T.
        (matmul_lhs_realloc):  Handle A2B2T.
        (check_conjg_variable):  Rename to
        (check_conjg_transpose_variable):  and also count TRANSPOSE.
        (inline_matmul_assign):  Handle A2B2T.

2016-01-17  Thomas Koenig  <tkoenig@gcc.gnu.org>

        PR fortran/66094
        * gfortran.dg/inline_matmul_13.f90:  New test.
        * gfortran.dg/matmul_bounds_8.f90:  New test.
        * gfortran.dg/matmul_bounds_9.f90:  New test.
        * gfortran.dg/matmul_bounds_10.f90:  New test.

Index: frontend-passes.c
===================================================================
--- frontend-passes.c	(Revision 232467)
+++ frontend-passes.c	(Arbeitskopie)
@@ -106,7 +106,7 @@ static int var_num = 1;
 
 /* What sort of matrix we are dealing with when inlining MATMUL.  */
 
-enum matrix_case { none=0, A2B2, A2B1, A1B2 };
+enum matrix_case { none=0, A2B2, A2B1, A1B2, A2B2T };
 
 /* Keep track of the number of expressions we have inserted so far 
    using create_var.  */
@@ -2080,7 +2080,7 @@ inline_limit_check (gfc_expr *a, gfc_expr *b, enum
   gfc_typespec ts;
   gfc_expr *cond;
 
-  gcc_assert (m_case == A2B2);
+  gcc_assert (m_case == A2B2 || m_case == A2B2T);
 
   /* Calculation is done in real to avoid integer overflow.  */
 
@@ -2240,6 +2240,18 @@ matmul_lhs_realloc (gfc_expr *c, gfc_expr *a, gfc_
       cond = build_logical_expr (INTRINSIC_OR, ne1, ne2);
       break;
 
+    case A2B2T:
+      ar->start[0] = get_array_inq_function (GFC_ISYM_SIZE, a, 1);
+      ar->start[1] = get_array_inq_function (GFC_ISYM_SIZE, b, 1);
+
+      ne1 = build_logical_expr (INTRINSIC_NE,
+				get_array_inq_function (GFC_ISYM_SIZE, c, 1),
+				get_array_inq_function (GFC_ISYM_SIZE, a, 1));
+      ne2 = build_logical_expr (INTRINSIC_NE,
+				get_array_inq_function (GFC_ISYM_SIZE, c, 2),
+				get_array_inq_function (GFC_ISYM_SIZE, b, 1));
+      cond = build_logical_expr (INTRINSIC_OR, ne1, ne2);
+
     case A2B1:
       ar->start[0] = get_array_inq_function (GFC_ISYM_SIZE, a, 1);
       cond = build_logical_expr (INTRINSIC_NE,
@@ -2708,7 +2720,7 @@ has_dimen_vector_ref (gfc_expr *e)
 
 /* If handed an expression of the form
 
-   CONJG(A)
+   TRANSPOSE(CONJG(A))
 
    check if A can be handled by matmul and return if there is an uneven number
    of CONJG calls.  Return a pointer to the array when everything is OK, NULL
@@ -2715,9 +2727,10 @@ has_dimen_vector_ref (gfc_expr *e)
    otherwise. The caller has to check for the correct rank.  */
 
 static gfc_expr*
-check_conjg_variable (gfc_expr *e, bool *conjg)
+check_conjg_transpose_variable (gfc_expr *e, bool *conjg, bool *transpose)
 {
   *conjg = false;
+  *transpose = false;
 
   do
     {
@@ -2733,6 +2746,8 @@ static gfc_expr*
 
 	  if (e->value.function.isym->id == GFC_ISYM_CONJG)
 	    *conjg = !*conjg;
+	  else if (e->value.function.isym->id == GFC_ISYM_TRANSPOSE)
+	    *transpose = !*transpose;
 	  else return NULL;
 	}
       else
@@ -2789,7 +2804,7 @@ inline_matmul_assign (gfc_code **c, int *walk_subt
   int i;
   gfc_code *if_limit = NULL;
   gfc_code **next_code_point;
-  bool conjg_a, conjg_b;
+  bool conjg_a, conjg_b, transpose_a, transpose_b;
 
   if (co->op != EXEC_ASSIGN)
     return 0;
@@ -2809,12 +2824,12 @@ inline_matmul_assign (gfc_code **c, int *walk_subt
   changed_statement = NULL;
 
   a = expr2->value.function.actual;
-  matrix_a = check_conjg_variable (a->expr, &conjg_a);
-  if (matrix_a == NULL)
+  matrix_a = check_conjg_transpose_variable (a->expr, &conjg_a, &transpose_a);
+  if (transpose_a || matrix_a == NULL)
     return 0;
 
   b = a->next;
-  matrix_b = check_conjg_variable (b->expr, &conjg_b);
+  matrix_b = check_conjg_transpose_variable (b->expr, &conjg_b, &transpose_b);
   if (matrix_b == NULL)
     return 0;
 
@@ -2828,10 +2843,28 @@ inline_matmul_assign (gfc_code **c, int *walk_subt
     return 0;
 
   if (matrix_a->rank == 2)
-    m_case = matrix_b->rank == 1 ? A2B1 : A2B2;
+    {
+      if (matrix_b->rank == 1)
+	m_case = A2B1;
+      else
+	{
+	  if (transpose_b)
+	    m_case = A2B2T;
+	  else
+	    m_case = A2B2;
+	}
+    }
   else
-    m_case = A1B2;
+    {
+      /* Vector * Transpose(B) not handled yet.  */
+      if (transpose_b)
+	m_case = none;
+      else
+	m_case = A1B2;
+    }
 
+  if (m_case == none)
+    return 0;
 
   ns = insert_block ();
 
@@ -3002,6 +3035,36 @@ inline_matmul_assign (gfc_code **c, int *walk_subt
 	  *next_code_point = test;
 	  next_code_point = &test->next;
 	}
+
+      if (m_case == A2B2T)
+	{
+	  c1 = get_array_inq_function (GFC_ISYM_SIZE, expr1, 1);
+	  a1 = get_array_inq_function (GFC_ISYM_SIZE, matrix_a, 1);
+	  test = runtime_error_ne (c1, a1, "Incorrect extent in return array in "
+				   "MATMUL intrinsic for dimension 1: "
+				   "is %ld, should be %ld");
+
+	  *next_code_point = test;
+	  next_code_point = &test->next;
+
+	  c2 = get_array_inq_function (GFC_ISYM_SIZE, expr1, 2);
+	  b1 = get_array_inq_function (GFC_ISYM_SIZE, matrix_b, 1);
+	  test = runtime_error_ne (c2, b1, "Incorrect extent in return array in "
+				   "MATMUL intrinsic for dimension 2: "
+				   "is %ld, should be %ld");
+	  *next_code_point = test;
+	  next_code_point = &test->next;
+
+	  a2 = get_array_inq_function (GFC_ISYM_SIZE, matrix_a, 2);
+	  b2 = get_array_inq_function (GFC_ISYM_SIZE, matrix_b, 2);
+
+	  test = runtime_error_ne (b2, a2, "Incorrect extent in argument B in "
+				   "MATMUL intrnisic for dimension 2: "
+				   "is %ld, should be %ld");
+	  *next_code_point = test;
+	  next_code_point = &test->next;
+
+	}
     }
 
   *next_code_point = assign_zero;
@@ -3050,6 +3113,39 @@ inline_matmul_assign (gfc_code **c, int *walk_subt
 
       break;
 
+    case A2B2T:
+      inline_limit_check (matrix_a, matrix_b, m_case);
+
+      u1 = get_size_m1 (matrix_b, 1);
+      u2 = get_size_m1 (matrix_a, 2);
+      u3 = get_size_m1 (matrix_a, 1);
+
+      do_1 = create_do_loop (gfc_copy_expr (zero), u1, NULL, &co->loc, ns);
+      do_2 = create_do_loop (gfc_copy_expr (zero), u2, NULL, &co->loc, ns);
+      do_3 = create_do_loop (gfc_copy_expr (zero), u3, NULL, &co->loc, ns);
+
+      do_1->block->next = do_2;
+      do_2->block->next = do_3;
+      do_3->block->next = assign_matmul;
+
+      var_1 = do_1->ext.iterator->var;
+      var_2 = do_2->ext.iterator->var;
+      var_3 = do_3->ext.iterator->var;
+
+      list[0] = var_3;
+      list[1] = var_1;
+      cscalar = scalarized_expr (co->expr1, list, 2);
+
+      list[0] = var_3;
+      list[1] = var_2;
+      ascalar = scalarized_expr (matrix_a, list, 2);
+
+      list[0] = var_1;
+      list[1] = var_2;
+      bscalar = scalarized_expr (matrix_b, list, 2);
+
+      break;
+
     case A2B1:
       u1 = get_size_m1 (matrix_b, 1);
       u2 = get_size_m1 (matrix_a, 1);

! { dg-do run }
! { dg-options "-fno-backtrace -fbounds-check -fno-realloc-lhs" }
! { dg-shouldfail "Fortran runtime error: Incorrect extent in return array in MATMUL intrinsic for dimension 2: is 2, should be 3" }
program main
  real, dimension(3,2) :: a
  real, dimension(3,2) :: b
  real, dimension(:,:), allocatable :: ret
  allocate (ret(3,3))
  a = 1.0
  b = 2.3
  ret = matmul(a,transpose(b))  ! This is OK
  deallocate(ret)
  allocate(ret(3,2))
  ret = matmul(a,transpose(b))  ! This should throw an error.
end program main
! { dg-output "Fortran runtime error: Incorrect extent in return array in MATMUL intrinsic for dimension 2: is 2, should be 3" }

! { dg-do run }
! { dg-options "-fbounds-check -ffrontend-optimize" }
! { dg-shouldfail "Fortran runtime error: Incorrect extent in argument B in MATMUL intrnisic for dimension 2: is 1, should be 2" }
module x
  implicit none
contains
  subroutine mmul(c, a, b)
    real, dimension(:,:), intent(in) :: a,b
    real, dimension(:,:), intent(out) :: c
    c = matmul(a,transpose(b))
  end subroutine mmul
end module x

program main
  use x
  integer, parameter :: n = 3, m=4, cnt=2
  real, dimension(n,cnt) :: a
  real, dimension(m,cnt-1) :: b
  real, dimension(n,m) :: c
  a = 1.0
  b = 2.3
  call mmul(c,a,b)
end program main

! { dg-do run }
! { dg-options "-fno-backtrace -fbounds-check -fno-realloc-lhs" }
! { dg-shouldfail "Fortran runtime error: Incorrect extent in return array in MATMUL intrinsic for dimension 1: is 4, should be 3" }
program main
  real, dimension(3,2) :: a
  real, dimension(3,2) :: b
  real, dimension(:,:), allocatable :: ret
  allocate (ret(3,3))
  a = 1.0
  b = 2.3
  ret = matmul(a,transpose(b))  ! This is OK
  deallocate(ret)
  allocate(ret(4,3))
  ret = matmul(a,transpose(b))  ! This should throw an error.
end program main
! { dg-output "Fortran runtime error: Incorrect extent in return array in MATMUL intrinsic for dimension 1: is 4, should be 3" }