开启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