/* 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. */
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. */
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,
/* 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
{
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
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;
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;
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 ();
*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;
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);
--- /dev/null
+! { 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" } }