Ошибка SIGFPE с обработкой gfortran 4.8.5

Question

Я использую программное обеспечение для вычислительной гидродинамики, которое скомпилировано с версией gfortran 4.8.5 на Ubuntu 16.04 LTS. Программное обеспечение может быть скомпилировано с одинарной или двойной точностью и с опцией оптимизации -O3. Поскольку у меня нет необходимых вычислительных ресурсов для запуска программного обеспечения CFD с двойной точностью, я компилирую его с одинарной точностью и следующими параметрами

ffpe-trap=invalid,zero,overflow

Я получаю ошибку SIGFPE в строке кода, содержащей функцию asin-

 INTEGER, PARAMETER :: sp = SELECTED_REAL_KIND( 6, 37) !< single precision
 INTEGER, PARAMETER :: wp = sp
 REAL(KIND=wp) zsm(:,:)

ela(i,j) = ASIN(zsm(ip,jp))

Другими словами, функция обратного греха и этот код являются частью дважды вложенного цикла FOR с индексами jp и ip. В настоящее время сотрудники программного обеспечения не могут мне помочь по ряду других причин, и поэтому я пытаюсь отладить это самостоятельно. Ошибка SIGFPE наблюдается только при компиляции с одинарной точностью, а не при компиляции с двойной точностью.

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

  if (zsm(ip,jp) >= 1.0 .or. zsm(ip,jp) <= -1.0) then
       print *,zsm(ip,jp),ip,jp
  end if

РЕДАКТИРОВАТЬ Я посмотрел на этот ответ Неожиданное поведение asin в R, и мне интересно, смогу ли я сделать что-то подобное в fortran, то есть с помощью функции max. Если оно опускается ниже -1 или больше 1, тогда закруглите его надлежащим образом. Как я могу сделать это с помощью gfortran, используя функцию max?

На моем рабочем столе следующая программа выполняется без проблем (т.е. она имеет возможность правильно обрабатывать подписанные нули), и поэтому я предполагаю, что ошибка SIGFPE возникает с аргументом больше 1 или меньше -1.

 program testa

 real a,x

 x = -0.0000
 a = asin(x)
 print *,a
 end program testa

Answer 1

У нас есть функции min и max в Fortran, поэтому я думаю, что мы можем использовать тот же метод, что и на связанной странице, т.е. asin( max(-1.0,min(1.0,x) ), Я попробовал следующий тест с gfortran-4.8 и 7.1:

program main
    implicit none
    integer, parameter :: sp = selected_real_kind( 6, 37 )
    integer, parameter :: wp = sp
    ! integer, parameter :: wp = kind( 0.0 )
    ! integer, parameter :: wp = kind( 0.0d0 )
    real(wp) :: x, a

    print *, "Input x"
    read(*,*) x

    print *, "x =", x
    print *, "equal to 1 ? :", x == 1.0_wp
    print *, asin( x )
    print *, asin( max( -1.0_wp, min( 1.0_wp, x ) ) )
end

который дает с wp = sp (или же wp = kind(0.0) на моем компьютере)

$ ./a.out
 Input x
1.00000001
 x =   1.00000000    
 equal to 1 ? : T
   1.57079625             (<- 1.5707964 for gfortran-4.8)
   1.57079625    

$ ./a.out
 Input x
1.0000001
 x =   1.00000012    
 equal to 1 ? : F
              NaN
   1.57079625 

и с wp = kind(0.0d0)

$ ./a.out
 Input x
1.0000000000000001
 x =   1.0000000000000000     
 equal to 1 ? : T
   1.5707963267948966
   1.5707963267948966     

$ ./a.out
 Input x
1.000000000000001     
 x =   1.0000000000000011     
 equal to 1 ? : F
                       NaN
   1.5707963267948966 

---

Если необходимо изменить много asin(x) и программа опирается на препроцессор C или Fortran, может быть удобно определить некоторый макрос, например

#define clamp(x) max(-1.0_wp,min(1.0_wp,x))

и использовать его как asin( clamp(x) ), Если мы хотим удалить такую ​​модификацию, мы можем просто изменить определение зажима () как #define clamp(x) (x), Другим подходом может быть определение некоторых asin2(x) функция, которая ограничивает x в [-1,1] и заменить встроенный asin от asin2 (либо как макрос, либо как функция Фортрана).