优化击败稳健性措施
Optimization defeats robustness measures
我目前正在研究用于数组求和的稳健方法,并实现了 Shewchuk 在 "Adaptive Precision Floating-Point Arithmetic and Fast Robust Geometric Predicates" 中发布的算法。虽然实施的算法在 gfortran
中按预期工作,但 ifort
优化了对策。
为了提供一些背景信息,这是我的代码:
module test_mod
contains
function shewchukSum( array ) result(res)
implicit none
real,intent(in) :: array(:)
real :: res
integer :: xIdx, yIdx, i, nPartials
real :: partials(100), hi, lo, x, y
nPartials = 0
do xIdx=1,size(array)
i = 0
x = array(xIdx)
! Calculate the partial sums
do yIdx=1,nPartials
y = partials(yIdx)
hi = x + y
if ( abs(x) < abs(y) ) then
lo = x - (hi - y)
else
lo = y - (hi - x)
endif
x = hi
! If a round-off error occured, store it. Exact comparison intended
if ( lo == 0. ) cycle
i = i + 1 ; partials(i) = lo
enddo ! yIdx
nPartials = i + 1 ; partials( nPartials ) = x
enddo ! xIdx
res = sum( partials(:nPartials) )
end function
end module
调用测试程序为
program test
use test_mod
implicit none
print *, sum([1.e0, 1.e16, 1.e0, -1.e16])
print *,shewchukSum([1.e0, 1.e16, 1.e0, -1.e16])
end program
使用 gfortran
的编译为所有优化级别生成正确的结果:
./a.out
0.00000000
2.00000000
然而,ifort
会为 -O0
:
以上的所有优化生成零
./a.out
0.00000000
0.00000000
我尝试调试代码并深入到汇编级别,发现 ifort
正在优化 lo
的计算和 if ( lo == 0. ) cycle
之后的操作。
是否有可能强制ifort
对所有优化级别执行完整操作?此加法是计算的关键部分,我希望它尽快 运行。
相比之下,-O2
处的 gfortran
执行此代码的速度大约是 -O0
处的 ifort
的八到十倍(针对长度 >100k 的数组进行测量)。
当谈到浮点运算时,ifort 的默认值通常是为了性能而不是严格的正确性。
有许多选项可以控制浮点行为。使用 ifort 16 和选项 -assume protect_parens
即使在更高的优化级别下我也能得到预期的行为。
此外,还有选项 -fp-model precise -fp-model source
(后者表示 -assume protect_parens
),您可能也会感兴趣。 -fp-model
的默认值是 fast=1
其中
allows value-unsafe optimizations
当然,这些可能会对性能产生影响,因此围绕浮点行为的其他选项也值得考虑。
可以在 Intel publication 中找到更多详细信息。
我目前正在研究用于数组求和的稳健方法,并实现了 Shewchuk 在 "Adaptive Precision Floating-Point Arithmetic and Fast Robust Geometric Predicates" 中发布的算法。虽然实施的算法在 gfortran
中按预期工作,但 ifort
优化了对策。
为了提供一些背景信息,这是我的代码:
module test_mod
contains
function shewchukSum( array ) result(res)
implicit none
real,intent(in) :: array(:)
real :: res
integer :: xIdx, yIdx, i, nPartials
real :: partials(100), hi, lo, x, y
nPartials = 0
do xIdx=1,size(array)
i = 0
x = array(xIdx)
! Calculate the partial sums
do yIdx=1,nPartials
y = partials(yIdx)
hi = x + y
if ( abs(x) < abs(y) ) then
lo = x - (hi - y)
else
lo = y - (hi - x)
endif
x = hi
! If a round-off error occured, store it. Exact comparison intended
if ( lo == 0. ) cycle
i = i + 1 ; partials(i) = lo
enddo ! yIdx
nPartials = i + 1 ; partials( nPartials ) = x
enddo ! xIdx
res = sum( partials(:nPartials) )
end function
end module
调用测试程序为
program test
use test_mod
implicit none
print *, sum([1.e0, 1.e16, 1.e0, -1.e16])
print *,shewchukSum([1.e0, 1.e16, 1.e0, -1.e16])
end program
使用 gfortran
的编译为所有优化级别生成正确的结果:
./a.out
0.00000000
2.00000000
然而,ifort
会为 -O0
:
./a.out
0.00000000
0.00000000
我尝试调试代码并深入到汇编级别,发现 ifort
正在优化 lo
的计算和 if ( lo == 0. ) cycle
之后的操作。
是否有可能强制ifort
对所有优化级别执行完整操作?此加法是计算的关键部分,我希望它尽快 运行。
相比之下,-O2
处的 gfortran
执行此代码的速度大约是 -O0
处的 ifort
的八到十倍(针对长度 >100k 的数组进行测量)。
当谈到浮点运算时,ifort 的默认值通常是为了性能而不是严格的正确性。
有许多选项可以控制浮点行为。使用 ifort 16 和选项 -assume protect_parens
即使在更高的优化级别下我也能得到预期的行为。
此外,还有选项 -fp-model precise -fp-model source
(后者表示 -assume protect_parens
),您可能也会感兴趣。 -fp-model
的默认值是 fast=1
其中
allows value-unsafe optimizations
当然,这些可能会对性能产生影响,因此围绕浮点行为的其他选项也值得考虑。
可以在 Intel publication 中找到更多详细信息。