Fortran 实现一阶自动微分
本帖最后由 Transpose 于 2021-6-2 20:51 编辑一阶自动微分程序，重载了+,-,*,/,**以及简单的数学函数sin,cos,log,exp程序使用前初始化 x%dv=1.0,因为数学上有 dx/dx=1
数据类型使用了浮点数，所以不支持整数的运算，例如
f(x)=x+2 !no
f(x)=x+2.d0 !yes
f(x)=2.d0+x !yes
希望大家纠错
module deriv_mod
!! autoderiv
implicit none
integer,parameter::db=selected_real_kind(15, 307)
type deriv
real(kind=db)::v
!! value
real(kind=db)::dv
!! 1st deriv
contains
generic::operator(+)=> add_dd,add_dn,add_nd
generic::operator(-)=> sub_dd,sub_dn,sub_nd
generic::operator(*)=> mult_dd,mult_dn,mult_nd
generic::operator(/)=> div_dd,div_dn,div_nd
generic::operator(**) => power_dd,power_dn,power_nd
procedure,private,pass(this)::add_dd,add_dn,add_nd
procedure,private,pass(this)::sub_dd,sub_dn,sub_nd
procedure,private,pass(this)::mult_dd,mult_dn,mult_nd
procedure,private,pass(this)::div_dd,div_dn,div_nd
procedure,private,pass(this)::power_dd,power_dn,power_nd
end type
interface sin
module procedure sin_d
end interface
interface cos
module procedure cos_d
end interface
interface exp
module procedure exp_d
end interface
interface log
module procedure log_d
end interface
contains
function add_dd(this, x) result(y)
!! f(x)+h(x)
implicit none
class(deriv), intent(in)::this
class(deriv), intent(in)::x
type(deriv)::y
y%v = this%v+x%v
y%dv = this%dv+x%dv
end function add_dd
function add_dn(this, x) result(y)
!! f(x)+a
implicit none
class(deriv), intent(in)::this
real(kind=db), intent(in)::x
type(deriv)::y
y%v = this%v+x
y%dv = this%dv
end function add_dn
function add_nd(x, this) result(y)
!! a+f(x)
implicit none
class(deriv), intent(in)::this
real(kind=db), intent(in)::x
type(deriv)::y
y%v = this%v+x
y%dv = this%dv
end function add_nd
function mult_dd(this, x) result(y)
!! f(x)*h(x)
implicit none
class(deriv), intent(in)::this
class(deriv), intent(in)::x
type(deriv)::y
y%v = this%v*x%v
y%dv = this%dv*x%v+this%v*x%dv
end function mult_dd
function mult_dn(this, x) result(y)
!! f(x)*a
implicit none
class(deriv), intent(in)::this
real(kind=db), intent(in)::x
type(deriv)::y
y%v = this%v*x
y%dv = this%dv*x
end function mult_dn
function mult_nd(x, this) result(y)
!! a*f(x)
implicit none
class(deriv), intent(in)::this
real(kind=db), intent(in)::x
type(deriv)::y
y%v = this%v*x
y%dv = this%dv*x
end function mult_nd
function sub_dd(this, x) result(y)
!! f(x)-h(x)
implicit none
class(deriv), intent(in)::this
class(deriv), intent(in)::x
type(deriv)::y
y%v = this%v-x%v
y%dv = this%dv-x%dv
end function sub_dd
function sub_dn(this, x) result(y)
!! f(x)-a
implicit none
class(deriv), intent(in)::this
real(kind=db), intent(in)::x
type(deriv)::y
y%v = this%v-x
y%dv = this%dv
end function sub_dn
function sub_nd(x ,this) result(y)
!! a-f(x)
implicit none
class(deriv), intent(in)::this
real(kind=db), intent(in)::x
type(deriv)::y
y%v = x-this%v
y%dv = -this%dv
end function sub_nd
function div_dd(this, x) result(y)
!! f(x)/h(x)
implicit none
class(deriv), intent(in)::this
class(deriv), intent(in)::x
type(deriv)::y
y%v = this%v/x%v
y%dv = (this%dv*x%v-this%v*x%dv)/(x%v)**2
end function div_dd
function div_dn(this, x) result(y)
!! f(x)/a
implicit none
class(deriv), intent(in)::this
real(kind=db), intent(in)::x
type(deriv)::y
y%v = this%v/x
y%dv = this%dv/x
end function div_dn
function div_nd(x, this) result(y)
!! a/f(x)
implicit none
class(deriv), intent(in)::this
real(kind=db), intent(in)::x
type(deriv)::y
y%v = x/this%v
y%dv =-this%dv*x/(this%v)**2
end function div_nd
function power_dd(this, x) result(y)
!! f(x)**h(x)
implicit none
class(deriv), intent(in)::this
class(deriv), intent(in)::x
type(deriv)::y
!y%v = (this%v)**(x%v)
!y%dv = (x%dv*log(this%v)+x%v/this%v*this%dv)*y%v
y=exp(x*log(this))
end function power_dd
function power_dn(this, x) result(y)
!! f(x)**a
implicit none
class(deriv), intent(in)::this
real(kind=db), intent(in)::x
type(deriv)::y
!y%v = (this%v)**(x)
!y%dv = x*(this%v)**(x-1.d0)*(this%dv)
y=exp(x*log(this))
end function power_dn
function power_nd(x, this) result(y)
!! a**f(x)
implicit none
class(deriv), intent(in)::this
real(kind=db), intent(in)::x
type(deriv)::y
!y%v = x**(this%v)
!y%dv= log(x)*this%dv*y%v
y=exp(this*log(x))
end function power_nd
function exp_d(x) result(z)
!! exp(x)
class(deriv), intent(in)::x
type(deriv)::z
z%v = exp(x%v)
z%dv = exp(x%v)*(x%dv)
end function exp_d
function sin_d(x) result(z)
!! sin(x)
class(deriv), intent(in)::x
type(deriv)::z
z%v = sin(x%v)
z%dv = cos(x%v)*(x%dv)
end function sin_d
function cos_d(x) result(z)
!! cos(x)
class(deriv), intent(in)::x
type(deriv)::z
z%v = cos(x%v)
z%dv = -sin(x%v)*(x%dv)
end function cos_d
function log_d(x) result(z)
!! log(x)
class(deriv), intent(in)::x
type(deriv)::z
z%v = log(x%v)
z%dv = 1._db/(x%v)*(x%dv)
end function log_d
end module
program main
use deriv_mod
implicit none
type(deriv),external::func
type(deriv)::x,y
! set dx=1, for dx/dx=1
x=deriv(3.d0,1.d0)
y=func(x)
write(*,*)"f(3.d0)=",y%v
write(*,*)"df/dx|(3.d0)=",y%dv
end program
function func(x) result(f1)
use deriv_mod
type(deriv),intent(in)::x
type(deriv)::f1
f1= x**sin(x)/log(x)-exp(cos(x))+3.d0
end function func
输出结果
f(3.d0)= 3.6913070496183309
df/dx|(3.d0)=-1.3760726245735238
数值严格结果
f(3)=3.691307049618331
df/dx(3)=-1.376072624573524
参考文献
. Modern Fortran in Practice,Arjen Markus, Deltares, Delft, The Netherlands, Chapter 3
这个东西有点意思~~很有直接写数学式子的感觉~~发挥了For tran这个名字的“特点” 这是个很蛋疼的事情……
我干过，比你这个更蛋疼的
用矩阵的形式表达偏导数，也就是说数据类型可以自定义有任何多个自变量的求一阶偏导，
支持单精度和双精度的求导，类型支持和单精度，双精度以及整型的混合运算
我干完这个之后看到神经网络里的反向自动微分，然后就把自己写的东西扔进了垃圾堆
主要是像在Fortran里实现NN里这种，用流程图实现的反向自动微分基本不可能，先天劣势太大了……
而实际大量的使用过程中，没有反向自动微分，太容易中间量爆炸了，内存和算力都HOLD不住 liudy02 发表于 2021-6-3 17:27
这是个很蛋疼的事情……
我干过，比你这个更蛋疼的
用矩阵的形式表达偏导数，也就是说数据类型可以自定义有 ...
另外，想实现高阶导数对Fortran也还得自己一阶一阶定义
用python等语言中人家做好的AI库里的自动微分功能他不香么……
页:
[1]