Choosing the default summation algorithm for Flang's SUM intrinsic

Subprojects Flang
1

I am working on issue #176282, which concerns differences in the numerical result of the SUM intrinsic between -O0 and -O1.

Current behavior

At present, Flang uses different implementations of SUM depending on the optimization level.

At -O0, SUM is lowered to a call to a Flang runtime function. For floating-point arrays, the runtime implementation uses Kahan summation.

At -O1 and higher, SUM is generally open-coded as a loop, using straightforward sequential accumulation rather than Kahan summation.

Consequently, changing the optimization level can change not only the generated code but also the summation algorithm itself.

I think it would be preferable for Flang to use a consistent algorithm across optimization levels, so I have been investigating which implementation should be used by default.

Experiment

I created an experimental change that makes the open-coded implementation used at -O1 and higher perform Kahan summation instead of straightforward accumulation.

The measurements were performed using Flang built from llvm-project commit 511286d5cffd, with this experimental change applied on top.

First, I measured the execution time of a small test program below whose runtime is dominated by calls to the SUM intrinsic. With the experimental implementation, the program was approximately four to five times slower than with straightforward accumulation.

program main
 integer,parameter :: N=76
 real(kind=4),parameter :: val=1._4
 real(kind=4),dimension(1:N,1:N,1:N) :: a3

 a3 = val
 do i=1, N-1
    do j=1, N
       do k=1, N
          a3(i+1,j,k) = a3(i,j,k) + val
       enddo
    enddo
 enddo
 write(1002) a3
 rewind 1002
 call sub
end program main

subroutine sub
  integer,parameter :: N=76
  real(kind=4),dimension(1:N,1:N,1:N) :: a
  real(kind=4) :: result
  integer :: c1, c2, rate
  real(8) :: elapsed

  call system_clock(count_rate=rate)

  read(1002) a
  call system_clock(c1)
  do i = 1, 5000
    result = sum(a)
  enddo
  call system_clock(c2)

  elapsed = real(c2 - c1, kind=8) / real(rate, kind=8)
  print *,'   SUM: ',result
  print *,'  Time: ',elapsed
end subroutine sub

Next, I measured the performance impact on SPEC CPU 2017 using an NVIDIA Grace CPU. I selected only benchmarks that use the SUM intrinsic internally.

Category Benchmark speedup
SPECrate INT 548.exchange2_r -0.022%
SPECrate FP 521.wrf_r -0.158%
527.cam4_r -0.437%
SPECspeed INT 648.exchange2_s 0.088%
SPECspeed FP 621.wrf_s 0.500%
627.cam4_s -11.830%
628.pop2_s 0.954%

Most of the benchmarks showed little or no performance degradation. However, the performance of 627.cam4_s regressed by approximately 11.83%. Further analysis indicated that the main cause of the regression was the increased execution time of the SUM operations in the benchmark.

Comparison with Other Compilers

I also examined the behavior and generated assembly of other Fortran compilers, including gfortran and ifx. In the cases I tested, both compilers consistently used straightforward accumulation for the SUM intrinsic at every optimization level from -O0 through -O3.

Possible approach

Kahan summation generally provides better numerical accuracy than straightforward accumulation, so it is a reasonable algorithm for implementing SUM when accuracy is the primary concern.

However, an overhead of more than 10% in an application that uses SUM extensively may be unacceptable to many users. For this reason, I am hesitant to make Kahan summation the default implementation at all optimization levels.

One possible approach would be:

  • Use straightforward accumulation as the default implementation of SUM, consistently across optimization levels.
  • Provide a compiler option that requests a compensated summation algorithm, such as Kahan summation, for users who prefer improved numerical accuracy and are willing to accept the performance cost.

This would eliminate the current algorithmic difference between -O0 and optimized builds while still giving users a way to request a more accurate implementation.

Does this seem like a reasonable direction? I would appreciate feedback on both the choice of the default algorithm and the idea of providing a compiler option for compensated summation.

2

Is the inline Kahan summation vectorized at higher optimization levels? If the compiler does not vectorize it, this may cause the slowdown. Could you please update the inline code to make it vectorizable? This change may improve the performance of the inline Kahan summation.

It is a good plan to let users choose the summation algorithm. As you said, this approach gives users the final choice. We can separately decide on the default algorithm. We can also decide whether to use different algorithms at different optimization levels.

1 Like
3

As it was said at today’s meeting, there might be some interference between Kahan summation and FastMath. I think some FastMath flags may allow LLVM to optimize Kahan summation into a regular summation. So if we let users force Kahan summation in the inlined SUM code, we have to define how it works with separately specified FastMath behavior (e.g. we can drop certain FastMath flags on the inlined code, etc.).

4

I agree. I have done a similar thing for MOD on REAL types, where the default code calls the precise implementation in the Fortran runtime and -ffast-math relaxes that to a less precise code-gen’ed version. We could do something very similar: use Kahan summation by default and have -ffast-sum to use a less accurate but faster algorithm (and maybe have -ffast-math set that that flag, too)

5

I have found this paper (from the collection at the Manchester University Numerical Linear Algebra Group) very helpful. It does not address reproducibility issues, only accuracy versus performance. In the light of the paper, I had prepared a timing harness to test ideas.

time_compensated.f90.txt (12.5 KB)

The kernels I am timing (nanoseconds per element are displayed, with accuracy error below them, for a variety of array sizes) are

  • Rx[1,4,8,16,32,128], performing the recursive sum with as many accumulators as indicated in the decimal number (the partial accumulators are added up with intrinsic SUM).
  • Cx[1,4,8,16,32,128], performing the compensated (“Kahan”) sum, again with as many accumulators as indicated in the decimal number (the partial accumulators are added up with Cx1).
  • COMBI, performing the recursive sum with as many accumulators as indicated in the decimal number (the partial accumulators are added up with Cx1).
  • SUM, performing the sum with intrinsic SUM.

You can verify that, in most cases, there is a Fortran kernel that outperforms the SUM intrinsic (in both speed and accuracy, when compiled/linked with -O3). This does not look good for SUM!

When adding -ffast-math, SUM is good at performance, but note that for an extra 10% of runtime, one can gain 2 extra decimal places of accuracy by using Cx16 or COMBI.

6

SUM intrinsic can handle non-contiguous arrays. You can imagine that in the case where each element of the array is hitting a different cache line, you probably have enough time to do the extra flops for compensated summation while waiting for the cache to update and you should get extra accuracy for “free”.

(We should probably be using non-temporal load instructions in such a case, since polluting the cache will give us no advantage).