arrays - 对于单精度数组，如何使 Fortran SUM 命令结果超过 2^24

Question

为了检查内存分配，我们用单位值填充单精度数组并使用 SUM 进行查询。和 DOT_PRODUCT命令。这些内在函数在 16777216 (= 2^24) 之后停止计数。我们如何获得这些命令来计算数十亿个元素？我们宁愿避免 DO循环。对于更高精度的阵列，这不是问题。

program allocator

  use iso_fortran_env
  implicit NONE

  integer, parameter :: sp    = selected_real_kind ( REAL32 )
  integer, parameter :: xlint = selected_int_kind  ( INT64 )

  integer ( xlint )                            :: n = 100000000
  real    ( sp ), allocatable, dimension ( : ) :: array

  integer   ( xlint )     :: alloc_status = 0
  character ( len = 255 ) :: alloc_msg = ""

!   ALLOCATE
    allocate ( array ( 1 : n ), stat = alloc_status, errmsg = alloc_msg )
    if ( alloc_status /= 0 ) print *, 'allocation error on ', n, ' elements: stat = ', alloc_status, ', errmsg = ', alloc_msg

!   POPULATE
    array = 1.0_sp
    write ( *, '( "number of elements allocated = ", g0 )' )    n
    write ( *, '( "sum of elements              = ", g0 )' )    sum ( array )
    write ( *, '( "dot product                  = ", g0, / )' ) dot_product ( array, array )

!   DEALLOCATE
    deallocate ( array, stat = alloc_status, errmsg = alloc_msg )
    if ( alloc_status /= 0 ) print *, 'deallocation error on ', n, ' elements: stat = ', alloc_status, ', errmsg = ', alloc_msg

    write ( *, '( "compiler version = ", A )' ) compiler_version()
    write ( *, '( "compiler options = ", A )' ) trim ( compiler_options() )

end program allocator

输出:

number of elements allocated = 100000000
sum of elements              = 16777216.
dot product                  = 16777216.

compiler version = GCC version 4.6.2 20111019 (prerelease)
compiler options = -fPIC -mmacosx-version-min=10.6.8 -mtune=core2

Answer 1

这是由于单精度实数的精度有限......

由于您的有效“数字”只有 24 位，因此您的分辨率为 1/2**24 = 1/16777216。换句话说，您无法将 1/1677721 的加法解析为 1，或者在您的情况下

16777216 + 1 = 16777216

为了能够解决 sum 所需的此操作和 dot_product (即使使用简单循环计算)，您(至少)还需要一点精度:

program allocator

  use iso_fortran_env
  implicit NONE

  integer, parameter :: sp    = REAL32
  integer, parameter :: xlint = INT64

  integer ( xlint )                            :: n = 100000000
  real    ( sp ), allocatable, dimension ( : ) :: array
  real    ( REAL64 )                           :: s
  integer ( xlint )                            :: i

  integer   ( xlint )     :: alloc_status = 0
  character ( len = 255 ) :: alloc_msg = ""

!   ALLOCATE
    allocate ( array ( 1 : n ), stat = alloc_status, errmsg = alloc_msg )
    if ( alloc_status /= 0 ) print *, 'allocation error on ', n, ' elements: stat = ', alloc_status, ', errmsg = ', alloc_msg

!   POPULATE
    array = 1.0_sp
    write ( *, '( "number of elements allocated = ", g0 )' )    n
    write ( *, '( "sum of elements              = ", g0 )' )    sum ( array )
    write ( *, '( "dot product                  = ", g0, / )' ) dot_product ( array, array )

    ! Calculate the sum using a double precision float
    s = real( array(1), REAL64 )
    do i=2,n
      s = s + real( array(i), REAL64 )
    enddo ! i
    write ( *, '( "sum of elements              = ", g0 )' )    s
    ! Calculate the dot product using a double precision float
    s = real( array(1), REAL64 )**2
    do i=2,n
      s = s + real( array(i), REAL64 )**2
    enddo ! i
    write ( *, '( "dot product                  = ", g0, / )' ) s

!   DEALLOCATE
    deallocate ( array, stat = alloc_status, errmsg = alloc_msg )
    if ( alloc_status /= 0 ) print *, 'deallocation error on ', n, ' elements: stat = ', alloc_status, ', errmsg = ', alloc_msg

    write ( *, '( "compiler version = ", A )' ) compiler_version()
    write ( *, '( "compiler options = ", A )' ) trim ( compiler_options() )

end program allocator

输出:

number of elements allocated = 100000000
sum of elements              = 16777216.0
dot product                  = 16777216.0

sum of elements              = 100000000.00000000
dot product                  = 100000000.00000000

compiler version = GCC version 4.8.4 20140605 (prerelease)
compiler options = -cpp -iprefix /home/elias/opt/sde/bin/../lib/gcc/x86_64-unknown-linux-gnu/4.8.4/ -mtune=generic -march=x86-64 -O0 -Wall -Wextra