Fortran 程序失败取决于子例程调用之前的写入语句
Fortran program fails depending on a write statement before subroutine call
自从我使用 Fortran 以来已经有很多年了，所以也许我遗漏了一个基本问题，但就是这样。我什至不知道如何正确描述这个问题，所以我提前为缺乏描述性信息道歉。
我正在编写一些 Fortran 模块来补充使用 f2py 的 Python 程序。一切似乎都运行良好，但我在一个子程序中遇到了一些奇怪的错误。我无法在一个小示例程序中重现这个问题，所以我从模块中剥离了相关的子程序并生成了一个小的测试主程序。主程序为：
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
注意我所有的数组都是从0开始的。我这样做是因为我通常先使用numpy开发Python中的方法，然后用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
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
同样，数组从 0 开始，上限通常取决于子例程的输入。这个子程序调用了其他的，但是他们没有显示。
编译命令及结果如下图。你可以看到第一个结果是假的。使用 -fno-backtrace 的第二个结果是正确的结果。第三个结果和第一个一样编译，但是在调用子程序之前插入了一条write语句，结果是正确的。
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>
出于某种原因，在调用子例程之前添加一条写语句使其成为 "work." 我不完全熟悉 -fno-backtrace 选项，但它也使它成为 "work"。我在使用 f2py 编译时添加了这个选项，但我仍然得到奇怪的结果，但我猜是一次一件事。在Python中，我将循环调用此子程序 10 次，输入相同，10 次结果中有 8 次是正确的，但 2 次是假的，但我跑题了...
感谢您的帮助和任何建议。
更新 1：
子程序CURVE_DERIVS_ALG1如下所示。它也调用其他子例程，但为简洁起见未显示。我还使用 -fbounds-check 进行了编译，并得到了上面显示的相同的虚假结果。
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
更新 2：
很抱歉 post，但是整个模块都在下面 post 编辑，以防有人想尝试 运行 它与上面的主程序一起使用。
! 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
在子程序CURVE_DERIVS_ALG1中，伪参数CK好像没有初始化，请问可以检查一下它的初始值吗？如果我在进入循环之前将其设置为 0.0d0 ，代码似乎可以正常工作，但我不确定这个初始值是否可以。 （另请注意，如果给出 INTENT(OUT)，则必须定义所有元素。）
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, :)
...
另一个潜在的问题是
IF (U .EQ. UK(N + 1)) THEN
比较两个浮点数。虽然这个程序似乎不满足这个条件，但重写它可能更安全，例如
IF ( abs( U - UK(N + 1) ) < 1.0d-10 ) THEN !! threshold depends on your need...
编辑：要自动检测上述 CK 错误，编译为
可能会有用
gfortran -finit-real=snan -ffpe-trap=invalid evaluate.f90 main.f90
这给出（在 Linux x86_64 上使用 gfortran4.8）
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 以来已经有很多年了，所以也许我遗漏了一个基本问题，但就是这样。我什至不知道如何正确描述这个问题，所以我提前为缺乏描述性信息道歉。
我正在编写一些 Fortran 模块来补充使用 f2py 的 Python 程序。一切似乎都运行良好，但我在一个子程序中遇到了一些奇怪的错误。我无法在一个小示例程序中重现这个问题，所以我从模块中剥离了相关的子程序并生成了一个小的测试主程序。主程序为：
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
注意我所有的数组都是从0开始的。我这样做是因为我通常先使用numpy开发Python中的方法，然后用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
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
同样，数组从 0 开始，上限通常取决于子例程的输入。这个子程序调用了其他的，但是他们没有显示。
编译命令及结果如下图。你可以看到第一个结果是假的。使用 -fno-backtrace 的第二个结果是正确的结果。第三个结果和第一个一样编译，但是在调用子程序之前插入了一条write语句，结果是正确的。
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>
出于某种原因，在调用子例程之前添加一条写语句使其成为 "work." 我不完全熟悉 -fno-backtrace 选项，但它也使它成为 "work"。我在使用 f2py 编译时添加了这个选项，但我仍然得到奇怪的结果，但我猜是一次一件事。在Python中，我将循环调用此子程序 10 次，输入相同，10 次结果中有 8 次是正确的，但 2 次是假的，但我跑题了...
感谢您的帮助和任何建议。
更新 1：
子程序CURVE_DERIVS_ALG1如下所示。它也调用其他子例程，但为简洁起见未显示。我还使用 -fbounds-check 进行了编译，并得到了上面显示的相同的虚假结果。
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
更新 2： 很抱歉 post，但是整个模块都在下面 post 编辑，以防有人想尝试 运行 它与上面的主程序一起使用。
! 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
在子程序CURVE_DERIVS_ALG1中，伪参数CK好像没有初始化，请问可以检查一下它的初始值吗？如果我在进入循环之前将其设置为 0.0d0 ，代码似乎可以正常工作，但我不确定这个初始值是否可以。 （另请注意，如果给出 INTENT(OUT)，则必须定义所有元素。）
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, :)
...
另一个潜在的问题是
IF (U .EQ. UK(N + 1)) THEN
比较两个浮点数。虽然这个程序似乎不满足这个条件，但重写它可能更安全，例如
IF ( abs( U - UK(N + 1) ) < 1.0d-10 ) THEN !! threshold depends on your need...
编辑：要自动检测上述 CK 错误，编译为
gfortran -finit-real=snan -ffpe-trap=invalid evaluate.f90 main.f90
这给出（在 Linux x86_64 上使用 gfortran4.8）
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 ()