This patch fixes the matmul bug, which was reported this morning.
Seemingly, the code never worked correctly.
Regarding the result patch: Mikael and I came independently to the same
result thus it must be right (or not?); additionally, I tested and found
out that matrix-scalar/scalar-matrix was mishandled.
Build and regtested on x86-64-gnu-linux.
OK for the trunk and the 4.6 and 4.7 branches? (I think we can ignore
4.5, can we?)
Tobias
PS: There is another bug, which prevents the simplification in some
cases (see PR), as it only results in a missed optimization (or reject
valid in constant expressions), I defer this to another patch.
2013-02-14 Tobias Burnus <bur...@net-b.de>
Mikael Morin <mik...@gcc.gnu.org>
PR fortran/56318
* simplify.c (gfc_simplify_matmul): Fix result shape
and matmul result.
2013-02-14 Tobias Burnus <bur...@net-b.de>
PR fortran/56318
* gcc/testsuite/gfortran.dg/matmul_9.f90: New.
diff --git a/gcc/fortran/simplify.c b/gcc/fortran/simplify.c
index f7401e9..a0909a3 100644
--- a/gcc/fortran/simplify.c
+++ b/gcc/fortran/simplify.c
@@ -3850,7 +3850,7 @@ gfc_simplify_matmul (gfc_expr *matrix_a, gfc_expr *matrix_b)
if (matrix_a->rank == 1 && matrix_b->rank == 2)
{
result_rows = 1;
- result_columns = mpz_get_si (matrix_b->shape[0]);
+ result_columns = mpz_get_si (matrix_b->shape[1]);
stride_a = 1;
stride_b = mpz_get_si (matrix_b->shape[0]);
@@ -3860,7 +3860,7 @@ gfc_simplify_matmul (gfc_expr *matrix_a, gfc_expr *matrix_b)
}
else if (matrix_a->rank == 2 && matrix_b->rank == 1)
{
- result_rows = mpz_get_si (matrix_b->shape[0]);
+ result_rows = mpz_get_si (matrix_a->shape[0]);
result_columns = 1;
stride_a = mpz_get_si (matrix_a->shape[0]);
stride_b = 1;
@@ -3873,7 +3873,7 @@ gfc_simplify_matmul (gfc_expr *matrix_a, gfc_expr *matrix_b)
{
result_rows = mpz_get_si (matrix_a->shape[0]);
result_columns = mpz_get_si (matrix_b->shape[1]);
- stride_a = mpz_get_si (matrix_a->shape[1]);
+ stride_a = mpz_get_si (matrix_a->shape[0]);
stride_b = mpz_get_si (matrix_b->shape[0]);
result->rank = 2;
diff --git a/gcc/testsuite/gfortran.dg/matmul_9.f90 b/gcc/testsuite/gfortran.dg/matmul_9.f90
new file mode 100644
index 0000000..bf2a299
--- /dev/null
+++ b/gcc/testsuite/gfortran.dg/matmul_9.f90
@@ -0,0 +1,47 @@
+! { dg-do run }
+! { dg-options "-fdump-tree-original" }
+!
+! PR fortran/56318
+!
+! Contributed by Alberto Luaces
+!
+SUBROUTINE mass_matrix
+ DOUBLE PRECISION,PARAMETER::m1=1.d0
+ DOUBLE PRECISION,DIMENSION(3,2),PARAMETER::A1=reshape([1.d0,0.d0, 0.d0, &
+ 0.d0,1.d0, 0.d0],[3,2])
+ DOUBLE PRECISION,DIMENSION(2,2),PARAMETER::Mel=reshape([1.d0/3.d0, 0.d0, &
+ 0.d0, 1.d0/3.d0],[2,2])
+
+ DOUBLE PRECISION,DIMENSION(3,3)::MM1
+
+ MM1=m1*matmul(A1,matmul(Mel,transpose(A1)))
+ !print '(3f8.3)', MM1
+ if (any (abs (MM1 &
+ - reshape ([1.d0/3.d0, 0.d0, 0.d0, &
+ 0.d0, 1.d0/3.d0, 0.d0, &
+ 0.d0, 0.d0, 0.d0], &
+ [3,3])) > epsilon(1.0d0))) &
+ call abort ()
+END SUBROUTINE mass_matrix
+
+program name
+ implicit none
+ integer, parameter :: A(3,2) = reshape([1,2,3,4,5,6],[3,2])
+ integer, parameter :: B(2,3) = reshape([3,17,23,31,43,71],[2,3])
+ integer, parameter :: C(3) = [-5,-7,-21]
+ integer, parameter :: m1 = 1
+
+! print *, matmul(B,C)
+ if (any (matmul(B,C) /= [-1079, -1793])) call abort()
+! print *, matmul(C,A)
+ if (any (matmul(C,A) /= [-82, -181])) call abort()
+! print '(3i5)', m1*matmul(A,B)
+ if (any (m1*matmul(A,B) /= reshape([71,91,111, 147,201,255, 327,441,555],&
+ [3,3]))) &
+ call abort()
+ call mass_matrix
+end program name
+
+! { dg-final { scan-tree-dump-times "matmul" 0 "original" } }
+! { dg-final { cleanup-tree-dump "original" } }
+
