From mboxrd@z Thu Jan  1 00:00:00 1970
Return-Path: <gcc-bugs-return-326695-listarch-gcc-bugs=gcc.gnu.org@gcc.gnu.org>
Received: (qmail 22488 invoked by alias); 2 Aug 2010 13:39:27 -0000
Received: (qmail 22423 invoked by uid 48); 2 Aug 2010 13:39:15 -0000
Date: Mon, 02 Aug 2010 13:39:00 -0000
Message-ID: <20100802133915.22422.qmail@sourceware.org>
X-Bugzilla-Reason: CC
References: <bug-45159-13404@http.gcc.gnu.org/bugzilla/>
Subject: [Bug fortran/45159] Unneccessary temporaries
In-Reply-To: <bug-45159-13404@http.gcc.gnu.org/bugzilla/>
Reply-To: gcc-bugzilla@gcc.gnu.org
To: gcc-bugs@gcc.gnu.org
From: "dominiq at lps dot ens dot fr" <gcc-bugzilla@gcc.gnu.org>
Mailing-List: contact gcc-bugs-help@gcc.gnu.org; run by ezmlm
Precedence: bulk
List-Id: <gcc-bugs.gcc.gnu.org>
List-Archive: <http://gcc.gnu.org/ml/gcc-bugs/>
List-Post: <mailto:gcc-bugs@gcc.gnu.org>
List-Help: <mailto:gcc-bugs-help@gcc.gnu.org>
Sender: gcc-bugs-owner@gcc.gnu.org
X-SW-Source: 2010-08/txt/msg00068.txt.bz2



------- Comment #7 from dominiq at lps dot ens dot fr  2010-08-02 13:39 -------
> I fail to see why a scalar temporary is not enough:
>
> do j = 1, N
>   do i = 1, N
>     tmp = 0.0
>     do k = 1, N
 >      tmp = a(i,k)*b(k,j)
>     end do
>     a(i,j) = tmp
>   end do
> end do
>
> Note: it won't work with a scalar if one reverses the i and the j loop.

Well, try:


real :: a(3,3), b(3,3), c(3,3), tmp
integer :: i, j, k, N=3
a = reshape([(i,i=1,9)],[3,3])
b = reshape([(10-i,i=1,9)],[3,3])
c = matmul(a, b)
do j = 1, N
  do i = 1, N
    tmp = 0.0
    do k = 1, N
      tmp = tmp + a(i,k)*b(k,j)
    end do
    a(i,j) = tmp
  end do
end do
print *, a
print *, c
if(any(a/=c)) call abort()
print *, 'OK'
end

The problem is that in order to compute a(1,2) you need a(1,1)*b(1,2), so you
cannot have already written a(1,1). Note that the effective way to compute the
matrix product on modern CPUs is to exchange the i and k loops and use a rank 1
tmp:

real :: a(3,3), b(3,3), c(3,3), tmp(3)
integer :: i, j, k, N=3
a = reshape([(i,i=1,9)],[3,3])
b = reshape([(10-i,i=1,9)],[3,3])
c = matmul(a, b)
do j = 1, N
  tmp = 0.0
  do k = 1, N
    do i = 1, N
      tmp(i) = tmp(i) + a(i,k)*b(k,j)
    end do
  end do
  b(:,j) = tmp
end do
if(any(b/=c)) call abort()
print *, 'OK'
end

and this work with a rank 1 temporary.

If one want to enter this kind of reduced temporaries, another candidate is
a=transpose(a) which requires a scalar temporary only.


-- 


http://gcc.gnu.org/bugzilla/show_bug.cgi?id=45159