Как получить результат команды Fortran SUM, превышающий 2 ^ 24 для массивов с одинарной точностью

Чтобы проверить распределение памяти, мы заполняем массивы одинарной точности единичными значениями и опрашиваем их с помощью команд 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

arrays sum fortran dot-product fortran95
person dantopa    schedule 14.08.2014    source источник
comment
Вы, наверное, знаете, почему это число является результатом, но почему вы против повторных циклов? Мне не ясно, спрашиваете ли вы, как получить большее число, или важно проверить распределение. Не могли бы вы уточнить (и если последнее, то добавить то, что вы хотите проверить)?   -  person francescalus    schedule 15.08.2014
comment
@francescalus Почему это число является результатом?   -  person PetrH    schedule 15.08.2014
comment
Точно так же, как ваше n не подходит для INT32, сумма массива не подходит для REAL32... Вы можете использовать REAL(SUM(DBLE(массив))), если хотите, чтобы результат был числом одинарной точности.   -  person Etienne Pellegrini    schedule 15.08.2014
comment
Если вы скомпилируете с -ffpe-trap=precision, вы увидите, что вы получаете ошибку времени выполнения Program received signal SIGFPE: Floating-point exception - erroneous arithmetic operation., потому что вы переполняете вещественное число с одинарной точностью в своем sum   -  person casey    schedule 15.08.2014
comment
@francescalus Наше приложение предназначено для диагностического инструмента на наших фермах Cray. Цель состоит в том, чтобы использовать естественные параллельные инструменты в Fortran, чтобы избежать использования MPI, OpenMP и т. д.   -  person dantopa    schedule 17.08.2014
comment
@Этьен Пеллегрини. Умная мысль на самом деле. Спасибо.   -  person dantopa    schedule 17.08.2014
comment
@кейси. Очень хорошо. Спасибо!   -  person dantopa    schedule 17.08.2014
comment
@francescalus Мы просматриваем все типы данных (целые, действительные, комплексные) и все доступные точности и захватываем все большие и большие куски памяти. Мы отслеживаем все шаги в процессе выделения-заполнения-освобождения и сообщаем сумму и скалярные произведения в диагностическую сводку. Мы собираем информацию о времени и проверяем производительность по мере того, как новые машины подключаются к сети.   -  person dantopa    schedule 17.08.2014
comment
Вы четко знаете, какие тесты хотите запустить, и у вас есть предложения, как это лучше сделать. Однако нужно ли вам делать суммирование/точечные произведения или будет достаточно простой проверки, что все элементы равны 1 (и т. д.)?   -  person francescalus    schedule 17.08.2014
comment
В качестве второстепенного момента: когда вы дойдете до очень большого n, вы хотите указать целочисленный литерал как тип, отличный от значения по умолчанию (100000000000000_xlint).   -  person francescalus    schedule 17.08.2014
comment
@francescalus Я считаю, что объявление в коде и метод, который вы показываете, эквивалентны. Что отличает эти методы? По поводу здравой идеи протестировать все элементы = 1: Попытка перетащить толпу ФОРТРАН 77 в век XXI.   -  person dantopa    schedule 18.08.2014
comment
@dantopa Без завершающего _xlint выражение имеет целочисленный вид по умолчанию, независимо от того факта, что оно является инициализатором для целого числа вида xlint.   -  person francescalus    schedule 18.08.2014
comment
Я не говорю, что вы ошибаетесь, когда смотрите на суммы и т. Д., Но с философской точки зрения вы исходите из того, работает ли мое распределение памяти? правильно ли мой численный анализ? ситуация.   -  person francescalus    schedule 18.08.2014
comment
@francescalus Спасибо за пристальное внимание к декларации. Что же касается философской точки зрения, то вы дали точное определение.   -  person dantopa    schedule 18.08.2014

Ответы (1)


arrow_upward
3
arrow_downward

Это связано с ограниченной точностью с одинарной точностью...

Поскольку у вас есть только 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
person Alexander Vogt    schedule 14.08.2014
comment
Действительно, но с простыми циклами можно использовать лучший алгоритм суммирования или, по крайней мере, дополнительные биты точности не требуют временного массивного массива. - person francescalus; 15.08.2014
comment
Верно, другими словами, 1/2**24 — это машинный эпсилон с одинарной точностью. Понятно. - person PetrH; 15.08.2014
comment
Но для меня его код работал с флагами ifort 14 и по умолчанию. Какое волшебство тогда проделал компилятор? - person PetrH; 15.08.2014
comment
@PetrH Возможно, промежуточные вычисления не с одинарной точностью; не (1.+(1.+(1.+(1.+.....))); оптимизация; что-то другое. - person francescalus; 15.08.2014
comment
@PetrH, что бы он ни делал, вы не можете рассчитывать на его поведение. Стандарт предусматривает, что результат sum() имеет тот же тип и вид, что и массив, который он сокращает (13.7.161 p4). - person casey; 15.08.2014
comment
@casey Да, но с real(real32) можно довольно близко подойти к 1e8. Насколько близко вы подходите, зависит от качества реализации - и, как вы говорите, не определено. - person francescalus; 15.08.2014
comment
@francescalus Я вижу, урок в том, что компилятор может выбирать, как выполнять промежуточные вычисления. Так что это может даже зависеть от таких вещей, как аппаратное обеспечение. Или параллельное вычисление суммы тоже может сделать это. - person PetrH; 15.08.2014