问题描述
我已经和Fortran合作了好几年了，所以也许我错过了一个根本性的问题，但是现在它已经结束了。我甚至不知道如何正确描述这个问题，所以我提前为缺乏描述性信息而道歉。我正在编写一些Fortran模块来补充Python程序使用f2py。似乎一切正常，但我在一个子程序中遇到了一些奇怪的错误。我无法在一个小样本程序中复制这个问题，所以我从模块中除去了相关的子程序并生成了一个小的测试主程序。主程序是：
$ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ b INTEGER :: N = 8，P = 2，D，I，J
双精度:: U，UK（0:11），CPW（0：8,0：3），CK（0：1） ，0：3）
D = 1
U = 0.45
UK =（/0.D0,0.D0,0.D0,0.25D0,0.25D0,0.5 D0,0.5D0,0.75D0，&
0.75D0,1.D0,1.D0,1.D0 /）
CPW（0，:) =（/1.D0 （1，:) =（/.707D0，.707D0,0.D0，.707D0 /）
CPW（2，：0.D0,0.D0,1.D0 /）
CPW ）=（/0.D0,1.D0,0.D0,1.D0 /）
CPW（3，:) =（/-.707D0,707D0,0.D0,707D0 /）
CPW（4，:) =（/-1.D0,0.D0,0.0D0.1D0 /）
CPW（5，:) =（/-.707D0,-.707D0 ，0.D0，.707D0 /）
CPW（6，:) =（/0.D0，-1.D0，0.D0,1.D0 /）
CPW（7，:) =（/.707D0，-707D0,0.D0，.707D0 /）
CPW（8，:) =（/1.D0,0.D0,0.D0,1.D0 /）
！第一个和第二个结果被注释掉了。
WRITE（*，*）FOO.BAR
CALL RAT_CURVE_DERIVS（N，P，UK，CPW，U，D，CK）
WRITE （*，'（100G15.5）'）（CK（I，J），J = 0,2）
END DO
END PROGRAM
请注意，我所有的数组都是从0.我这样做是因为我通常首先使用numpy开发方法，然后在Fortran中重写，而对于程序整体而言，将数组的值设置为0而不是1是更自然的。在实际的程序中，所有变量指定的主程序来自Python。
EVALUATE中的子程序RAT_CURVE_DERIVS是：
在子程序 CURVE_DERIVS_ALG1 ，伪参数 CK 似乎没有初始化，那么你可以检查它的初始值吗？如果我在进入循环之前将它设置为 0.0d0 ，代码似乎可以正常工作，但我不确定这个初始值是否正确。 （请注意，如果给出 INTENT（OUT），则必须定义所有元素。）
<$ p $ （0，D，CK）
...
DOUBLE PRECISION，INTENT（OUT）:: CK（0：D，P，UK，CPW，U，D，CK） ，0：3）
...
CALL FIND_SPAN（N，P，U，UK，SPAN）
CALL DERS_BASIS_FUNS（SPAN，U，P，DU，UK，M，NDERS）
CK（:, :) :) = 0.0d0 !! < --- here
DO K = 0，DU
DO J = 0，P
CK（K，:) = CK（K，:) + NDERS（ K，J）* CPW（SPAN - P + J，:)
...
另一个潜在的问题是：
$ b $ pre $ IF（U .EQ。UK（N + 1））THEN
比较两个浮点数。尽管这个条件在这个程序中似乎没有被满足，但重写它可能是比较安全的，例如
IF（abs（U- UK（N + 1））
编辑：要自动检测 CK 的上述错误，编译为
gfortran -finit-real = snan -ffpe-trap = invalid evaluate.f90 main.f90
（在Linux x86_64上使用gfortran4.8）
pre $ $ $ $ c $编程接收信号8（SIGFPE）：浮点异常。
Backtrace为这个错误：
＃0 0x00000039becac5f4在等待（）从/lib64/libc.so.6
＃1 0x00000039c501400d在?? （）从/usr/lib64/libgfortran.so.3
＃2 0x00000039c501582e在?? （）from /usr/lib64/libgfortran.so.3
＃3 0x00000039c50146ca in ?? （）从/usr/lib64/libgfortran.so.3
＃4<信号处理程序调用>
＃5 0x0000000000401787在__evaluate_MOD_curve_derivs_alg1（）< ---这里
＃在__evaluate_MOD_rat_curve_derivs 6 0x0000000000400fce（）
＃7 0x0000000000402b26在MAIN__（）
＃8 0x0000000000402cbb在main（）
It's been a number of years since I've worked with Fortran, so maybe I'm missing a fundamental issue, but here it goes. I'm not even sure how to properly describe this issue, so I apologize in advance for a lack of descriptive information.
I'm writing some Fortran modules to supplement a Python program using f2py. Everything seems to be working fine, but I am encountering some strange errors in one subroutine. I couldn't replicate the issue in a small sample program, so I stripped out the relevant subroutines from the module and generated a small test main program. The main program is:
PROGRAM MAIN
USE EVALUATE
IMPLICIT NONE
INTEGER :: N=8, P=2, D, I, J
DOUBLE PRECISION :: U, UK(0:11), CPW(0:8, 0:3), CK(0:1, 0:3)
D = 1
U = 0.45
UK = (/0.D0, 0.D0, 0.D0, 0.25D0, 0.25D0, 0.5D0, 0.5D0, 0.75D0, &
0.75D0, 1.D0, 1.D0, 1.D0 /)
CPW(0, :) = (/1.D0, 0.D0, 0.D0, 1.D0 /)
CPW(1, :) = (/.707D0, .707D0, 0.D0, .707D0 /)
CPW(2, :) = (/0.D0, 1.D0, 0.D0, 1.D0 /)
CPW(3, :) = (/-.707D0, .707D0, 0.D0, .707D0 /)
CPW(4, :) = (/-1.D0, 0.D0, 0.D0, 1.D0 /)
CPW(5, :) = (/-.707D0, -.707D0, 0.D0, .707D0 /)
CPW(6, :) = (/0.D0, -1.D0, 0.D0, 1.D0 /)
CPW(7, :) = (/.707D0, -.707D0, 0.D0, .707D0 /)
CPW(8, :) = (/1.D0, 0.D0, 0.D0, 1.D0 /)
! This is commented out for the first and second results.
WRITE(*,*) "FOO.BAR"
CALL RAT_CURVE_DERIVS(N, P, UK, CPW, U, D, CK)
WRITE(*,*) "WRITING RESULTS"
DO I = 0, D
WRITE(*, '(100G15.5)') (CK(I, J), J = 0, 2)
END DO
END PROGRAM
Note that all my arrays start at 0. I am doing this since I usually develop the methods in Python first using numpy and then rewrite in Fortran, and for the program as a whole, it's more natural to start arrays at 0 rather than 1. In the actual program all the variables specified in the main program come from Python.
The subroutine RAT_CURVE_DERIVS in EVALUATE is:
SUBROUTINE RAT_CURVE_DERIVS(N, P, UK, CPW, U, D, CK)
IMPLICIT NONE
!F2PY INTENT(IN) N, P, UK, CPW, U, D
!F2PY INTENT(OUT) CK
!F2PY DEPEND(N, P) UK
!F2PY DEPEND(N) CPW
!F2PY DEPEND(D) CK
INTEGER, INTENT(IN) :: N, P, D
DOUBLE PRECISION, INTENT(IN) :: U, UK(0:N + P + 1), CPW(0:N, 0:3)
DOUBLE PRECISION, INTENT(OUT) :: CK(0:D, 0:2)
INTEGER :: I, K, J, X
DOUBLE PRECISION :: BC, V(0:2), CDERS(0:D, 0:3)
DOUBLE PRECISION :: ADERS(0:D, 0:2), WDERS(0:D)
CALL CURVE_DERIVS_ALG1(N, P, UK, CPW, U, D, CDERS)
ADERS = CDERS(:, 0:2)
WDERS = CDERS(:, 3)
DO K = 0, D
V = ADERS(K, :)
DO I = 1, K
CALL BINOMIAL(K, I, BC)
V = V - BC * WDERS(I) * CK(K - I, :)
END DO
CK(K, :) = V / WDERS(0)
END DO
END SUBROUTINE RAT_CURVE_DERIVS
Again the arrays start at 0 and the upper bound usually depends on an input to the subroutine. This subroutine calls others, but they are not shown.
The compile commands and results are shown below. You can see the first results are bogus. The second results using -fno-backtrace are the correct results. The third results are compiled as the first, but a write statements is inserted before the call to the subroutine, and the results are correct.
C:\Users\Trevor\Documents\Temp>gfortran evaluate.f90 main.f90
C:\Users\Trevor\Documents\Temp>a.exe
WRITING RESULTS
-0.16453-170 0.19209E-33 0.69763E+58
0.70809E-65 -0.82668E+72 -Infinity
C:\Users\Trevor\Documents\Temp>gfortran evaluate.f90 main.f90 -fno-backtrace
C:\Users\Trevor\Documents\Temp>a.exe
WRITING RESULTS
-0.95586 0.29379 0.0000
-1.8340 -5.9662 0.0000
C:\Users\Trevor\Documents\Temp>gfortran evaluate.f90 main.f90
C:\Users\Trevor\Documents\Temp>a.exe
FOO.BAR
WRITING RESULTS
-0.95586 0.29379 0.0000
-1.8340 -5.9662 0.0000
C:\Users\Trevor\Documents\Temp>
For some reason, adding a write statement before calling the subroutine makes it "work." I am not completely familiar with the -fno-backtrace option, but it makes it "work" too. I added this option when compiling using f2py, and I still get strange results, but one thing at a time I guess. In Python, I will call this subroutine 10 times in a loop with the same inputs, and 8 out of 10 result will be correct, but 2 will be bogus, but I digress...
Thanks for the help and any suggestions.
UPDATE 1:
The subroutine CURVE_DERIVS_ALG1 is shown below. It too calls other subroutines, but they are not shown for brevity. I also compiled with -fbounds-check and got the same bogus results shown above.
SUBROUTINE CURVE_DERIVS_ALG1(N, P, UK, CPW, U, D, CK)
IMPLICIT NONE
!F2PY INTENT(IN) N, P, UK, CPW, U, D
!F2PY INTENT(OUT) CK
!F2PY DEPEND(N, P) UK
!F2PY DEPEND(N) CPW
!F2PY DEPEND(D) CK
INTEGER, INTENT(IN) :: N, P, D
DOUBLE PRECISION, INTENT(IN) :: U, UK(0:N + P + 1), CPW(0:N, 0:3)
DOUBLE PRECISION, INTENT(OUT) :: CK(0:D, 0:3)
INTEGER :: DU, K, SPAN, J, M
DOUBLE PRECISION :: NDERS(0:MIN(D,P), 0:P)
DU = MIN(D, P)
M = N + P + 1
CALL FIND_SPAN(N, P, U, UK, SPAN)
CALL DERS_BASIS_FUNS(SPAN, U, P, DU, UK, M, NDERS)
DO K = 0, DU
DO J = 0, P
CK(K, :) = CK(K, :) + NDERS(K, J) * CPW(SPAN - P + J, :)
END DO
END DO
END SUBROUTINE CURVE_DERIVS_ALG1
UPDATE 2:Sorry for the long post, but the entire module is posted below in case anyone wants to try and run it with the main program above.
! FORTRAN source for geometry.tools.evaluate
MODULE EVALUATE
CONTAINS
SUBROUTINE FIND_SPAN(N, P, U, UK, MID)
IMPLICIT NONE
!F2PY INENT(IN) N, P, U, UK
!F2PY INTENT(OUT) MID
!F2PY DEPEND(N, P) UK
INTEGER, INTENT(IN) :: N, P
DOUBLE PRECISION, INTENT(IN) :: U
DOUBLE PRECISION, INTENT(IN) :: UK(0:N + P + 1)
INTEGER, INTENT(OUT) :: MID
INTEGER :: LOW, HIGH
! SPECIAL CASE
IF (U .EQ. UK(N + 1)) THEN
MID = N
RETURN
END IF
LOW = P
HIGH = N + 1
MID = (LOW + HIGH) / 2
DO WHILE ((U .LT. UK(MID)) .OR. (U .GE. UK(MID + 1)))
IF (U .LT. UK(MID)) THEN
HIGH = MID
ELSE
LOW = MID
END IF
MID = (LOW + HIGH) / 2
END DO
END SUBROUTINE FIND_SPAN
SUBROUTINE BASIS_FUNS(I, U, P, UK, M, N)
IMPLICIT NONE
!F2PY INTENT(IN) I, U, P, UK, M
!F2PY INTENT(OUT) N
!F2PY DEPEND(M) UK
INTEGER, INTENT(IN) :: I, P, M
DOUBLE PRECISION, INTENT(IN) :: U
DOUBLE PRECISION, INTENT(IN) :: UK(0:M)
DOUBLE PRECISION, INTENT(OUT) :: N(0:P)
INTEGER :: J, R
DOUBLE PRECISION :: TEMP, SAVED
DOUBLE PRECISION :: LEFT(0:P), RIGHT(0:P)
N(0) = 1.D0
DO J = 1, P
LEFT(J) = U - UK(I + 1 - J)
RIGHT(J) = UK(I + J) - U
SAVED = 0.D0
DO R = 0, J - 1
TEMP = N(R) / (RIGHT(R + 1) + LEFT(J - R))
N(R) = SAVED + RIGHT(R + 1) * TEMP
SAVED = LEFT(J - R) * TEMP
END DO
N(J) = SAVED
END DO
END SUBROUTINE BASIS_FUNS
SUBROUTINE DERS_BASIS_FUNS(I, U, P, N, UK, M, DERS)
IMPLICIT NONE
!F2PY INTENT(IN) I, U, P, N, UK, M
!F2PY INTENT(OUT) DERS
!F2PY DEPEND(M) UK
INTEGER, INTENT(IN) :: I, P, N, M
DOUBLE PRECISION, INTENT(IN) :: U
DOUBLE PRECISION, INTENT(IN) :: UK(0:M)
DOUBLE PRECISION, INTENT(OUT) :: DERS(0:N, 0:P)
INTEGER :: J, K, R, J1, J2, RK, PK, S1, S2
DOUBLE PRECISION :: SAVED, TEMP, NDU(0:P, 0:P), LEFT(0:P), &
RIGHT(0:P), A(0:1, 0:P), D
NDU(0, 0) = 1.D0
DO J = 1, P
LEFT(J) = U - UK(I + 1 - J)
RIGHT(J) = UK(I + J) - U
SAVED = 0.D0
DO R = 0, J - 1
NDU(J, R) = RIGHT(R + 1) + LEFT(J - R)
TEMP = NDU(R, J - 1) / NDU(J, R)
NDU(R, J) = SAVED + RIGHT(R + 1) * TEMP
SAVED = LEFT(J - R) * TEMP
END DO
NDU(J, J) = SAVED
END DO
DO J = 0, P
DERS(0, J) = NDU(J, P)
END DO
DO R = 0, P
S1 = 0
S2 = 1
A(0, 0) = 1.D0
DO K = 1, N
D = 0.D0
RK = R - K
PK = P - K
IF (R .GE. K) THEN
A(S2, 0) = A(S1, 0) / NDU(PK + 1, RK)
D = A(S2, 0) * NDU(RK, PK)
END IF
IF (RK .GE. -1) THEN
J1 = 1
ELSE
J1 = -RK
END IF
IF (R - 1 .LE. PK) THEN
J2 = K - 1
ELSE
J2 = P - R
END IF
DO J = J1, J2
A(S2, J) = (A(S1, J) - A(S1, J - 1)) / &
NDU(PK + 1, RK + J)
D = D + A(S2, J) * NDU(RK + J, PK)
END DO
IF (R .LE. PK) THEN
A(S2, K) = -A(S1, K - 1) / NDU(PK + 1, R)
D = D + A(S2, K) * NDU(R, PK)
END IF
DERS(K, R) = D
J = S1
S1 = S2
S2 = J
END DO
END DO
R = P
DO K = 1, N
DO J = 0, P
DERS(K, J) = DERS(K, J) * R
END DO
R = R * (P - K)
END DO
END SUBROUTINE DERS_BASIS_FUNS
SUBROUTINE CURVE_DERIVS_ALG1(N, P, UK, CPW, U, D, CK)
IMPLICIT NONE
!F2PY INTENT(IN) N, P, UK, CPW, U, D
!F2PY INTENT(OUT) CK
!F2PY DEPEND(N, P) UK
!F2PY DEPEND(N) CPW
!F2PY DEPEND(D) CK
INTEGER, INTENT(IN) :: N, P, D
DOUBLE PRECISION, INTENT(IN) :: U, UK(0:N + P + 1), CPW(0:N, 0:3)
DOUBLE PRECISION, INTENT(OUT) :: CK(0:D, 0:3)
INTEGER :: DU, K, SPAN, J, M
DOUBLE PRECISION :: NDERS(0:MIN(D,P), 0:P)
DU = MIN(D, P)
M = N + P + 1
CALL FIND_SPAN(N, P, U, UK, SPAN)
CALL DERS_BASIS_FUNS(SPAN, U, P, DU, UK, M, NDERS)
DO K = 0, DU
DO J = 0, P
CK(K, :) = CK(K, :) + NDERS(K, J) * CPW(SPAN - P + J, :)
END DO
END DO
END SUBROUTINE CURVE_DERIVS_ALG1
SUBROUTINE RAT_CURVE_DERIVS(N, P, UK, CPW, U, D, CK)
IMPLICIT NONE
!F2PY INTENT(IN) N, P, UK, CPW, U, D
!F2PY INTENT(OUT) CK
!F2PY DEPEND(N, P) UK
!F2PY DEPEND(N) CPW
!F2PY DEPEND(D) CK
INTEGER, INTENT(IN) :: N, P, D
DOUBLE PRECISION, INTENT(IN) :: U, UK(0:N + P + 1), CPW(0:N, 0:3)
DOUBLE PRECISION, INTENT(OUT) :: CK(0:D, 0:2)
INTEGER :: I, K, J, X
DOUBLE PRECISION :: BC, V(0:2), CDERS(0:D, 0:3)
DOUBLE PRECISION :: ADERS(0:D, 0:2), WDERS(0:D)
CALL CURVE_DERIVS_ALG1(N, P, UK, CPW, U, D, CDERS)
ADERS = CDERS(:, 0:2)
WDERS = CDERS(:, 3)
DO K = 0, D
V = ADERS(K, :)
DO I = 1, K
CALL BINOMIAL(K, I, BC)
V = V - BC * WDERS(I) * CK(K - I, :)
END DO
CK(K, :) = V / WDERS(0)
END DO
END SUBROUTINE RAT_CURVE_DERIVS
SUBROUTINE BINOMIAL(N, K, BC)
IMPLICIT NONE
!F2PY INTENT(IN) N, K
!F2PY INTENT(OUT) BC
INTEGER, INTENT(IN) :: N, K
DOUBLE PRECISION, INTENT(OUT) :: BC
INTEGER :: I, KK
IF ((K .LT. 0) .OR. ( K .GT. N)) THEN
BC = 0.D0
RETURN
END IF
IF ((K .EQ. 0) .OR. ( K .EQ. N)) THEN
BC = 1.D0
RETURN
END IF
KK = MIN(K, N - K)
BC = 1.D0
DO I = 0, KK - 1
BC = BC * DBLE(N - I) / DBLE(I + 1)
END DO
END SUBROUTINE BINOMIAL
END MODULE
In the subroutine CURVE_DERIVS_ALG1, the dummy argument CK seems not initialized, so could you check its initial value? If I set it to 0.0d0 before entering the loop, the code seems to work fine, but I am not sure if this initial value is OK. (Please also note that if INTENT(OUT) is given, all the elements must be defined.)
SUBROUTINE CURVE_DERIVS_ALG1(N, P, UK, CPW, U, D, CK)
...
DOUBLE PRECISION, INTENT(OUT) :: CK(0:D, 0:3)
...
CALL FIND_SPAN(N, P, U, UK, SPAN)
CALL DERS_BASIS_FUNS(SPAN, U, P, DU, UK, M, NDERS)
CK(:,:) = 0.0d0 !! <--- here
DO K = 0, DU
DO J = 0, P
CK(K, :) = CK(K, :) + NDERS(K, J) * CPW(SPAN - P + J, :)
...
Another potential issue is
IF (U .EQ. UK(N + 1)) THEN
which compares two floating-point numbers. Although this condition seems not met in this program, it is probably safer to rewrite this as, e.g.
IF ( abs( U - UK(N + 1) ) < 1.0d-10 ) THEN !! threshold depends on your need...
EDIT: To detect the above error of CK automatically, it may be useful to compile as
gfortran -finit-real=snan -ffpe-trap=invalid evaluate.f90 main.f90
which gives (with gfortran4.8 on Linux x86_64)
Program received signal 8 (SIGFPE): Floating-point exception. Backtrace for this error: #0 0x00000039becac5f4 in wait () from /lib64/libc.so.6 #1 0x00000039c501400d in ?? () from /usr/lib64/libgfortran.so.3 #2 0x00000039c501582e in ?? () from /usr/lib64/libgfortran.so.3 #3 0x00000039c50146ca in ?? () from /usr/lib64/libgfortran.so.3 #4 <signal handler called> #5 0x0000000000401787 in __evaluate_MOD_curve_derivs_alg1 () <--- here #6 0x0000000000400fce in __evaluate_MOD_rat_curve_derivs () #7 0x0000000000402b26 in MAIN__ () #8 0x0000000000402cbb in main ()
这篇关于子程序调用之前，Fortran程序根据写入语句而失败的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持！