问题描述
我正在尝试学习 Fortran(不幸的是，这对我的研究小组来说是必需品) - 我为自己设定的任务之一是将 Numerical Recipes 书中的一个必要函数(关联勒让德多项式)打包成符合 fortran 03 的模块.原程序(f77)有如下形式的一些错误处理:
I'm trying to learn Fortran (unfortunately a necessity for my research group) - one of the tasks I set myself was to package one of the necessary functions (Associated Legendre polynomials) from the Numerical Recipes book into a fortran 03 compliant module. The original program (f77) has some error handling in the form of the following:
if(m.lt.0.or.m.gt.1.or.abs(x).gt.1)pause 'bad arguments in plgndr'
自 f77 以来，Pause 似乎已被弃用，因为使用此行会给我一个编译错误，所以我尝试了以下操作:
Pause seems to have been deprecated since f77 as using this line gives me a compiling error, so I tried the following:
module sha_helper
implicit none
public :: plgndr, factorial!, ylm
contains
! numerical recipes Associated Legendre Polynomials rewritten for f03
function plgndr(l,m,x) result(res_plgndr)
integer, intent(in) :: l, m
real, intent(in) :: x
real :: res_plgndr, fact, pll, pmm, pmmp1, somx2
integer :: i,ll
if (m.lt.0.or.m.gt.l.or.abs(x).gt.1) then
write (*, *) "bad arguments to plgndr, aborting", m, x
res_plgndr=-10e6 !return a ridiculous value
else
pmm = 1.
if (m.gt.0) then
somx2 = sqrt((1.-x)*(1.+x))
fact = 1.
do i = 1, m
pmm = -pmm*fact*somx2
fact = fact+2
end do
end if
if (l.eq.m) then
res_plgndr = pmm
else
pmmp1 = x*(2*m+1)*pmm
if(l.eq.m+1) then
res_plgndr = pmmp1
else
do ll = m+2, l
pll = (x*(2*ll-1)*pmmp1-(ll+m-1)*pmm)/(ll-m)
pmm = pmmp1
pmmp1 = pll
end do
res_plgndr = pll
end if
end if
end if
end function plgndr
recursive function factorial(n) result(factorial_result)
integer, intent(in) :: n
integer, parameter :: RegInt_K = selected_int_kind(20) !should be enough for the factorials I am using
integer (kind = RegInt_K) :: factorial_result
if (n <= 0) then
factorial_result = 1
else
factorial_result = n * factorial(n-1)
end if
end function factorial
! function ylm(l,m,theta,phi) result(res_ylm)
! integer, intent(in) :: l, m
! real, intent(in) :: theta, phi
! real :: res_ylm, front_block
! real, parameter :: pi = 3.1415926536
! front_block = sqrt((2*l+1)*factorial(l-abs(m))/(4*pi*))
! end function ylm
end module sha_helper
else 之后的主代码可以工作，但是如果我执行我的主程序并使用错误值调用该函数，则程序会在执行打印语句之前冻结.我知道 print 语句是问题所在，因为注释掉它可以让函数正常执行，返回 -10e6 作为值.理想情况下，我希望程序在给出用户可读的错误消息后崩溃，因为给 plgndr 函数提供错误的值是程序的致命错误.程序 sha_lmc 正在使用函数 plgndr.目前所做的只是读取一些数组，然后打印一个 plgndr 的值进行测试(早期).模块 sha_helper 中的函数 ylm 也没有完成，因此被注释掉了.该代码使用 gfortran sha_helper.f03 sha_lmc.f03 -o sha_lmc 编译，并且gfortran --版本GNU Fortran (GCC) 4.8.2
The main code after the else works, but if I execute my main program and call the function with bad values, the program freezes before executing the print statement. I know that the print statement is the problem, as commenting it out allows the function to execute normally, returning -10e6 as the value. Ideally, I would like the program to crash after giving a user readable error message, as giving bad values to the plgndr function is a fatal error for the program. The function plgndr is being used by the program sha_lmc. Currently all this does is read some arrays and then print a value of plgndr for testing (early days). The function ylm in the module sha_helper is also not finished, hence it is commented out. The code compiles using gfortran sha_helper.f03 sha_lmc.f03 -o sha_lmc, andgfortran --versionGNU Fortran (GCC) 4.8.2
!Spherical Harmonic Bayesian Analysis testbed for Lagrangian Dynamical Monte Carlo
program sha_analysis
use sha_helper
implicit none
!Analysis Parameters
integer, parameter :: harm_order = 6
integer, parameter :: harm_array_length = (harm_order+1)**2
real, parameter :: coeff_lo = -0.1, coeff_hi = 0.1, data_err = 0.01 !for now, data_err fixed rather than heirarchical
!Monte Carlo Parameters
integer, parameter :: run = 100000, burn = 50000, thin = 100
real, parameter :: L = 1.0, e = 1.0
!Variables needed by the program
integer :: points, r, h, p, counter = 1
real, dimension(:), allocatable :: x, y, z
real, dimension(harm_array_length) :: l_index_list, m_index_list
real, dimension(:,:), allocatable :: g_matrix
!Open the file, allocate the x,y,z arrays and read the file
open(1, file = 'Average_H_M_C_PcP_boschi_1200.xyz', status = 'old')
read(1,*) points
allocate(x(points))
allocate(y(points))
allocate(z(points))
print *, "Number of Points: ", points
readloop: do r = 1, points
read(1,*) x(r), y(r), z(r)
end do readloop
!Set up the forwards model
allocate(g_matrix(harm_array_length,points))
!Generate the l and m values of spherical harmonics
hloop: do h = 0, harm_order
ploop: do p = -h,h
l_index_list(counter) = h
m_index_list(counter) = p
counter = counter + 1
end do ploop
end do hloop
print *, plgndr(1,2,0.1)
!print *, ylm(1,1,0.1,0.1)
end program sha_analysis
推荐答案
您的程序执行所谓的递归 IO - 对 plgndr 的初始调用位于 IO 语句的输出项列表中(打印语句)[将输出定向到控制台] - 在该函数内部，您还尝试执行另一个 IO 语句[输出到控制台].这是不允许的 - 请参阅 F2003 的 9.11p2 和 p3 或 F2008 的 9.12p2.
Your program does what is known as recursive IO - the initial call to plgndr is in the output item list of an IO statement (a print statement) [directing output to the console] - inside that function you then also attempt to execute another IO statement [that outputs to the console]. This is not permitted - see 9.11p2 and p3 of F2003 or 9.12p2 of F2008.
解决办法是把函数调用和主程序中的io语句分开，即
A solution is to separate the function invocation from the io statement in the main program, i.e.
REAL :: a_temporary
...
a_temporary = plgndr(1,2,0.1)
PRINT *, a_temporary
F2008 中的其他替代方案(但不是 F2003 - 因此第一段中的 [ ] 部分)包括将函数的输出定向到不同的逻辑单元(注意 WRITE (*, ... 和 PRINT ... 引用相同的单位).
Other alternatives in F2008 (but not F2003 - hence the [ ] parts in the first paragraph) include directing the output from the function to a different logical unit (note that WRITE (*, ... and PRINT ... reference the same unit).
在 F2008 中，您还可以将 WRITE 语句替换为带有消息的 STOP 语句(该消息必须是常量 - 这不会让您报告有问题的值).
In F2008 you could also replace the WRITE statement with a STOP statement with a message (the message must be a constant - which wouldn't let you report the problematic values).
无意中调用递归 IO 的可能性是某些编程风格不鼓励在函数中执行 IO 的部分原因.
The potential for inadvertently invoking recursive IO is part of the reason that some programming styles discourage conducting IO in functions.
这篇关于从 Fortran 模块中定义的函数打印到标准输出的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持！