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と 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

上記の有益なコメントスレッドから逸脱する申し訳ありません。ここで私の2セントは問題です - あなたのコードをgfortranで再コンパイルしました（Codeblocksを使用して、重要ではありません）。使用可能なコンパイラフラグを有効にしても問題が特定できないと私は同意します。当初は、バウンドチェックの問題（配列の境界の外にあるメモリブロックを指していますが、初期化されていません）があるように見えましたが、そうではありませんでした。

は、しかし - 私は、コンパイラのために、次のフラグをオン：

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

をコードブロックでは、あなたがProject>Build Options>Other Compiler Optionsでこれを行うと、（私はそれが他のIDEで同様であろうと想像）これを手動で追加することができます。それは、誤った浮動小数点演算（ゼロ除算、変数の型で指定された値の範囲外の数など）が発生したときにエラーを発生させることです。このような問題は、あなたのプログラムで何か問題を引き起こす可能性があります。例えば、数値計算の処理を中断してデノーマル・ナンバーを丸めたり表現するまでの間、パフォーマンスが大幅に低下することがあります（zeroNaN s）。

デバッグツールを使用すると、プログラム内でデバッグツールを使用して修正できます。プログラムの数値エラーを捜すためのこのフラグと他の戦略の詳細については、this GFortran manual sectionを参照してください。幸運の狩りミス！