y_my_own_da.f90

This is a toy "DA" package which treats third order polynomials in two variables. Download by clicking here

It purpose is to illustrate the fundamentals of forward automatic differentiation which is what Berz DA-package does.

It is also written to explain how operator overloading works.

The actual original package of Berz is far more complex and sophisticated. In addition its overloading requires certain technical tricks to deal with the fact that Berz's original Fortran77 deals with integer pointers to a big array where he stores all his polynomials. (FPP overloads the LBL version of Berz's package, and not the thing in COSY-INFINITY

Berz's original paper in Particle Accelerator remains a good description of his DA-package. Later papers contain too many abstractions with little or no examples. However I can provide example programs to illustrate Berz's abstractions ("tower of ideals, nilpotency, λ depth, contractions, Levi-Cevita fields, etc...") Many of these concepts are abstractions of trivialities and serve little purpose in optics even though they may be of general mathematical interest. Therefore do not get distracted by this verbiage. We are doing simple algebra on polynomials and representing functions by polynomials on the computer.

MODULE MY_OWN_DA
IMPLICIT NONE
INTEGER,PARAMETER:: DP=8
PRIVATE INPUT_REAL_IN_MY_TAYLOR
INTEGER :: N_TPSA_EXP = 10
INTEGER :: MY_ORDER = 0

TYPE MY_TAYLOR
REAL(DP) A(0:9)
Taylor = a₀+a₁x₁+a ₂x₂+a₃x₁²+ a₄x₂²+a₅x₁ x₂+a₆x₁³ +a₇x₂³+a₈x₁²x₂+a₉x₁x₂²
END TYPE MY_TAYLOR

INTERFACE ASSIGNMENT (=)
MODULE PROCEDURE INPUT_MY_TAYLOR_IN_REAL → Taylor=real
MODULE PROCEDURE INPUT_REAL_IN_MY_TAYLOR → real=Taylor
END INTERFACE
N.B. In this toy DA, Taylor=Taylor is not needed because Fortran90 does automatically the correct job. This is not true when we overload a pointer as in the case of Berz's package.

INTERFACE OPERATOR (*)
Taylor=Taylor*Taylor is interesting click below
           MODULE PROCEDURE MUL → Taylor * Taylor
           MODULE PROCEDURE DMULSC → Taylor * Real(dp)
           MODULE PROCEDURE DSCMUL → Real(dp) * Taylor
END INTERFACE

INTERFACE OPERATOR (/)
Taylor=Taylor/Taylor is also interesting because first case separating TPSA from DA part to create a nilpotent operator! click below
           MODULE PROCEDURE DIV → Taylor / Taylor
           MODULE PROCEDURE DDIVSC → Taylor / Real(dp)
           MODULE PROCEDURE DSCDIV → Real(dp) / Taylor
           MODULE PROCEDURE IDIVSC → Taylor / Integer → added because useful in example code
END INTERFACE

INTERFACE OPERATOR (+)
           MODULE PROCEDURE ADD → Taylor + Taylor
           MODULE PROCEDURE UNARYADD → +Taylor
           MODULE PROCEDURE DADDSC → Taylor + Real(dp)
           MODULE PROCEDURE DSCADD → Real(dp) + Taylor
END INTERFACE

INTERFACE OPERATOR (-)
           MODULE PROCEDURE SUBS → Taylor - Taylor
           MODULE PROCEDURE UNARYSUB → -Taylor
           MODULE PROCEDURE DSUBSC → Taylor - Real(dp)
           MODULE PROCEDURE DSCSUB → Real(dp) - Taylor
END INTERFACE

INTERFACE OPERATOR (**)
MODULE PROCEDURE POW → Taylor ** Integer
END INTERFACE

Overloading standard procedures

INTERFACE EXP
MODULE PROCEDURE DEXPT → exp(Taylor)
END INTERFACE

INTERFACE LOG
MODULE PROCEDURE DLOGT → log(Taylor)
END INTERFACE

INTERFACE SQRT
MODULE PROCEDURE DSQRTT → sqrt(Taylor)
END INTERFACE

INTERFACE COS
MODULE PROCEDURE DCOST → cos(Taylor)
END INTERFACE

INTERFACE SIN
MODULE PROCEDURE DSINT → sin(Taylor)
END INTERFACE

INTERFACE TPSA_EXP
MODULE PROCEDURE TPSA_EXPT → exp(Taylor) using an infinite sum
END INTERFACE

INTERFACE PRINT
MODULE PROCEDURE PRINT_MY_TAYLOR → Fancy print for human eye
MODULE PROCEDURE PRINT_MY_REAL
END INTERFACE

User defined operator

INTERFACE OPERATOR (.MONO.)
MODULE PROCEDURE MONOT → Creates the monomial r x_i
END INTERFACE

INTERFACE OPERATOR (.SUB.)
MODULE PROCEDURE SUBT → Creates the monomial r x_i
END INTERFACE

INTERFACE ALLOC
MODULE PROCEDURE ALLOC_MY_TAYLOR → for compatibility with FPP
END INTERFACE

INTERFACE KILL
MODULE PROCEDURE ALLOC_MY_TAYLOR → for compatibility with FPP
END INTERFACE

INTERFACE INIT
MODULE PROCEDURE INIT_FAKE → for compatibility with FPP
END INTERFACE

INTERFACE MORPH
MODULE PROCEDURE MORPHT → for compatibility with FPP
END INTERFACE

CONTAINS

SUBROUTINE ALLOC_MY_TAYLOR( S2 )
      IMPLICIT NONE
      TYPE (MY_TAYLOR), INTENT (INOUT) :: S2

      S2%A=0.0_DP

END SUBROUTINE ALLOC_MY_TAYLOR

SUBROUTINE INIT_FAKE(NO,ND,NP,NDPT )
      IMPLICIT NONE
      INTEGER NO,ND,NP,NDPT


END SUBROUTINE INIT_FAKE

FUNCTION MORPHT( T )
      IMPLICIT NONE
      TYPE (MY_TAYLOR) MORPHT
      TYPE (MY_TAYLOR), INTENT(IN) :: T

      MORPHT= T

END FUNCTION MORPHT

FUNCTION MONOT( R, I )
      IMPLICIT NONE
      TYPE (MY_TAYLOR) MONOT
      REAL(DP), INTENT(IN) :: R
      INTEGER, INTENT(IN) :: I

      MONOT=0.0_DP
      MONOT%A(I)= R

END FUNCTION MONOT

FUNCTION SUBT( T, I )
      IMPLICIT NONE
      TYPE (MY_TAYLOR) SUBT
      TYPE (MY_TAYLOR), INTENT(IN) :: T
      INTEGER, INTENT(IN) :: I(2)
      INTEGER K

      SUBT=0.0_DP

      K=I(1)+10*I(2)

      SELECT CASE(K)
      CASE(0)
      SUBT=T%A(0)
      CASE(1)
      SUBT=T%A(1)
      CASE(10)
      SUBT=T%A(2)
      CASE(2)
      SUBT=T%A(3)
      CASE(20)
      SUBT=T%A(4)
      CASE(11)
      SUBT=T%A(5)
      CASE(3)
      SUBT=T%A(6)
      CASE(30)
      SUBT=T%A(7)
      CASE(12)
      SUBT=T%A(8)
      CASE(21)
      SUBT=T%A(9)

      END SELECT

END FUNCTION SUBT

FUNCTION ADD( S1, S2 )
      IMPLICIT NONE
      TYPE (MY_TAYLOR) ADD
      TYPE (MY_TAYLOR), INTENT (IN) :: S1, S2

      ADD%A=S1%A + S2%A

END FUNCTION ADD

FUNCTION DADDSC( S1, SC )
      IMPLICIT NONE
      TYPE (MY_TAYLOR) DADDSC
      TYPE (MY_TAYLOR), INTENT (IN) :: S1
      REAL(DP), INTENT (IN) :: SC

      DADDSC=S1
      DADDSC%A(0)= S1%A(0) + SC

END FUNCTION DADDSC

FUNCTION DSCADD( SC , S1)
      IMPLICIT NONE
      TYPE (MY_TAYLOR) DSCADD
      TYPE (MY_TAYLOR), INTENT (IN) :: S1
      REAL(DP), INTENT (IN) :: SC

      DSCADD=S1
      DSCADD%A(0)= S1%A(0) + SC

END FUNCTION DSCADD

FUNCTION UNARYADD( S1 )
      IMPLICIT NONE
      TYPE (MY_TAYLOR) UNARYADD
      TYPE (MY_TAYLOR), INTENT (IN) :: S1

      UNARYADD=S1

END FUNCTION UNARYADD

FUNCTION SUBS( S1, S2 )
      IMPLICIT NONE
      TYPE (MY_TAYLOR) SUBS
      TYPE (MY_TAYLOR), INTENT (IN) :: S1, S2

      SUBS%A=S1%A - S2%A

END FUNCTION SUBS

FUNCTION UNARYSUB( S1 )
      IMPLICIT NONE
      TYPE (MY_TAYLOR) UNARYSUB
      TYPE (MY_TAYLOR), INTENT (IN) :: S1

      UNARYSUB%A=-S1%A

END FUNCTION UNARYSUB

FUNCTION DSUBSC( S1, SC )
      IMPLICIT NONE
      TYPE (MY_TAYLOR) DSUBSC
      TYPE (MY_TAYLOR), INTENT (IN) :: S1
      REAL(DP), INTENT (IN) :: SC

      DSUBSC=S1
      DSUBSC%A(0)= S1%A(0) - SC

END FUNCTION DSUBSC

FUNCTION DSCSUB( SC , S1 )
      IMPLICIT NONE
      TYPE (MY_TAYLOR) DSCSUB
      TYPE (MY_TAYLOR), INTENT (IN) :: S1
      REAL(DP), INTENT (IN) :: SC

      DSCSUB=-S1
      DSCSUB=SC + DSCSUB

END FUNCTION DSCSUB

FUNCTION MUL( S1, S2 )
      IMPLICIT NONE
      TYPE (MY_TAYLOR) MUL
      TYPE (MY_TAYLOR), INTENT (IN) :: S1, S2

T₁ =a₀+a₁x₁+a₂x₂+a₃x₁²+ a₄x₂²+a₅x₁ x₂+a₆x₁³ +a₇x₂³+a₈x₁²x₂+a₉x₁x₂² = (a₀,a₁,a₂,a₃,a₄,a₅,a₆,a₇,a₈,a₉)
T₂ =b₀+b₁x₁+b₂x₂+b₃x₁²+ b₄x₂²+b₅x₁ x₂+b₆x₁³ +b₇x₂³+b₈x₁²x₂+b₉x₁x₂² = (b₀,b₁,b₂,b₃,b₄,b₅,b₆,b₇,b₈,b₉)

T₁ * T₂ =a₀b₀+{a₀b₁+b₀ a₁}x₁+{a₀b₂+b₀ a₂}x₂+{a₀b₃+b₀ a₃+a₁²x₁²}x₁² +{a₀b₄+b₀a₄+a₂²}x₂²+{a₀b₅+b₀a₅+a₁b₂+b₁a₂}x₁x₂+ ...
= (a₀b₀ , a₀b₁+b₀a₁, a₀b₂+b₀a₂ , a₀b₃+b₀a₃+a₁²x₁² , a₀b₄+b₀a₄+a₂² , a₀b₅+b₀a₅+a₁b₂+b₁a₂ ,...)
Back up

      MUL%A(0)=S1%A(0)*S2%A(0)
      MUL%A(1)=S1%A(1)*S2%A(0)+S1%A(0)*S2%A(1)
      MUL%A(2)=S1%A(2)*S2%A(0)+S1%A(0)*S2%A(2)
      MUL%A(3)=S1%A(3)*S2%A(0)+S1%A(0)*S2%A(3)+S1%A(1)*S2%A(1)
      MUL%A(4)=S1%A(4)*S2%A(0)+S1%A(0)*S2%A(4)+S1%A(2)*S2%A(2)
      MUL%A(5)=S1%A(5)*S2%A(0)+S1%A(0)*S2%A(5)+S1%A(1)*S2%A(2)+S1%A(2)*S2%A(1)
      MUL%A(6)=S1%A(6)*S2%A(0)+S1%A(0)*S2%A(6)+S1%A(1)*S2%A(3)+S1%A(3)*S2%A(1)
      MUL%A(7)=S1%A(7)*S2%A(0)+S1%A(0)*S2%A(7)+S1%A(2)*S2%A(4)+S1%A(4)*S2%A(2)
      MUL%A(8)=S1%A(8)*S2%A(0)+S1%A(0)*S2%A(8)+S1%A(1)*S2%A(5)+S1%A(5)*S2%A(1)+S1%A(2)*S2%A(3)+S1%A(3)*S2%A(2)
      MUL%A(9)=S1%A(9)*S2%A(0)+S1%A(0)*S2%A(9)+S1%A(1)*S2%A(4)+S1%A(4)*S2%A(1)+S1%A(2)*S2%A(5)+S1%A(5)*S2%A(2)

END FUNCTION MUL

FUNCTION DMULSC( S1, SC )
      IMPLICIT NONE
      TYPE (MY_TAYLOR) DMULSC
      TYPE (MY_TAYLOR), INTENT (IN) :: S1
      REAL(DP), INTENT (IN) :: SC

      DMULSC%A= S1%A*SC

END FUNCTION DMULSC

FUNCTION DSCMUL( SC ,S1 )
      IMPLICIT NONE
      TYPE (MY_TAYLOR) DSCMUL
      TYPE (MY_TAYLOR), INTENT (IN) :: S1
      REAL(DP), INTENT (IN) :: SC

      DSCMUL%A= S1%A*SC

END FUNCTION DSCMUL

FUNCTION DDIVSC( S1, SC )
      IMPLICIT NONE
      TYPE (MY_TAYLOR) DDIVSC
      TYPE (MY_TAYLOR), INTENT (IN) :: S1
      REAL(DP), INTENT (IN) :: SC

      DDIVSC%A= S1%A/SC

END FUNCTION DDIVSC

FUNCTION IDIVSC( S1, SC )
      IMPLICIT NONE
      TYPE (MY_TAYLOR) IDIVSC
      TYPE (MY_TAYLOR), INTENT (IN) :: S1
      INTEGER, INTENT (IN) :: SC

      IDIVSC%A= S1%A/SC

END FUNCTION IDIVSC

FUNCTION POW( S1,N )
      IMPLICIT NONE
      TYPE (MY_TAYLOR) POW , T
      TYPE (MY_TAYLOR), INTENT (IN) :: S1
      INTEGER, INTENT (IN) :: N
     INTEGER I

      POW=1.D0

      IF(N<=0) THEN
      T=1.D0/S1
      DO I=1,-N
      POW=POW*T
      ENDDO
     ELSE
      DO I=1,N
      POW=POW*S1
      ENDDO

     ENDIF

END FUNCTION POW

FUNCTION INV( S1 )
      IMPLICIT NONE
      TYPE (MY_TAYLOR) INV,T
      TYPE (MY_TAYLOR), INTENT (IN) :: S1
T₁^-1 = (a₀+a₁x₁+a₂x₂ +a₃x₁²+ a₄x₂²+a₅x₁x₂+ a₆x₁³ +a₇x₂³+a₈x₁²x₂+a₉x₁x₂²)^-1
= a₀^-1(1 + a_1/a₀x₁+a_2/a₀x₂+a_3/a₀x₁²+ a_4/a₀x₂²+a_5/a₀x₁x₂ +a_6/a₀x₁³ +a_7/a₀x₂³+a_8/a₀x₁²x₂+a_9/a₀x₁x₂²)^-1
= a₀^-1 (1 + T)^-1 where T is nilpotent, i.e., T⁴=0!
Therefore (1+T)^-1= 1-T+T²-T³ exactly in our algebra!
Back up
     T=S1/S1%A(0)
     T%A(0)=0.D0

     INV=1.D0-T+T*T-T*T*T
      INV=INV/S1%A(0)

END FUNCTION INV

FUNCTION DIV( S1, S2 )
      IMPLICIT NONE
      TYPE (MY_TAYLOR) DIV
      TYPE (MY_TAYLOR), INTENT (IN) :: S1, S2
Go to function INV
      DIV=INV(S2)
      DIV=S1*DIV
END FUNCTION DIV

FUNCTION DSCDIV( SC , S1)
      IMPLICIT NONE
      TYPE (MY_TAYLOR) DSCDIV
      TYPE (MY_TAYLOR), INTENT (IN) :: S1
      REAL(DP), INTENT (IN) :: SC

      DSCDIV=INV(S1)
      DSCDIV= SC * DSCDIV

END FUNCTION DSCDIV

SUBROUTINE INPUT_REAL_IN_MY_TAYLOR( S2, S1 )
      IMPLICIT NONE
      REAL(DP), INTENT (IN) :: S1
      TYPE (MY_TAYLOR), INTENT (INOUT) :: S2

      S2%A=0.0_DP
      S2%A(0)=S1

END SUBROUTINE INPUT_REAL_IN_MY_TAYLOR

SUBROUTINE INPUT_MY_TAYLOR_IN_REAL( S2, S1 )
      IMPLICIT NONE
      TYPE (MY_TAYLOR), INTENT (IN) :: S1
      REAL(DP), INTENT (INOUT) :: S2

      S2=S1%A(0)

END SUBROUTINE INPUT_MY_TAYLOR_IN_REAL

FUNCTION DEXPT( S1 )
      IMPLICIT NONE
      TYPE (MY_TAYLOR) DEXPT,T,TT
      TYPE (MY_TAYLOR), INTENT (IN) :: S1
      INTEGER I,NO
     T=S1
     T%A(0)=0.0_DP

      NO=3

     DEXPT=1.0_DP
     TT=1.0_DP

      DO I=1,NO
      TT=TT*T/I
      DEXPT=DEXPT + TT
      ENDDO

      DEXPT=DEXPT*EXP(S1%A(0))

END FUNCTION DEXPT

FUNCTION DLOGT( S1 )
      IMPLICIT NONE
      TYPE (MY_TAYLOR) DLOGT,T
      TYPE (MY_TAYLOR), INTENT (IN) :: S1

     T=S1/S1%A(0)
     T%A(0)=0.0_DP

     DLOGT=T-T**2/2.0_DP+T**3/3.0_DP
      DLOGT=DLOGT+ LOG(S1%A(0))

END FUNCTION DLOGT

FUNCTION DSQRTT( S1 )
      IMPLICIT NONE
      TYPE (MY_TAYLOR) DSQRTT,T
      TYPE (MY_TAYLOR), INTENT (IN) :: S1

     T=S1/S1%A(0)
     T%A(0)=0.0_DP

     DSQRTT=1.0_DP + T/2.0_DP - T**2/8.0_DP+T**3/16.0_DP
      DSQRTT=DSQRTT* SQRT(S1%A(0))

END FUNCTION DSQRTT

FUNCTION DCOST(S1)
      IMPLICIT NONE
      TYPE (MY_TAYLOR) DCOST,T
      TYPE (MY_TAYLOR), INTENT (IN) :: S1

     T=S1
     T%A(0)=0.0_DP

     DCOST=COS(S1%A(0))*(1.0_DP-T**2/2.0_DP)-SIN(S1%A(0))*(T-T**3/6.0_DP)

END FUNCTION DCOST

FUNCTION DSINT(S1)
      IMPLICIT NONE
      TYPE (MY_TAYLOR) DSINT,T
      TYPE (MY_TAYLOR), INTENT (IN) :: S1

     T=S1
     T%A(0)=0.0_DP

     DSINT=SIN(S1%A(0))*(1.0_DP-T**2/2.0_DP)+COS(S1%A(0))*(T-T**3/6.0_DP)

END FUNCTION DSINT

FUNCTION TPSA_EXPT( S1 )
      IMPLICIT NONE
      TYPE (MY_TAYLOR) TPSA_EXPT,T
      TYPE (MY_TAYLOR), INTENT (IN) :: S1
      INTEGER I


      T=1.0_DP
      TPSA_EXPT=1.0_DP

      DO I=1,N_TPSA_EXP
      T=S1*T/I
      TPSA_EXPT=TPSA_EXPT+T
      ENDDO

END FUNCTION TPSA_EXPT

SUBROUTINE PRINT_MY_REAL( S , MF)
      IMPLICIT NONE
      REAL(DP) S
      INTEGER MF

     101 FORMAT(3(3X,A7,E20.13))
     100 FORMAT(3X,A7,E20.13)

      WRITE(MF,100) " (0,0) ", S

END SUBROUTINE PRINT_MY_REAL

SUBROUTINE PRINT_MY_TAYLOR( S , MF)
      IMPLICIT NONE
      TYPE (MY_TAYLOR) S
      INTEGER MF

     101 FORMAT(3(3X,A7,E20.13))
     100 FORMAT(3X,A7,E20.13)

      WRITE(MF,100) " (0,0) ", S%A(0)
      IF(MY_ORDER>=1) WRITE(MF,101) " (1,0) ", S%A(1)
      IF(MY_ORDER>=1) WRITE(MF,101) " (0,1) ", S%A(2)
      IF(MY_ORDER>=2) WRITE(MF,100) " (2,0) ", S%A(3)
      IF(MY_ORDER>=2) WRITE(MF,100) " (0,2) ", S%A(4)
      IF(MY_ORDER>=2) WRITE(MF,100) " (1,1) ", S%A(5)
      IF(MY_ORDER>=3) WRITE(MF,100) " (3,0) ", S%A(6)
      IF(MY_ORDER>=3) WRITE(MF,100) " (0,3) ", S%A(7)
      IF(MY_ORDER>=3) WRITE(MF,100) " (2,1) ", S%A(8)
      IF(MY_ORDER>=3) WRITE(MF,100) " (1,2) ", S%A(9)

END SUBROUTINE PRINT_MY_TAYLOR

END MODULE MY_OWN_DA