在Fortran中使用递归函数时出现随机大浮点错误和NaN

Question

我正在尝试创建一个使用b-样条线和递归函数的铝箔截面。问题是它会返回一些非常大的浮点数或NaN格式的值，有时甚至是无穷大的。哪个值在重新编译后似乎会发生变化，有时由于某些原因会比其他值更多。它正常返回的值是正确的。

有没有人有任何建议问题可能是什么？任何帮助都是非常感谢的。我知道我这里的代码很长，但我认为一切都是必要的。首先是程序，然后是带有子例程和函数的模块。

正确答案应该是：xb = 0 0.0147 0.0529 0.1063 0.1675 0.2448

 `0.3564    0.5202    0.7231    0.9047    1.0000` and `yb = 0   -0.0088   -0.0166   -0.0227   -0.0264   -0.0272 -0.0245   -0.0178   -0.0084   -0.0009         0`

module splines
     IMPLICIT NONE

     interface CoxdeBoor
        module procedure CoxdeBoor
     end interface CoxdeBoor


CONTAINS    

SUBROUTINE Splinefoil(N, degree, params, pts, xb, yb)
    INTEGER :: degree, N, i, pts
    REAL :: Pxb(pts+2), Pyb(pts+2)
    REAL :: tegap
    REAL, dimension (:):: yb, xb, params
    tegap=0.00; 

    Pxb(1:2)=0
    Pyb(1)=0
    DO i=1, pts-1
        !x-coordinates of bottom
        Pxb(i+2)=params(((2*i)-1)+2*pts)
        !y-coordinates of bottom
        Pyb(i+1)=params(2*i-2+2*pts)
    END DO
    Pyb(pts+1)=params(4*pts-2)
    Pxb(pts+2)=1
    Pyb(pts+2)=-0.5*tegap

    call Bspline(N,degree,Pxb, Pyb, pts,xb,yb)
END SUBROUTINE Splinefoil 

SUBROUTINE Bspline(N, degree, Px, Py, pts, x, y)
    INTEGER :: degree, N, j, i, n1, pts, l
    INTEGER, parameter :: len=10
    REAL :: B(pts+2), Qx(N/2+1), Qy(N/2+1), t(N/2+1), dt
    REAL, dimension(1:len) :: T2
    REAL, dimension (:):: Px, Py, y, x 
    n1=pts+1

    DO i=1, degree+n1+2
        IF (i<(degree+2)) THEN
            T2(i)=0
        ELSE IF (((degree+2) <= i) .and. (i < (n1 +2))) THEN
            T2(i)=(i-degree-1)/real((n1+1-degree))
        ELSE
            T2(i)=1
        END IF
    END DO

    t(1)=0
    dt=1/real((N/2))
    DO j=1, N/2
        t(j+1)=t(j) + dt 
    END DO
    Qx(1)=0
    DO l=1, N/2+1
        DO i=0, n1
            B(i+1)= CoxdeBoor(degree,i+1, T2,t(l))
            x(l)=x(l) + B(i+1) * Px(i+1)
        END DO
    END DO

    DO l=1, N/2+1
        DO i=0, n1
            B(i+1)= CoxdeBoor(degree,i+1, T2,t(l))
            y(l)=y(l) + B(i+1) * Py(i+1)
        END DO
    END DO  
END SUBROUTINE Bspline

RECURSIVE FUNCTION CoxdeBoor(i, j, x, t) result(val)
    INTEGER :: i, j, m
    REAL :: val, t, t1, t2
    REAL, dimension(:) :: x
    m=10
    IF (i==0) THEN
        IF (x(j)<=t .and. t<x(j+1)) THEN
            val=1
        ELSE IF (x(j)<=t .and. t==x(j+1) .and. x(j+1)==1) THEN
            val=1
        ELSE
            val=0
        END IF
        ELSE
        IF (x(j)<x(j+i)) THEN
            t1 =CoxdeBoor(i-1,j,x,t)
            val=(t - x(j)) / (x(j+i) - x(j)) * t1
        ELSE
            val=0
         END IF
        IF (j<m) THEN
            IF (x(j+1)<x(j+i+1)) THEN
                t2 =CoxdeBoor(i-1,j+1,x,t)
                 val=val + (x(j+i+1) - t) / (x(j+i+1) - x(j+1)) * t2
            END IF
        END IF
     END IF
END FUNCTION CoxdeBoor

END MODULE splines

Program TEST
     USE splines
     implicit none
     INTEGER, parameter :: N=20, pts=4
     INTEGER :: degree, i
     REAL :: params(14), xb(N/2+1), yb(N/2+1), xt(N/2+1), yt(N/2+1), x(N+1), y(N+1)

     params=(/0.010021287940814, 0.069234038308141, 0.039810312675194, 0.602154240414724, 0.027370571639440, 0.705186051614965,&
         0.015116139770247, -0.010144644178286, 0.119067997228366, -0.028919194962962, 0.338791094291084, -0.028735107857216,&
         0.965914604459008, 0.004397962157839/)

     degree =3

     call Splinefoil(N, degree, params, pts, xt, yt)

     write(*,*) xt
     write(*,*) yt
     OPEN(4, FILE='foil3.dat')
     WRITE(4, '(2F10.4)') (xt(i), yt(i), i=1, 11)
     CLOSE(4)
END PROGRAM TEST

Answer 1

很抱歉偏离了上面有用的评论帖子。这是我在这个问题上的两点看法--我用gfortran重新编译了你的代码(使用代码块，这并不重要)。我同意启用可用的编译器标志无助于查明问题。最初，它看起来像是存在边界检查问题(您引用的是数组边界之外的内存块，该块未初始化)，但事实并非如此。

但是-我为编译器打开了以下标志：

-ffpe-trap=invalid,zero,overflow,underflow,inexact,denormal

在代码块中，您可以在Project>Build Options>Other Compiler Options中完成此操作，并手动添加此代码(我想在其他IDE中也是如此)。它所做的是在遇到错误的浮点操作(除以零、变量类型指定的值范围之外的数字等)时引发错误。这样的问题可能会在你的程序中导致许多问题--从显著降低性能(例如，denormal会中断计算过程，直到CPU找出如何舍入/表示非正规化数字)到给你完全错误的结果(zero给你NaNs)。

使用调试工具，您现在可以在程序中检查并修复这些问题。有关此标志的更多信息以及其他用于查找程序中的数字错误的策略，请查看this GFortran manual section。祝你好运，错误搜索！