pr0*_*amR 0 floating-point gradient double-precision fortran95
我在 Fortran 95 中编写了一个基本算法,使用通过称为理查森外推法的过程增强的中心差来计算函数的梯度(代码中规定了一个示例)。
function f(n,x)
! The scalar multivariable function to be differentiated
integer :: n
real(kind = kind(1d0)) :: x(n), f
f = x(1)**5.d0 + cos(x(2)) + log(x(3)) - sqrt(x(4))
end function f
!=====!
!=====!
!=====!
program gradient
!==============================================================================!
! Calculates the gradient of the scalar function f at x=0using a finite !
! difference approximation, with a low order Richardson extrapolation. !
!==============================================================================!
parameter (n = 4, M = 25)
real(kind = kind(1d0)) :: x(n), xhup(n), xhdown(n), d(M), r(M), dfdxi, h0, h, gradf(n)
h0 = 1.d0
x = 3.d0
! Loop through each component of the vector x and calculate the appropriate
! derivative
do i = 1,n
! Reset step size
h = h0
! Carry out M successive central difference approximations of the derivative
do j = 1,M
xhup = x
xhdown = x
xhup(i) = xhup(i) + h
xhdown(i) = xhdown(i) - h
d(j) = ( f(n,xhup) - f(n,xhdown) ) / (2.d0*h)
h = h / 2.d0
end do
r = 0.d0
do k = 3,M r(k) = ( 64.d0*d(k) - 20.d0*d(k-1) + d(k-2) ) / 45.d0
if ( abs(r(k) - r(k-1)) < 0.0001d0 ) then
dfdxi = r(k)
exit
end if
end do
gradf(i) = dfdxi
end do
! Print out the gradient
write(*,*) " "
write(*,*) " Grad(f(x)) = "
write(*,*) " "
do i = 1,n
write(*,*) gradf(i)
end do
end program gradient
Run Code Online (Sandbox Code Playgroud)
在单精度下它运行良好并给了我不错的结果。但是,当我尝试更改为代码中所示的双精度时,在尝试编译时出现错误,声称赋值语句
d(j) = ( f(n,xhup) - f(n,xhdown) ) / (2.d0*h)
Run Code Online (Sandbox Code Playgroud)
正在产生类型不匹配real(4)/real(8)。我尝试了几种不同的双精度声明,并在代码中附加了每个适当的双精度常量d0,每次都会得到相同的错误。我有点困惑该函数如何f产生单精度数字。
问题是 f 没有在主程序中明确定义,因此它被隐式假定为单精度,即 gfortran 的 real(4) 类型。
我完全同意 High Performance Mark 的评论,您确实应该implicit none在所有 Fortran 代码中使用它,以确保所有对象都被明确声明。这样,您将获得有关 f 未明确定义的更合适的错误消息。
另外,您还可以考虑两件事:
在模块中定义函数并将该模块导入主程序中。最好仅在模块内定义所有子例程/函数,以便编译器在调用函数时可以对参数的数量和类型进行额外检查。
您可以(再次在模块中)引入一个精度常数并在必须指定实数类型的任何地方使用它。以下面的示例为例,仅更改行
integer, parameter :: dp = kind(1.0d0)
Run Code Online (Sandbox Code Playgroud)
进入
integer, parameter :: dp = kind(1.0)
Run Code Online (Sandbox Code Playgroud)
您会将所有实际变量从双精度更改为单精度。另请注意_dp文字常量的后缀而不是d0后缀,这也会自动调整其精度。
module accuracy
implicit none
integer, parameter :: dp = kind(1.0d0)
end module accuracy
module myfunc
use accuracy
implicit none
contains
function f(n,x)
integer :: n
real(dp) :: x(n), f
f = 0.5_dp * x(1)**5 + cos(x(2)) + log(x(3)) - sqrt(x(4))
end function f
end module myfunc
program gradient
use myfunc
implicit none
real(dp) :: x(n), xhup(n), xhdown(n), d(M), r(M), dfdxi, h0, h, gradf(n)
:
end program gradient
Run Code Online (Sandbox Code Playgroud)