[Patch, fortran] PR35740 - a = conjg(transpose(a)) still gives wrong results, see bug 31994

["Paul Richard Thomas" <paul dot richard dot thomas at gmail dot com>]

:ADDPATCH fortran:

trans-array.c contains an in-line TRANSPOSE routine that interchanges
the dimensions of the descriptor but does not copy the data.  There is
also a library function to implement transpose but this involves a
copy.

 a = transpose(conjg(a)) works correctly but  a = conjg(transpose(a))
does not.  This happens because the lines in trans-intrinsic. that
call the in-line function

      if (se->ss && se->ss->useflags)
	{
	  gfc_conv_tmp_array_ref (se);
	  gfc_advance_se_ss_chain (se);
	}
      else
	gfc_conv_array_transpose (se, expr->value.function.actual->expr);

do not work because the call to gfc_conv_tmp_array_ref (se) is never
made, since CONJG is elemental.  In addition, a temporary should be
made to avoid the dependency between the lhs and the rhs, generated by
the transpose.  Therefore, there is no loss in using the library call,
for this case.  Non-variable expressions are also getting messed up
because the offset was not calculated.  Again, this is fixed by the
library call.

The fix is to make sure that actual argument expressions of elemental
procedures never get marked as non-copying intrinsics.  This is
accomplished in resolve.c.

As a consequence, matmul (transpose (conjg (a)), b)) is going to be
more efficient that matmul (conjg (transpose (a)), b) because the
first can use the inline trick, whereas the second uses the library.
So watch your tildes and stars:)

Bootstrapped and regtested on x86_ia64/FC8 - OK for trunk and 4.3?

Paul

2008-03-30  Paul Thomas  <pault@gcc.gnu.org>

	PR fortran/35740
	* resolve.c (resolve_function, resolve_call): If the procedure
	is elemental do not look for noncopying intrinsics.

2008-03-30  Paul Thomas  <pault@gcc.gnu.org>

	PR fortran/35740
	* gfortran.dg/transpose_conjg_1.f90: New test.

Index: gcc/fortran/resolve.c
===================================================================
*** gcc/fortran/resolve.c	(revision 133709)
--- gcc/fortran/resolve.c	(working copy)
*************** resolve_function (gfc_expr *expr)
*** 2374,2380 ****
        gfc_expr_set_symbols_referenced (expr->ts.cl->length);
      }
  
!   if (t == SUCCESS)
      find_noncopying_intrinsics (expr->value.function.esym,
  				expr->value.function.actual);
  
--- 2374,2385 ----
        gfc_expr_set_symbols_referenced (expr->ts.cl->length);
      }
  
!   if (t == SUCCESS
! 	&& !((expr->value.function.esym
! 		&& expr->value.function.esym->attr.elemental)
! 			||
! 	     (expr->value.function.isym
! 		&& expr->value.function.isym->elemental)))
      find_noncopying_intrinsics (expr->value.function.esym,
  				expr->value.function.actual);
  
*************** resolve_call (gfc_code *c)
*** 2845,2851 ****
    if (resolve_elemental_actual (NULL, c) == FAILURE)
      return FAILURE;
  
!   if (t == SUCCESS)
      find_noncopying_intrinsics (c->resolved_sym, c->ext.actual);
    return t;
  }
--- 2850,2856 ----
    if (resolve_elemental_actual (NULL, c) == FAILURE)
      return FAILURE;
  
!   if (t == SUCCESS && !(c->resolved_sym && c->resolved_sym->attr.elemental))
      find_noncopying_intrinsics (c->resolved_sym, c->ext.actual);
    return t;
  }
Index: gcc/testsuite/gfortran.dg/transpose_conjg_1.f90
===================================================================
*** gcc/testsuite/gfortran.dg/transpose_conjg_1.f90	(revision 0)
--- gcc/testsuite/gfortran.dg/transpose_conjg_1.f90	(revision 0)
***************
*** 0 ****
--- 1,37 ----
+ !  { dg-do run }
+ !  Tests the fix for PR35740, where the trick of interchanging the descriptor
+ !  dimensions to implement TRANSPOSE did not work if it is an argument of
+ !  an elemental function - eg. CONJG.  The fix forces a library call for such
+ !  cases.  During the diagnosis of the PR, it was found that the scalarizer was
+ !  completely thrown if the argument of TRANSPOSE was a non-variable
+ !  expression; eg a + c below.  This is also fixed by the library call.
+ !
+ !  Contributed by Dominik Muth <dominik.muth@gmx.de>
+ !
+ program main
+   implicit none
+   complex, dimension(2,2) :: a,b,c,d
+   a(1,1) = (1.,1.)
+   a(2,1) = (2.,2.)
+   a(1,2) = (3.,3.)
+   a(2,2) = (4.,4.)
+ !
+   b = a
+   b = conjg(transpose(b))
+   d = a
+   d = transpose(conjg(d))
+   if (any (b /= d)) call abort ()
+ !
+   d = matmul (b,  a )
+   if (any (d /= matmul (transpose(conjg(a)), a))) call abort ()
+   if (any (d /= matmul (conjg(transpose(a)), a))) call abort ()
+ !
+   c = (0.0,1.0)
+   b = conjg(transpose(a + c))
+   d = transpose(conjg(a + c))
+   if (any (b /= d)) call abort ()
+ !
+   d = matmul (b,  a + c)
+   if (any (d /= matmul (transpose(conjg(a + c)), a + c))) call abort ()
+   if (any (d /= matmul (conjg(transpose(a + c)), a + c))) call abort ()
+  END program main