Как избежать занижения при умножении действительных чисел?

Я хотел бы знать, если есть правильный способ умножения или квадрата числа с плавающей запятой (или удваивать) вместе без ошибки занижения при компиляции моего кода на Фортране, как

gfortran -ffpe-trap=invalid,zero,overflow,underflow ...

Я знаю что underflow опция не всегда хорошая идея, но мне интересно, возможно ли сделать умножение с этой опцией. На самом деле, в следующем примере я знаю, что может произойти потеря значения, но, возможно, я не знаю другого случая в моих кодах. Вот почему я хотел бы сохранить эту опцию, если это возможно.

Вот пример, где я вычисляю вектор u для каждого индекса x,y матрицы; 2 значения, составляющие векторы тезисов, находятся между 0 и 1. Затем я вычисляю квадрат его нормы.

Очень логично, что из-за этой квадратной операции у меня будут заниженные значения. Поскольку эти очень маленькие значения можно считать нулем для меня. Есть ли способ не underflow лучше, чем с помощью if сравнение?

implicit none
double  :: u(100,100,2), uSqr(100,100)
integer :: x,y

DO x= 1, 100
    DO y = 1, 100
                CALL Poisin( u(x,y,:), x, y )
    ENDDO
ENDDO

uSqr = u(:,:,1)*u(:,:,1) + u(:,:,2) * u(:,:,2) ! where comes the underflow errors
Источник

2 ответа

Решение

У вас есть ответ, который рассматривает конкретный способ избежать чрезмерного снижения цен при определенных обстоятельствах. Это использует hypot функция. Это частично ответ: если вы хотите избежать недополнения, возможно, есть способ переписать алгоритм, чтобы избежать этого.

Для более общих случаев (например, в этом вопросе), где требуется точный контроль флагов исключений, это не подходит. Однако компиляторы часто предлагают интерфейсы для процедур обработки исключений.

Один из переносимых способов сделать это - использовать средство IEEE на Fortran 2003. [Если вы используете gfortran, вам понадобится как минимум версия 5.0, но есть аналогичные доступные для компилятора способы.]

Fortran определяет исключения IEEE и флаги. Флаг может быть тихим или сигнальным. То, что вы хотите, это для части, где недополнение не является полезной диагностикой, чтобы не влиять на состояние флага недопущения после этого вычисления.

Флаг известен как IEEE_UNDERFLOW, Мы можем запросить и установить его статус с помощью вызовов подпрограмм. IEEE_GET_FLAG(IEEE_UNDERFLOW, value) а также IEEE_SET_FLAG(IEEE_UNDERFLOW, value), Если мы ожидаем, но не заботимся о недостаточном объеме, мы также хотим убедиться, что исключение не является остановкой. Подпрограмма IEEE_SET_HALTING_MODE(IEEE_UNDERFLOW, value) контролирует этот режим.

Итак, аннотированный пример.

  use, intrinsic :: ieee_arithmetic, only : IEEE_SELECTED_REAL_KIND
  use, intrinsic :: ieee_exceptions

  implicit none

  ! We want an IEEE kind, but this doesn't ensure support for underflow control
  integer, parameter :: rk=IEEE_SELECTED_REAL_KIND(6, 70)

  ! State preservation/restoration
  logical should_halt, was_flagged

  real(rk) x

  ! Get the original halting mode and signal state 
  call ieee_get_halting_mode(IEEE_UNDERFLOW, should_halt)
  call ieee_get_flag(IEEE_UNDERFLOW, was_flagged)

  ! Ensure we aren't going to halt on underflow
  call ieee_set_halting_mode(IEEE_UNDERFLOW, .FALSE.)

  ! The irrelevant computation   
  x=TINY(x)
  x=x**2

  ! And restore our old state
  call ieee_set_halting_mode(IEEE_UNDERFLOW, should_halt)
  call ieee_set_flag(IEEE_UNDERFLOW, was_flagged)

end program
Источник

Если это просто для того, чтобы избежать потери, возможно, вы хотите использовать hypot или тебе действительно нужен квадрат гипотенузы?

Реализация hypot следует избегать проблем переполнения sqrt(x**2+y**2),

Источник
Другие вопросы по тегам