2016-09-21 5 views
1

b- 스플라인 및 재귀 함수를 사용하여 호일 섹션을 만들려고합니다. 문제는 매우 큰 부동 소수점 수 또는 NaN 및 심지어 무한대의 값 중 일부를 반환한다는 것입니다. 어떤 값들은 다시 컴파일 한 후에 변경되는 것으로 보이며 때로는 어떤 이유로 인해 다른 것보다 더 많이 변경되는 것으로 보입니다. 정상적으로 반환되는 값은 정확합니다.Fortran에서 재귀 함수를 사용할 때 임의의 큰 부동 소수점 오류 및 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 
+1

어딘가에서 변수를 초기화하는 것을 잊어 버린 것 같아요. 아마도 여기 누군가가 당신을 위해 그것을 찾을 수 있지만 도움이 될만한 도구가 많이 있습니다. 시작하려면 컴파일러 런타임 검사를 사용 했습니까? (ifort check-uninit) 예 : – Ross

+0

모듈 당 하나의 '암시 적 없음'으로 충분하다는 것을 알고 계셨습니까? 또한 유닛 번호 '4'를 사용하는 것은 좋지 않습니다. 큰 숫자를 사용하십시오. –

+0

들여 쓰기를 조금 더 일관되게하려고했습니다. 그것은 여전히 ​​개선 될 수 있지만 너무 많은 작업이 될 것입니다. 나는 컴파일러에서 찾을 수있는 모든 디버깅 플래그로 컴파일 할 제안을한다. ('gfortran -g -fbacktrace -fcheck = all -Wall' 또는'gfotran -g -fsanitize = address, undefined'). –

답변

1

죄송합니다 위의 유용한 주석 스레드에서 벗어날 수 있습니다. 문제의 두 가지 핵심은 다음과 같습니다. 코드를 gfortran으로 다시 컴파일했습니다 (Codeblocks를 사용하면 중요하지 않음). 사용 가능한 컴파일러 플래그를 활성화해도 문제를 정확하게 지적하지 않는다는 것에 동의합니다. 처음에는 바운드 검사 문제 (배열의 바깥에있는 초기화되지 않은 메모리 블록을 가리키는 것)가있는 것처럼 보였습니다.하지만 그렇지 않았습니다.

그러나 - 나는 컴파일러에 대해 다음 플래그를 설정 :

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

을 Codeblocks에서 당신이 Project>Build Options>Other Compiler Options에서이 작업을 수행하고 (나는 그것이 다른 IDE에는 변함이 없을 것입니다 상상) 수동으로 추가 할 수 있습니다. 그것이하는 것은 오류가있는 부동 소수점 연산 (0으로 나누기, 변수 유형으로 지정된 값 범위 밖의 숫자 등)이 발생하면 오류가 발생한다는 것입니다. 이러한 문제로 인해 프로그램에서 여러 가지 문제가 발생할 수 있습니다 (예 : denormal은 계산 과정을 중단하여 CPU가 비정규 숫자를 반올림/나타낼 때까지 계산 함) 잘못된 결과를 완전히 줄 수 있습니다 (zero). NaN 에스).

디버깅 도구를 사용하면 이제 프로그램을 통해 디버깅 도구를 찾아보고 수정할 수 있습니다. 이 플래그와 프로그램의 수치 오류를 찾아내는 다른 전략에 대한 자세한 내용은 this GFortran manual section을 확인하십시오. 행운을 빌어 요!