[原创] Fortran 实现一阶自动微分

+,-,*,/,**

sin,cos,log,exp

x%dv=1.0

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