開啟fortran浮點數異常處理

fortran對於浮點數異常默認是忽略的，如果要開啟，有兩種方法，一種可以通過調用c函數去開啟這部分異常處理，並且需要將函數打包成庫，鏈接到fortran中，另一個是通過編譯選項設置。

第一種方法的代碼：

trapfpe.c

/*
* ref:
* https://gcc.gnu.org/onlinedocs/gcc-3.3.6/g77/Floating_002dpoint-Exception-Handling.html
* https://riptutorial.com/fortran/example/7149/calling-c-from-fortran
* https://stackoverflow.com/questions/17845931/calling-c-function-subroutine-in-fortran-code
* https://stackoverflow.com/questions/1202494/why-doesnt-attribute-constructor-work-in-a-static-library
*/

#define _GNU_SOURCE 1
#include <fenv.h>
#include <stdio.h>

// 在main函數開始前調用
// fortran對於浮點數異常默認是忽略的，這里只能通過調用c函數去開啟這部分異常處理
// 並且需要將函數打包成庫，鏈接到fortran中
static void __attribute__ ((constructor))
trapfpe ()
{
  /* Enable some exceptions.  At startup all exceptions are masked.  */
  printf("trapfpe: enabling exceptions\n");
  feenableexcept (FE_INVALID|FE_DIVBYZERO|FE_OVERFLOW);
}

void fun()
{
    printf("func\n");
}

// static void before(void) __attribute__((constructor));
// static void after(void) __attribute__((destructor));

// static void before()
// {
//     printf("before main\n");
// }

// static void after(void)
// {

//     printf("after main\n");

// }

main.f90

program main
    implicit none

! 必須顯式調用靜態庫中一個的函數，才能讓trapefpe生效
! 因此這里定義一個接口，讓fortran調用c函數
INTERFACE
  SUBROUTINE fun() BIND(C)
    USE, INTRINSIC :: ISO_C_BINDING, ONLY: C_INT
    IMPLICIT NONE
  END SUBROUTINE fun
END INTERFACE

    double precision :: a
    double precision :: b
    double precision :: c

    a = 1.0
    b = 0.0
    c = a / b
    print *, c
    call fun()
end program

編譯方法：

# 生成trapfe.o
gcc -c trapfpe.c
# 生成靜態庫文件libtrapfe.a
ar -r libtrapfpe.a trapfpe.o
# 將.目錄下的trapfe.a作為靜態庫，鏈接到main中
gfortran -g main.f90 -L. -ltrapfpe

運行結果：

➜  test ./a.out 
trapfpe: enabling exceptions

Program received signal SIGFPE: Floating-point exception - erroneous arithmetic operation.

Backtrace for this error:
#0  0x7fd098838d21 in ???
#1  0x7fd098837ef5 in ???
#2  0x7fd09851a0bf in ???
#3  0x5619e64951bf in MAIN__
        at /mnt/f/programing/test/main.f90:19
#4  0x5619e649526b in main
        at /mnt/f/programing/test/main.f90:22
[1]    990 floating point exception  ./a.out

可以看到，第19行出現了浮點數異常，因為這里做了除0操作。

第二種方法的代碼

fortran

! https://gcc.gnu.org/onlinedocs/gfortran/Debugging-Options.html
program main 
    implicit none 
    real a,b,c
    b = 1.0
    c = 0.0
    a = b / c 
    write(*,*) a
end program main

編譯：

gfortran -g main2.f90 -ffpe-trap=invalid,zero,overflow

輸出：

➜  fortran-exception-handle ./a.out 

Program received signal SIGFPE: Floating-point exception - erroneous arithmetic operation.

Backtrace for this error:
#0  0x7fa4d0733d21 in ???
#1  0x7fa4d0732ef5 in ???
#2  0x7fa4d05640bf in ???
#3  0x5620d91911af in MAIN__
        at /mnt/f/programing/tools-code/fortran-exception-handle/main2.f90:6
#4  0x5620d9191260 in main
        at /mnt/f/programing/tools-code/fortran-exception-handle/main2.f90:8
[1]    1161 floating point exception  ./a.out