re PR fortran/66094 (Handle transpose(A) in inline matmul)

2016-01-24  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-24  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.

From-SVN: r232774

7 files changed, 227 insertions, 10 deletions

diff --git a/gcc/fortran/ChangeLog b/gcc/fortran/ChangeLog
index 9bd4ecf..275acaf 100644
--- a/gcc/fortran/ChangeLog
+++ b/gcc/fortran/ChangeLog
@@ -1,3 +1,14 @@
+2016-01-24  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-21  Jerry DeLisle  <jvdelisle@gcc.gnu.org>
 
 	PR fortran/65996
diff --git a/gcc/fortran/frontend-passes.c b/gcc/fortran/frontend-passes.c
index 9fad41d..340fd6e 100644
--- a/gcc/fortran/frontend-passes.c
+++ b/gcc/fortran/frontend-passes.c
@@ -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 matrix_case m_case)
   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_expr *b,
       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,16 +2720,17 @@ 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
    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 @@ check_conjg_variable (gfc_expr *e, bool *conjg)
 
 	  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_subtrees,
   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_subtrees,
   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_subtrees,
     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_subtrees,
 	  *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_subtrees,
 
       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);
diff --git a/gcc/testsuite/ChangeLog b/gcc/testsuite/ChangeLog
index 0967d44..60d9412 100644
--- a/gcc/testsuite/ChangeLog
+++ b/gcc/testsuite/ChangeLog
@@ -1,3 +1,11 @@
+2016-01-24  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.
+
 2016-01-23  Tom de Vries  <tom@codesourcery.com>
 
 	PR tree-optimization/69426
diff --git a/gcc/testsuite/gfortran.dg/inline_matmul_13.f90 b/gcc/testsuite/gfortran.dg/inline_matmul_13.f90
new file mode 100644
index 0000000..e887ee9
--- /dev/null
+++ b/gcc/testsuite/gfortran.dg/inline_matmul_13.f90
@@ -0,0 +1,47 @@
+! { dg-do run }
+! { dg-options "-ffrontend-optimize -fdump-tree-original -Wrealloc-lhs" }
+! PR 66094: Check functionality for MATMUL(A, TRANSPSE(B))
+module x
+contains
+  subroutine mm1(a,b,c)
+    real, dimension(:,:), intent(in) :: a, b
+    real, dimension(:,:), intent(out) :: c
+    c = -42.
+    c = matmul(a, transpose(b))
+  end subroutine mm1
+end module x
+
+program main
+  use x
+  implicit none
+  integer, parameter :: n = 3, m=4, cnt=2
+  real, dimension(n,cnt) :: a
+  real, dimension(m,cnt) :: b
+  real, dimension(n,m) :: c, cres
+  real, dimension(:,:), allocatable :: calloc
+
+  data a / 2., -3., 5., -7., 11., -13./
+  data b /17., -23., 29., -31., 37., -39., 41., -47./
+  data cres / -225., 356., -396., 227., -360., 392., &
+       -229., 364., -388., 267., -424., 456./ 
+
+  c = matmul(a,transpose(b))
+  if (sum(c-cres)>1e-4) call abort
+  call mm1 (a, b, c)
+  if (sum(c-cres)>1e-4) call abort
+
+  ! Unallocated
+  calloc = matmul(a,transpose(b)) ! { dg-warning "Code for reallocating the allocatable array" }
+  if (any(shape(c) /= shape(calloc))) call abort
+  if (sum(calloc-cres)>1e-4) call abort
+  deallocate(calloc)
+
+  ! Allocated to wrong shape
+  allocate (calloc(10,10))
+  calloc = matmul(a,transpose(b)) ! { dg-warning "Code for reallocating the allocatable array" }
+  if (any(shape(c) /= shape(calloc))) call abort
+  if (sum(calloc-cres)>1e-4) call abort
+  deallocate(calloc)
+
+end program main
+! { dg-final { scan-tree-dump-times "_gfortran_matmul" 0 "original" } }
diff --git a/gcc/testsuite/gfortran.dg/matmul_bounds_10.f90 b/gcc/testsuite/gfortran.dg/matmul_bounds_10.f90
new file mode 100644
index 0000000..6608b49
--- /dev/null
+++ b/gcc/testsuite/gfortran.dg/matmul_bounds_10.f90
@@ -0,0 +1,16 @@
+! { 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" }
diff --git a/gcc/testsuite/gfortran.dg/matmul_bounds_8.f90 b/gcc/testsuite/gfortran.dg/matmul_bounds_8.f90
new file mode 100644
index 0000000..2764cf3
--- /dev/null
+++ b/gcc/testsuite/gfortran.dg/matmul_bounds_8.f90
@@ -0,0 +1,16 @@
+! { 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" }
diff --git a/gcc/testsuite/gfortran.dg/matmul_bounds_9.f90 b/gcc/testsuite/gfortran.dg/matmul_bounds_9.f90
new file mode 100644
index 0000000..5552e40
--- /dev/null
+++ b/gcc/testsuite/gfortran.dg/matmul_bounds_9.f90
@@ -0,0 +1,23 @@
+! { 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