Fortran纯函数中临时数组生成的异常警告问题

使用gfortran 8.1.0编译Fortran代码时，开启-O3 -Warray-temporaries选项遇到了临时数组处理的费解行为。以下是构造矩阵直积和的代码示例：

主程序代码 #

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

编译后出现两个警告，但rAC = -kron(3, A, 5, C)处的警告不符合预期：第一行r1 = -kron(3, A, 5, C) - kron(3, B, 5, D)中完全相同的-kron(3, A, 5, C)部分并未触发临时数组警告，两者行为不一致。

对应的Maux模块代码 #

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 8.1.0优化器的行为限制导致的：

• 在复合表达式r1 = -kron(...) - kron(...)中，编译器可以将第一个-kron(...)的计算结果直接写入目标数组r1的内存空间，无需额外创建临时数组——因为后续的减法操作可以直接在r1的内存上进行，优化器能识别这种“原地计算”的模式。
• 而在单独的赋值语句rAC = -kron(...)中，-kron(...)是一个带一元负号的数组表达式。gfortran 8.1.0的优化器无法将“函数返回数组”和“一元负号操作”合并为直接写入目标数组的操作，必须先创建临时数组存储kron的返回结果，再对临时数组取负后赋值给rAC，因此触发警告。

避免不必要临时数组的方法 #

1. 拆分赋值操作
   将带负号的赋值拆分为两步：先直接赋值kron的结果，再原地取负。这样编译器可以直接将kron的结果写入目标数组，无需临时数组：

   rAC = kron(3, A, 5, C)
   rAC = -rAC
   

2. 升级gfortran版本
   gfortran 9及以后的版本优化器对这类场景做了改进，能够识别并优化单独赋值中的一元负号操作，避免临时数组的创建。

3. 修改函数接口，添加原地计算选项
   可以给kron函数添加一个可选的输出参数，直接将计算结果写入目标数组（甚至支持带符号的输出），从根源上避免临时数组：

   pure subroutine kron_inplace(nA, A, nB, B, C, sign)
     integer, intent(in)         :: nA, nB
     double precision, intent(in)  :: A(nA, nA), B(nB, nB)
     double precision, intent(out) :: C(nA*nB, nA*nB)
     integer, intent(in), optional :: sign
     integer :: iA, jA, iB, jB
     real(8) :: s
   
     s = 1.0d0
     if (present(sign)) s = sign(s)
   
     forall (iA=1:nA, jA=1:nA, iB=1:nB, jB=1:nB)&
          C(iA + (iB-1) * nA, jA + (jB-1) * nA) = s * A(iA, jA) * B(iB, jB)
   end subroutine kron_inplace
   

   调用时直接传入目标数组和符号：

   call kron_inplace(3, A, 5, C, rAC, sign=-1)
   

内容的提问来源于stack exchange，提问作者yarchik