Bug 32962 - b = conjg(transpose(a)) is erroneous if b is an allocatable array
Summary: b = conjg(transpose(a)) is erroneous if b is an allocatable array
Status: RESOLVED FIXED
Alias: None
Product: gcc
Classification: Unclassified
Component: fortran (show other bugs)
Version: 4.2.1
: P3 normal
Target Milestone: 4.3.0
Assignee: Paul Thomas
URL:
Keywords: wrong-code
Depends on:
Blocks: 32834
  Show dependency treegraph
 
Reported: 2007-08-01 19:50 UTC by Wirawan Purwanto
Modified: 2007-08-14 13:46 UTC (History)
4 users (show)

See Also:
Host:
Target:
Build:
Known to work:
Known to fail: 4.3.0 4.2.1 4.1.3
Last reconfirmed: 2007-08-08 16:34:50


Attachments
Failing testcase with higher-dimensional b array (203 bytes, text/plain)
2007-08-01 20:00 UTC, Wirawan Purwanto
Details
Magically successful testcase! (195 bytes, text/plain)
2007-08-01 20:02 UTC, Wirawan Purwanto
Details

Note You need to log in before you can comment on or make changes to this bug.
Description Wirawan Purwanto 2007-08-01 19:50:30 UTC
It turns out that bug #31994 (http://gcc.gnu.org/bugzilla/show_bug.cgi?id=31994) has not been fully resolved!

I found another testcase that would fail gcc 4.2.1, namely if the destination of the expression is an allocatable array. Use the following testcase:

program testcase
  complex(kind=8), allocatable :: a(:,:), b(:,:)
  allocate(a(2,2))
  allocate(b(2,2))
  a(1,1) = 2
  a(2,1) = 3
  a(1,2) = 7
  a(2,2) = 11
  b = conjg(transpose(a))  ! This does NOT work with gcc 4.2.1

  print *, 'original = ', a
  print *, 'WRONG conjg(transpose(a)) = ', b
end program

OBSERVATIONS:
* The bug appears regardless whether a is allocatable or not.
* The bug does not appear if b is NOT an allocatable array, e.g. a complex(kind=8),dimension(2,2) in the testcase above.
* The bug does not appear if we reverse the order of expression to: b = transpose(conjg(a)) .
Comment 1 Wirawan Purwanto 2007-08-01 20:00:15 UTC
Created attachment 14006 [details]
Failing testcase with higher-dimensional b array

I found that the result of conjg(transpose(a)) is also wrong if b is a part of higher-dimensional array (see testcase):

b(:,:,i) = conjg(transpose(a))
Comment 2 Dominique d'Humieres 2007-08-01 20:01:38 UTC
Confirmed with 4.3.0 of today on PPC Darwin8:

 original =  (  2.00000000000000     ,  0.00000000000000     ) (  3.00000000000000     ,  0.00000000000000     ) (  7.00000000000000     ,  0.00000000000000     ) (  11.0000000000000     ,  0.00000000000000     )
 WRONG conjg(transpose(a)) =  (-1.506281629186761E-154, 1.506293009711558E-154) ( 1.382417208487872E+306,-1.382417208482782E+306) (-1.506293009711557E-154, 1.506293009711558E-154) (  2.00000000000000     , -0.00000000000000     )

Comment 3 Wirawan Purwanto 2007-08-01 20:02:29 UTC
Created attachment 14007 [details]
Magically successful testcase!

And by the way, if BOTH b and a are part of higher-dimensional arrays (see testcase), then

  b(:,:,1) = conjg(transpose(a(:,:,1)))

magically works again!

Wirawan
Comment 4 Thomas Koenig 2007-08-02 05:42:20 UTC
Confirmed.

Paul, I'm putting you on the CC list because you fixed
PR 31994, the original conjg(transpose()) bug.  Maybe
you can do something about this one :-)

Thomas
Comment 5 Tobias Burnus 2007-08-02 07:56:09 UTC
The essential point is that "b" is allocatable; "a" can be a parameter or a simple dimension(2,2) array. The kind does also not matter. If one uses only conjg or only transpose there is no error.
Also, if one reverts the order, i.e. uses transpose(conjg(a)), the result is correct.

program testcase
  implicit none
  complex, allocatable :: b(:,:)
  complex,parameter :: a(1,1) = Reshape([ 2 ],[1,1])
  allocate(b(1,1))
  b = conjg(transpose(a))
!  print *, a
!  print *, b
end program
Comment 6 Tobias Burnus 2007-08-02 16:32:50 UTC
For my testcase (with parameter) I find the following in the original dump:

(*D.1371)[(S.4 + D.1380) * D.1384 + D.1387] = CONJ_EXPR <(*(complex4[0:] *) parm.1.data)[S.4 * D.1383 + D.1386]>;

Rewritten this is:
  b.data[(0+1)*1+(-1)] = CONJ_EXPR< a.data[0*1+(-2)]>
or shorter:
  b.data[0] = a.data[-2]

For A:
  S.4 * D.1383 + D.1386
  = 0 * 1      + (-2)
  = -2
where
  S.4 = 0  loop variable
  D.1383 = 1  (a.dim[0].stride)
  D.1386 = a.offset + a.dim[1].stride * S.4
         = -2       + 1               * 0
and for B:
    (S.4 + D.1380) * D.1384 + D.1387
  = (0   + 1     ) * 1      + (-1)
  =  0
where
  S.4 = 0   loop variable
  D.1380 = b.lbound = 1
  D.1384 = b.dim[0].stride = 1
  D.1387 = a.offset + a.dim[1].stride * S.3
         = (-2)     + 1               * 0
  S.3 = 0  other loop variable
Comment 7 Tobias Burnus 2007-08-02 20:56:57 UTC
  CONJ_EXPR< a.data[0*1+(-2)]>
The offset = -2 is generated in gfc_conv_expr_descriptor:

      offset = gfc_index_zero_node;
      for (n = 0; n < ndim; n++) // ndim = 2
	      start = info->start[dim]; //  start = 0;
	  tmp = gfc_conv_array_lbound (desc, n); // "tmp = 1"
          // "tmp = 0 - 1 = -1"
	  tmp = fold_build2 (MINUS_EXPR, TREE_TYPE (tmp), start, tmp);
          // tmp = tmp * 1 = -1
	  tmp = fold_build2 (MULT_EXPR, TREE_TYPE (tmp), tmp, stride);
          // offset = offset + tmp  = offset + (-1)
	  offset = fold_build2 (PLUS_EXPR, TREE_TYPE (tmp), offset, tmp);

Still, I don't quite understand how it is supposed to work.
Comment 8 kargl 2007-08-03 04:29:10 UTC
Change severity to "normal".  Yes, I know the bug is "critical"
to Wirawan, but a Fortran bug is rarely considered "critical".
Comment 9 Paul Thomas 2007-08-04 20:13:09 UTC
(In reply to comment #4)
> Confirmed.
> 
> Paul, I'm putting you on the CC list because you fixed
> PR 31994, the original conjg(transpose()) bug.  Maybe
> you can do something about this one :-)
> 

Nice try but I think that it has nothing to do with the above PR. This looks to be a scalarizer problem to me.  I am away until Wednesday - I'll have alook when I get back.

Paul
Comment 10 Paul Thomas 2007-08-04 22:13:19 UTC
The problem occurs whenever the destination array is described by and array descriptor; eg.

  complex(kind=8) :: a(2,2), b(2,2)
  call testcase (a, b)
contains
  subroutine testcase (a, b)
    complex(kind=8) :: a(:,:), b(:,:)
    a(1,1) = 2
    a(2,1) = 3
    a(1,2) = 7
    a(2,2) = 11
    b = conjg(transpose(a))

    print *, 'original = ', a
    print *, 'WRONG conjg(transpose(a)) = ', b
    end subroutine
end program

also fails and no allocatable arrays are invoved.

Paul
Comment 11 Paul Thomas 2007-08-08 16:34:49 UTC
The patch below fixes this problem and all its variants.  It regtests, except for allocatable_function_1.f90.  This runs OK but the number of free's goes up from 10 to 11; hence the failure on exepecting 10 of them.  I need to understand why this is happening - I believe that it is because an extra call to SIZE is required for on of the scalarizer loops.  Also, I want to see if there is a neater way of fixing the problem.

In summary, the problem lies with the scalarizer and can be fixed in gfc_conv_loop_setup.

Paul

Index: /svn/trunk/gcc/fortran/trans-array.c
===================================================================
*** /svn/trunk/gcc/fortran/trans-array.c        (revision 127212)
--- /svn/trunk/gcc/fortran/trans-array.c        (working copy)
*************** gfc_conv_loop_setup (gfc_loopinfo * loop
*** 3220,3226 ****
          /* Criteria for choosing a loop specifier (most important first):
             doesn't need realloc
             stride of one
!            known stride
             known lower bound
             known upper bound
           */
--- 3220,3226 ----
          /* Criteria for choosing a loop specifier (most important first):
             doesn't need realloc
             stride of one
!            known stride (preferably a variable)
             known lower bound
             known upper bound
           */
*************** gfc_conv_loop_setup (gfc_loopinfo * loop
*** 3232,3237 ****
--- 3232,3241 ----
          else if (INTEGER_CST_P (info->stride[n])
                   && !INTEGER_CST_P (specinfo->stride[n]))
            loopspec[n] = ss;
+         else if (INTEGER_CST_P (info->stride[n])
+                  && ss->expr->expr_type == EXPR_VARIABLE
+                  && loopspec[n]->expr->expr_type != EXPR_VARIABLE)
+           loopspec[n] = ss;
          else if (INTEGER_CST_P (info->start[n])
                   && !INTEGER_CST_P (specinfo->start[n]))
            loopspec[n] = ss;
Comment 12 patchapp@dberlin.org 2007-08-12 11:45:22 UTC
Subject: Bug number PR32962

A patch for this bug has been added to the patch tracker.
The mailing list url for the patch is http://gcc.gnu.org/ml/gcc-patches/2007-08/msg00774.html
Comment 13 Paul Thomas 2007-08-13 06:16:17 UTC
Subject: Bug 32962

Author: pault
Date: Mon Aug 13 06:16:03 2007
New Revision: 127391

URL: http://gcc.gnu.org/viewcvs?root=gcc&view=rev&rev=127391
Log:
2007-08-13  Paul Thomas  <pault@gcc.gnu.org>

	PR fortran/32962
	* trans-array.c (gfc_conv_array_transpose): Set the offset
	of the destination to zero if the loop is zero based.

2007-08-13  Paul Thomas  <pault@gcc.gnu.org>

	PR fortran/32962
	* gfortran.dg/transpose_1.f90: New test.


Added:
    trunk/gcc/testsuite/gfortran.dg/transpose_1.f90
Modified:
    trunk/gcc/fortran/ChangeLog
    trunk/gcc/fortran/trans-array.c
    trunk/gcc/testsuite/ChangeLog

Comment 14 Paul Thomas 2007-08-13 06:18:55 UTC
Fixed on trunk

Paul