Generation of array temporaries for pure functions in fortran

user_4954079 0
Viewed 112

My question is about some puzzling behavior of gfortran (gcc version 8.1.0 (GCC)) in treating temporary arrays. The code below is compiled with -O3 -Warray-temporaries.

Consider the following possibility to construct sum of direct products of matrices:

program Pexam1
  use Maux, only: kron
  double precision:: A(3, 3), B(3, 3)
  double precision:: C(5, 5), D(5, 5)
  double precision:: rAC(15, 15), rBD(15, 15), r1(15, 15), r2(15, 15)
  
  r1 = -kron(3, A, 5, C) - kron(3, B, 5, D)
!                             1
! Warning: Creating array temporary at (1)

  rAC = -kron(3, A, 5, C)
!           1
! Warning: Creating array temporary at (1)   

  r2 = rAC - kron(3, B, 5, D)
end program Pexam1

There are two warnings, but I do not expect the second one (in rAC = -kron(3, A, 5, C)). Why? Simply in the first line r1 = -kron(3, A, 5, C) - kron(3, B, 5, D) the identical first part is evaluated just fine without a temporary array. What is the reason for this inconsistency? What should one keep in mind to avoid unnecessary creation of temporary arrays?

The module Maux is as follows

module Maux
 contains
 pure function kron(nA, A, nB, B) result (C)
    ! // direct product of matrices 
    integer, intent(in)         :: nA, nB
    double precision, intent(in)  :: A(nA, nA), B(nB, nB)
    double precision              :: C(nA*nB, nA*nB)
    integer :: iA, jA, iB, jB

    forall (iA=1:nA, jA=1:nA, iB=1:nB, jB=1:nB)&
         C(iA + (iB-1) * nA, jA + (jB-1) * nA) = A(iA, jA) * B(iB, jB)

  end function kron
end module Maux
gfortran fortran arrays
1 Answers

If you build this program, you'll notice an array temporary is only created in case you're evaluating kron with the minus sign:

 
 program Pexam1
  use Maux, only: kron
  double precision:: A(3, 3), B(3, 3)
  double precision:: C(5, 5), D(5, 5)
  double precision:: rAC(15, 15), rBD(15, 15), r1(15, 15), r2(15, 15)
  

  rAC = kron(3, A, 5, C) ! no temporary

  r2  = -kron(3, B, 5, D) ! temporary
end program Pexam1

Don't forget that every assignment in fortran (=) is a separate operation than what happens to the r.h.s. of the equal sign. So when kron is negative, the compiler needs to:

  1. evaluate kron(3,B,5,D), put it in a temporary (let's call it tmp)
  2. evaluate -tmp
  3. assign the result of this temporary operation (could have any type/kind) to r2

Apparently, gFortran is good enough that when you have step 1. only (positive assignment), the result of the operation is copied directly to rAC without any temporaries (I guess because rAC has the same rank, type, shape as the result of the call to kron).

user_47028700
Related