Fortran Module for Dual Numbers



! DUAL Dual numbers module for automatic differentiation.

!

! To do automatic differentiation we need a number d such that d^2 == 0.

! Then by a simple application of Taylor's theorem we have

!

!     f(x + d) = f(x) + f'(x) d

!

! with all higher-order terms vanishing. Therefore applying a function to a dual

! number x + d generates both the value f(x) *and* the derivative f'(x).

!

! We define a dual number module that overrides the basic arithmetic functions:

!

!     (x1 + y1 d) + (x2 + y2 d) = (x1 + x2) + (y1 + y2) d

!

!     (x1 + y1 d) * (x2 + y2 d) = (x1 * x2) + (x1 * y2 + y1 * x2) d

!

! We then make use of linearity of the derivative,

!

!     (f+g)'(x) = f'(x) + g'(x)

!

! Leibniz' rule for products,

!

!     (fg)'(x) = f(x) g'(x) + f'(x) g(x)

!

! and the chain rule for function composition.

!

!     (f.g)'(x) = g'(x) f'(g(x))

!

! the rule for composition is similar to multiplication of dual numbers:

!

!     f(g(x+d)) = f(g(x) + g'(x) d) = f(g(x)) + f'(g(x)) g'(x) d

!

!

! Author:Sukhbinder

!



        MODULE DUAL

        IMPLICIT NONE

        PRIVATE



        TYPE,PUBLIC:: DUAL_NUM

        REAL(8)::x_rl

        REAL(8)::x_di

        END TYPE DUAL_NUM



        PUBLIC OPERATOR (-)

        PUBLIC OPERATOR (+)

        PUBLIC OPERATOR (*)

        PUBLIC OPERATOR (/)

        PUBLIC ABS



        INTERFACE OPERATOR (-)

        MODULE PROCEDURE MINUS_DN

        END INTERFACE



        INTERFACE OPERATOR (+)

        MODULE PROCEDURE ADD_DN

        END INTERFACE



        INTERFACE OPERATOR (*)

        MODULE PROCEDURE MULT_DN

        END INTERFACE



        INTERFACE OPERATOR (/)

        MODULE PROCEDURE DIV_DN

        END INTERFACE



        INTERFACE OPERATOR (**)

        MODULE PROCEDURE POWER_DN

        END INTERFACE



        INTERFACE ABS

        MODULE PROCEDURE ABS_DN

        END INTERFACE





        INTERFACE SQRT

        MODULE PROCEDURE SQRT_DN

        END INTERFACE



        INTERFACE EXP

        MODULE PROCEDURE EXP_DN

        END INTERFACE



        INTERFACE LOG

        MODULE PROCEDURE LOG_DN

        END INTERFACE



        INTERFACE SIN

        MODULE PROCEDURE SIN_DN

        END INTERFACE



        INTERFACE COS

        MODULE PROCEDURE COS_DN

        END INTERFACE



        INTERFACE TAN

        MODULE PROCEDURE TAN_DN

        END INTERFACE

        CONTAINS





        ELEMENTAL FUNCTION POWER_DN(u,v) RESULT(res)

        TYPE (DUAL_NUM), INTENT(IN)::u

        TYPE (DUAL_NUM)::res

        INTEGER,INTENT(IN) :: v

        res%x_rl = (u%x_rl)**v

        res%x_di= v * u%x_di * (u%x_rl)**(v-1)

        END FUNCTION POWER_DN



        ELEMENTAL FUNCTION MINUS_DN(u,v) RESULT(res)

        TYPE (DUAL_NUM), INTENT(IN)::u,v

        TYPE (DUAL_NUM)::res

        res%x_rl = u%x_rl-v%x_rl

        res%x_di= u%x_di-v%x_di

        END FUNCTION MINUS_DN







        ELEMENTAL FUNCTION ADD_DN(u,v) RESULT(res)

        TYPE (DUAL_NUM), INTENT(IN)::u,v

        TYPE (DUAL_NUM)::res

        res%x_rl = u%x_rl+v%x_rl

        res%x_di= u%x_di+v%x_di

        END FUNCTION ADD_DN





        ELEMENTAL FUNCTION MULT_DN(u,v) RESULT(res)

        TYPE (DUAL_NUM), INTENT(IN)::u,v

        TYPE (DUAL_NUM)::res

        res%x_rl = u%x_rl*v%x_rl

        res%x_di= u%x_di*v%x_rl + u%x_rl*v%x_di

        END FUNCTION MULT_DN





        ELEMENTAL FUNCTION DIV_DN(u,v) RESULT(res)

        TYPE (DUAL_NUM), INTENT(IN)::u,v

        REAL(8)::tmp

        TYPE (DUAL_NUM)::res

        tmp=1.D0/v%x_rl

        res%x_rl = u%x_rl*tmp

        res%x_di =(u%x_di- res%x_rl*v%x_di)*tmp

        END FUNCTION DIV_DN



        ELEMENTAL FUNCTION ABS_DN(u) RESULT(res)

        TYPE (DUAL_NUM), INTENT(IN)::u

        TYPE (DUAL_NUM)::res

        res%x_rl = ABS(u%x_rl)

        res%x_di = SIGN(abs(u%x_di),u%x_rl)

        END FUNCTION ABS_DN





        ELEMENTAL FUNCTION SQRT_DN(u) RESULT(res)

        TYPE (DUAL_NUM), INTENT(IN)::u

        TYPE (DUAL_NUM)::res

        res%x_rl = SQRT(u%x_rl)

        res%x_di = u%x_di /(2* SQRT(u%x_rl))

        END FUNCTION SQRT_DN



        ELEMENTAL FUNCTION EXP_DN(u) RESULT(res)

        TYPE (DUAL_NUM), INTENT(IN)::u

        TYPE (DUAL_NUM)::res

        res%x_rl = EXP(u%x_rl)

        res%x_di = u%x_di * EXP(u%x_rl)

        END FUNCTION EXP_DN



        ELEMENTAL FUNCTION LOG_DN(u) RESULT(res)

        TYPE (DUAL_NUM), INTENT(IN)::u

        TYPE (DUAL_NUM)::res

        res%x_rl = LOG(u%x_rl)

        res%x_di = u%x_di / u%x_rl

        END FUNCTION LOG_DN



        ELEMENTAL FUNCTION SIN_DN(u) RESULT(res)

        TYPE (DUAL_NUM), INTENT(IN)::u

        TYPE (DUAL_NUM)::res

        res%x_rl = SIN(u%x_rl)

        res%x_di = u%x_di * COS(u%x_rl)

        END FUNCTION  SIN_DN



        ELEMENTAL FUNCTION COS_DN(u) RESULT(res)

        TYPE (DUAL_NUM), INTENT(IN)::u

        TYPE (DUAL_NUM)::res

        res%x_rl = COS(u%x_rl)

        res%x_di = -1*u%x_di *sin(u%x_rl)

        END FUNCTION COS_DN



        ELEMENTAL FUNCTION TAN_DN(u) RESULT(res)

        TYPE (DUAL_NUM), INTENT(IN)::u

        TYPE (DUAL_NUM)::res

        res%x_rl = TAN(u%x_rl)

        res%x_di = u%x_di * sec(u%x_rl)**2

        END FUNCTION TAN_DN





        END MODULE DUAL

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s