Fortran 实现一阶自动微分

Transpose

一阶自动微分程序，重载了

[Fortran] syntaxhighlighter_viewsource syntaxhighlighter_copycode

+,-,*,/,**

以及简单的数学函数

[Fortran] syntaxhighlighter_viewsource syntaxhighlighter_copycode

sin,cos,log,exp

程序使用前初始化

[Fortran] syntaxhighlighter_viewsource syntaxhighlighter_copycode

x%dv=1.0

,因为数学上有 dx/dx=1
数据类型使用了浮点数，所以不支持整数的运算，例如

[Fortran] syntaxhighlighter_viewsource syntaxhighlighter_copycode

f(x)=x+2    !no
f(x)=x+2.d0 !yes
f(x)=2.d0+x !yes

希望大家纠错

[Fortran] syntaxhighlighter_viewsource syntaxhighlighter_copycode

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

输出结果

[Fortran] syntaxhighlighter_viewsource syntaxhighlighter_copycode

 f(3.d0)=   3.6913070496183309
 df/dx|(3.d0)=  -1.3760726245735238

数值严格结果

[Fortran] syntaxhighlighter_viewsource syntaxhighlighter_copycode

f(3)=3.691307049618331
df/dx(3)=-1.376072624573524

参考文献

[1]. Modern Fortran in Practice,Arjen Markus, Deltares, Delft, The Netherlands, Chapter 3

vvt

这个东西有点意思~~很有直接写数学式子的感觉~~发挥了For tran这个名字的“特点”

liudy02

这是个很蛋疼的事情……
我干过，比你这个更蛋疼的
用矩阵的形式表达偏导数，也就是说数据类型可以自定义有任何多个自变量的求一阶偏导，
支持单精度和双精度的求导，类型支持和单精度，双精度以及整型的混合运算
我干完这个之后看到神经网络里的反向自动微分，然后就把自己写的东西扔进了垃圾堆
主要是像在Fortran里实现NN里这种，用流程图实现的反向自动微分基本不可能，先天劣势太大了……
而实际大量的使用过程中，没有反向自动微分，太容易中间量爆炸了，内存和算力都HOLD不住

liudy02

liudy02 发表于 2021-6-3 17:27
这是个很蛋疼的事情……
我干过，比你这个更蛋疼的
用矩阵的形式表达偏导数，也就是说数据类型可以自定义有 ...

另外，想实现高阶导数对Fortran也还得自己一阶一阶定义
用python等语言中人家做好的AI库里的自动微分功能他不香么……