子程序调用之前，Fortran程序根据写入语句而失败 [英] Fortran program fails depending on a write statement before subroutine call
问题描述
我正在编写一些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是：
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
双精度，INTENT（IN）:: U，UK（0：N + P + 1），CP W（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-1，:)
END DO
CK（K，:) = V / WDERS（0）
END DO
END SUBROUTINE RAT_CURVE_DERIVS
阵列从0开始，上限通常取决于子例程的输入。这个子程序调用其他程序，但不显示。
编译命令和结果如下所示。你可以看到第一个结果是假的。使用-fno-backtrace的第二个结果是正确的结果。第三个结果编译为第一个结果，但在调用子例程之前插入了写入语句，结果正确。
C：\ Users \Trevor\Documents\Temp> gfortran evaluate.f90 main.f90
C：\ Users \Trevor\Documents\Temp> a.exe
书写结果
-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-回溯
C：\Users\Trevor\Documents\Temp> A.EXE
结果写入
-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
书写结果
-0.95586 0.29379 0.0000
-1.8340 -5.9662 0.0000
C：\ Users \Trevor\Documents\Temp>
出于某种原因，在调用子例程之前添加写入语句会使其工作。我不完全熟悉-fno-backtrace选项，但它也使其工作。我在使用f2py编译时添加了这个选项，但我仍然得到了奇怪的结果，但是我一次只能看到一个。在Python中，我会用相同的输入在循环中调用这个子程序10次，并且10个结果中的8个将是正确的，但是2将是假的，但我离题了......
感谢您的帮助和任何建议。
更新1：
子程序CURVE_DERIVS_ALG1是如下所示。它也调用其他子程序，但为了简洁起见，它们未被显示。
SUBROUTINE CURVE_DERIVS_ALG1（N，P，UK， CPW，U，D，CK）
隐含无
！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
双精度，INTENT（IN）:: U，UK（0：N + P + 1），CPW（0：N，0：3）
双精度，意图（OUT）:: CK（0：D，0：3）
INTEGER :: DU，K，SPAN，J，M
双精度:: 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
更新2：
对不起，对于很长的文章，但整个模块发布在下面，以防任何人想要用上面的主程序运行它。
！ FORTRAN源代码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
！特殊情况
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
双精度，目标（上）:: U
双精度，目的（英）:: UK（0 ：M）
DOUBLE PRECISION，INTENT（OUT）:: N（0：P）
INTEGER :: J，R
双精度:: 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）=已存储
END DO
END SUBROUTINE BASIS_FUNS
SUBROUTINE DERS_BASIS_FUNS（I，U，P，N，UK，M，DERS）
隐式无
！F2PY INTENT（IN）I，U，P ，N，UK，M
！F2PY INTENT（OUT）DERS
！F2PY DEPEND（M）英国
整数，意图（上）:: 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
双精度:: 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 ，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（RGE.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 .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）
隐式无
！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（ （0：N + P + 1），CPW（0：N，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）
隐式无无
！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）
INT EGER :: I，K，J，X
DOUBLE PRECISION :: BC，V（0：2），CDERS（0：D，0：3）
双精度:: 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-1，:)
END DO
CK WDERS（0）
END DO
END SUBROUTINE RAT_CURVE_DERIVS
SUBROUTINE BINOMIAL（N，K，BC）
隐式无
！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-1）/ DBLE（I + 1）
END DO
END SUBROUTINE BINOMIAL
END MODULE
在子程序 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））<1.0d-10）THEN !!阈值取决于您的需要...
编辑：要自动检测 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程序根据写入语句而失败的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
