Second Order Adjoints with the NAG Fortran 95 Compiler Uwe Naumann, - - PowerPoint PPT Presentation

▶
second order adjoints with the nag fortran 95 compiler
SMART_READER_LITE
LIVE PREVIEW

Second Order Adjoints with the NAG Fortran 95 Compiler Uwe Naumann, - - PowerPoint PPT Presentation

Oct 30, 2023 667 likes •1.3k views

Second Order Adjoints with the NAG Fortran 95 Compiler Uwe Naumann, Michael Maier RWTH Aachen, Germany Bruce Christianson, Dmitrij Gendler, Jan Riehme University of Hertfordshire, UK riehme@stce.rwth-aachen.de May 22, 2007 Outline COMPAD-II

slide-1
SLIDE 1

Second Order Adjoints with the NAG Fortran 95 Compiler

Uwe Naumann, Michael Maier RWTH Aachen, Germany Bruce Christianson, Dmitrij Gendler, Jan Riehme University of Hertfordshire, UK riehme@stce.rwth-aachen.de May 22, 2007

slide-2
SLIDE 2

Outline

1

COMPAD-II

2

Tangent linear mode by parsetree manipulation

3

Adjoint code generation

4

Second Order Adjoints

6

Outlook

Uwe Naumann, Michael Maier RWTH Aachen, Germany[1ex] Bruce Christianson, Dmitrij Gendler, Jan Riehme University of Hertfordshire, Second Order Adjoints with the NAG Fortran 95 Compiler May 22, 2007 2 / 29

slide-3
SLIDE 3

COMPAD-II

Followup of COMPAD-I project (2002-2004) Funded be EPSRC and QinetiQ Duration : 2006/10 – 2008/10 2 full time (Dima, me) + 1 student (Micheal) Target (A): COMPAD : Compiler AD / Adjoints

robust adjoint compiler based on the NAGWare Fortran 95 compiler simple checkpointing of loops

Target (B) : ADAC : Automatic Differentiation of Assembler Code

differentiate a top-level function given as assembler code control file specifies independent and dependend variables coupling with NAGWare Fortran 95 compiler to get second order adjoints based on source transformation only but this is the sound of the future

Uwe Naumann, Michael Maier RWTH Aachen, Germany[1ex] Bruce Christianson, Dmitrij Gendler, Jan Riehme University of Hertfordshire, Second Order Adjoints with the NAG Fortran 95 Compiler May 22, 2007 3 / 29

slide-4
SLIDE 4

COMPAD-II

Followup of COMPAD-I project (2002-2004) Funded be EPSRC and QinetiQ Duration : 2006/10 – 2008/10 2 full time (Dima, me) + 1 student (Micheal) Target (A): COMPAD : Compiler AD / Adjoints

robust adjoint compiler based on the NAGWare Fortran 95 compiler simple checkpointing of loops

Target (B) : ADAC : Automatic Differentiation of Assembler Code

differentiate a top-level function given as assembler code control file specifies independent and dependend variables coupling with NAGWare Fortran 95 compiler to get second order adjoints based on source transformation only but this is the sound of the future

Uwe Naumann, Michael Maier RWTH Aachen, Germany[1ex] Bruce Christianson, Dmitrij Gendler, Jan Riehme University of Hertfordshire, Second Order Adjoints with the NAG Fortran 95 Compiler May 22, 2007 3 / 29

slide-5
SLIDE 5

COMPAD-II

Followup of COMPAD-I project (2002-2004) Funded be EPSRC and QinetiQ Duration : 2006/10 – 2008/10 2 full time (Dima, me) + 1 student (Micheal) Target (A): COMPAD : Compiler AD / Adjoints

robust adjoint compiler based on the NAGWare Fortran 95 compiler simple checkpointing of loops

Target (B) : ADAC : Automatic Differentiation of Assembler Code

differentiate a top-level function given as assembler code control file specifies independent and dependend variables coupling with NAGWare Fortran 95 compiler to get second order adjoints based on source transformation only but this is the sound of the future

Uwe Naumann, Michael Maier RWTH Aachen, Germany[1ex] Bruce Christianson, Dmitrij Gendler, Jan Riehme University of Hertfordshire, Second Order Adjoints with the NAG Fortran 95 Compiler May 22, 2007 3 / 29

slide-6
SLIDE 6

COMPAD-II — Basic setup

Modified NAGWare Fortran 95 compiler with builtin AD Provide a (couple of) COMPAD MODULE (s) defining a COMPAD TYPE Any source code compiled with the a -ad flag get differentiated, no special top-level routine, no changes required by the user passive driver program, explicitely using a COMPAD MODULE, all active variables have to be of COMPAD TYPE Two modes: Tangent linear mode by compiler generated overloading adjoint mode with compiler generated adjoint code (source trafo) Two flavors each selected by choosing a specific COMPAD MODULE: pure mode with scalar adjoints derived mode where values and adjoints are represented by derived types

Uwe Naumann, Michael Maier RWTH Aachen, Germany[1ex] Bruce Christianson, Dmitrij Gendler, Jan Riehme University of Hertfordshire, Second Order Adjoints with the NAG Fortran 95 Compiler May 22, 2007 4 / 29

slide-7
SLIDE 7

COMPAD-II — Basic setup

Modified NAGWare Fortran 95 compiler with builtin AD Provide a (couple of) COMPAD MODULE (s) defining a COMPAD TYPE Any source code compiled with the a -ad flag get differentiated, no special top-level routine, no changes required by the user passive driver program, explicitely using a COMPAD MODULE, all active variables have to be of COMPAD TYPE Two modes: Tangent linear mode by compiler generated overloading adjoint mode with compiler generated adjoint code (source trafo) Two flavors each selected by choosing a specific COMPAD MODULE: pure mode with scalar adjoints derived mode where values and adjoints are represented by derived types

Uwe Naumann, Michael Maier RWTH Aachen, Germany[1ex] Bruce Christianson, Dmitrij Gendler, Jan Riehme University of Hertfordshire, Second Order Adjoints with the NAG Fortran 95 Compiler May 22, 2007 4 / 29

slide-8
SLIDE 8

COMPAD-II — Basic setup

Modified NAGWare Fortran 95 compiler with builtin AD Provide a (couple of) COMPAD MODULE (s) defining a COMPAD TYPE Any source code compiled with the a -ad flag get differentiated, no special top-level routine, no changes required by the user passive driver program, explicitely using a COMPAD MODULE, all active variables have to be of COMPAD TYPE Two modes: Tangent linear mode by compiler generated overloading adjoint mode with compiler generated adjoint code (source trafo) Two flavors each selected by choosing a specific COMPAD MODULE: pure mode with scalar adjoints derived mode where values and adjoints are represented by derived types

Uwe Naumann, Michael Maier RWTH Aachen, Germany[1ex] Bruce Christianson, Dmitrij Gendler, Jan Riehme University of Hertfordshire, Second Order Adjoints with the NAG Fortran 95 Compiler May 22, 2007 4 / 29

slide-9
SLIDE 9

COMPAD MODULEs and COMPAD TYPEs

For pure mode COMPAD MODULE has to provide Derived type COMPAD TYPE with double precision component val Adjoint mode: COMPAD TYPE with additional drv (double) Overloaded operators and intrinsic functions for the COMPAD TYPE A COMPAD MODULE for the derived mode has to provide COMPAD TYPE with component val of arbitrary VALUE TYPE Overloaded operators and intrinsic functions for the COMPAD TYPE Adjoint mode:

COMPAD TYPE with additional component drv of arbitrary DERIV TYPE Overloaded operators and intrinsic funcitions for the DERIV TYPE

Uwe Naumann, Michael Maier RWTH Aachen, Germany[1ex] Bruce Christianson, Dmitrij Gendler, Jan Riehme University of Hertfordshire, Second Order Adjoints with the NAG Fortran 95 Compiler May 22, 2007 5 / 29

slide-10
SLIDE 10

COMPAD MODULEs and COMPAD TYPEs

For pure mode COMPAD MODULE has to provide Derived type COMPAD TYPE with double precision component val Adjoint mode: COMPAD TYPE with additional drv (double) Overloaded operators and intrinsic functions for the COMPAD TYPE A COMPAD MODULE for the derived mode has to provide COMPAD TYPE with component val of arbitrary VALUE TYPE Overloaded operators and intrinsic functions for the COMPAD TYPE Adjoint mode:

COMPAD TYPE with additional component drv of arbitrary DERIV TYPE Overloaded operators and intrinsic funcitions for the DERIV TYPE

Uwe Naumann, Michael Maier RWTH Aachen, Germany[1ex] Bruce Christianson, Dmitrij Gendler, Jan Riehme University of Hertfordshire, Second Order Adjoints with the NAG Fortran 95 Compiler May 22, 2007 5 / 29

slide-11
SLIDE 11

COMPAD-II — Current status

Overloading mode MINPACK-II , scalar valued problems only

working, that is gradients are correct need to replace statement functions (deprecated feature) with internal subprograms by hand

CNS - code from J¨

  • rg Stiller, TU Dresden

10000 lines of Fortran 95 code uses allocatable arrays and other Fortran 95 features derivatives checked against divided differences Runtime: factor 12 with gcc and factor 4 with Intel C compiler

Adjoint mode we have just started very simple test programs (Rosenbrock function) works just one MINPACK-II problems tested / works (Elastic plastic torsion)

Uwe Naumann, Michael Maier RWTH Aachen, Germany[1ex] Bruce Christianson, Dmitrij Gendler, Jan Riehme University of Hertfordshire, Second Order Adjoints with the NAG Fortran 95 Compiler May 22, 2007 6 / 29

slide-12
SLIDE 12

COMPAD-II — Current status

Overloading mode MINPACK-II , scalar valued problems only

working, that is gradients are correct need to replace statement functions (deprecated feature) with internal subprograms by hand

CNS - code from J¨

  • rg Stiller, TU Dresden

10000 lines of Fortran 95 code uses allocatable arrays and other Fortran 95 features derivatives checked against divided differences Runtime: factor 12 with gcc and factor 4 with Intel C compiler

Adjoint mode we have just started very simple test programs (Rosenbrock function) works just one MINPACK-II problems tested / works (Elastic plastic torsion)

Uwe Naumann, Michael Maier RWTH Aachen, Germany[1ex] Bruce Christianson, Dmitrij Gendler, Jan Riehme University of Hertfordshire, Second Order Adjoints with the NAG Fortran 95 Compiler May 22, 2007 6 / 29

slide-13
SLIDE 13

COMPAD-II — Current status

Overloading mode MINPACK-II , scalar valued problems only

working, that is gradients are correct need to replace statement functions (deprecated feature) with internal subprograms by hand

CNS - code from J¨

  • rg Stiller, TU Dresden

10000 lines of Fortran 95 code uses allocatable arrays and other Fortran 95 features derivatives checked against divided differences Runtime: factor 12 with gcc and factor 4 with Intel C compiler

Adjoint mode we have just started very simple test programs (Rosenbrock function) works just one MINPACK-II problems tested / works (Elastic plastic torsion)

Uwe Naumann, Michael Maier RWTH Aachen, Germany[1ex] Bruce Christianson, Dmitrij Gendler, Jan Riehme University of Hertfordshire, Second Order Adjoints with the NAG Fortran 95 Compiler May 22, 2007 6 / 29

slide-14
SLIDE 14

Tangent linear mode by parsetree manipulation

Parsetree manipulation Tangent linear mode by overloading

Uwe Naumann, Michael Maier RWTH Aachen, Germany[1ex] Bruce Christianson, Dmitrij Gendler, Jan Riehme University of Hertfordshire, Second Order Adjoints with the NAG Fortran 95 Compiler May 22, 2007 7 / 29

slide-15
SLIDE 15

TLM by parsetree manipulation

Specifying the -ad (or the -ad overload) flag of the compiler cause Generation of use - statements for a COMPAD MODULE Activate during parsing: change type of any floating point variables to COMPAD TYPE Note that we have now Association by address (in contrast to Association by name in first sourcetrafo prototype) Resolution of overloading for COMPAD TYPE is business as usual Special care has to be taken for Initialisation expressions: Create structure constructors DATA entries: encapsulate values within COMPAD TYPE constructors during semantic analysis Literal constants in calls must be encapsulated too Deactivation of floating point parameters

Uwe Naumann, Michael Maier RWTH Aachen, Germany[1ex] Bruce Christianson, Dmitrij Gendler, Jan Riehme University of Hertfordshire, Second Order Adjoints with the NAG Fortran 95 Compiler May 22, 2007 8 / 29

slide-16
SLIDE 16

TLM by parsetree manipulation

Specifying the -ad (or the -ad overload) flag of the compiler cause Generation of use - statements for a COMPAD MODULE Activate during parsing: change type of any floating point variables to COMPAD TYPE Note that we have now Association by address (in contrast to Association by name in first sourcetrafo prototype) Resolution of overloading for COMPAD TYPE is business as usual Special care has to be taken for Initialisation expressions: Create structure constructors DATA entries: encapsulate values within COMPAD TYPE constructors during semantic analysis Literal constants in calls must be encapsulated too Deactivation of floating point parameters

Uwe Naumann, Michael Maier RWTH Aachen, Germany[1ex] Bruce Christianson, Dmitrij Gendler, Jan Riehme University of Hertfordshire, Second Order Adjoints with the NAG Fortran 95 Compiler May 22, 2007 8 / 29

slide-17
SLIDE 17

TLM by parsetree manipulation

Specifying the -ad (or the -ad overload) flag of the compiler cause Generation of use - statements for a COMPAD MODULE Activate during parsing: change type of any floating point variables to COMPAD TYPE Note that we have now Association by address (in contrast to Association by name in first sourcetrafo prototype) Resolution of overloading for COMPAD TYPE is business as usual Special care has to be taken for Initialisation expressions: Create structure constructors DATA entries: encapsulate values within COMPAD TYPE constructors during semantic analysis Literal constants in calls must be encapsulated too Deactivation of floating point parameters

Uwe Naumann, Michael Maier RWTH Aachen, Germany[1ex] Bruce Christianson, Dmitrij Gendler, Jan Riehme University of Hertfordshire, Second Order Adjoints with the NAG Fortran 95 Compiler May 22, 2007 8 / 29

slide-18
SLIDE 18

TLM by parsetree manipulation – Example

subroutine F(x,y) double precision :: x, y y=sin(x) end subroutine

|-ROOT \-S SUBP |-SUBROUTINE | |-NAME=F | \-LIST | |-NAME=X | \-NAME=Y |-TYPE | |-DOUBLE P | |-NAME=X | \-NAME=Y |-ASGN | |-NAME=Y | \-FUNREF | |-NAME=SIN | \-NAME=X \-END |-ROOT \-S SUBP |-SUBROUTINE | |-NAME=F | \-LIST | |-NAME=X | \-NAME=Y |-USE | \-NAME=compad module |-TYPE | |-USERTYPE | | |-NAME=COMPAD TYPE | | \-COMMENT | | |-DOUBLE P | | |-ICONST=6 | | \-ICONST=0 | |-NAME=X | \-NAME=Y |-ASGN | |-NAME=Y | \-FUNREF | |-NAME=SIN | \-NAME=X \-END

Uwe Naumann, Michael Maier RWTH Aachen, Germany[1ex] Bruce Christianson, Dmitrij Gendler, Jan Riehme University of Hertfordshire, Second Order Adjoints with the NAG Fortran 95 Compiler May 22, 2007 9 / 29

slide-19
SLIDE 19

TLM by parsetree manipulation – Example

subroutine F(x,y) double precision :: x, y y=sin(x) end subroutine

|-ROOT \-S SUBP |-SUBROUTINE | |-NAME=F | \-LIST | |-NAME=X | \-NAME=Y |-TYPE | |-DOUBLE P | |-NAME=X | \-NAME=Y |-ASGN | |-NAME=Y | \-FUNREF | |-NAME=SIN | \-NAME=X \-END |-ROOT \-S SUBP |-SUBROUTINE | |-NAME=F | \-LIST | |-NAME=X | \-NAME=Y |-USE | \-NAME=compad module |-TYPE | |-USERTYPE | | |-NAME=COMPAD TYPE | | \-COMMENT | | |-DOUBLE P | | |-ICONST=6 | | \-ICONST=0 | |-NAME=X | \-NAME=Y |-ASGN | |-NAME=Y | \-FUNREF | |-NAME=SIN | \-NAME=X \-END

Uwe Naumann, Michael Maier RWTH Aachen, Germany[1ex] Bruce Christianson, Dmitrij Gendler, Jan Riehme University of Hertfordshire, Second Order Adjoints with the NAG Fortran 95 Compiler May 22, 2007 9 / 29

slide-20
SLIDE 20

TLM by parsetree manipulation – Example

subroutine F(x,y) double precision :: x, y y=sin(x) end subroutine

|-ROOT \-S SUBP |-SUBROUTINE | |-NAME=F | \-LIST | |-NAME=X | \-NAME=Y |-TYPE | |-DOUBLE P | |-NAME=X | \-NAME=Y |-ASGN | |-NAME=Y | \-FUNREF | |-NAME=SIN | \-NAME=X \-END |-ROOT \-S SUBP |-SUBROUTINE | |-NAME=F | \-LIST | |-NAME=X | \-NAME=Y |-USE | \-NAME=compad module |-TYPE | |-USERTYPE | | |-NAME=COMPAD TYPE | | \-COMMENT | | |-DOUBLE P | | |-ICONST=6 | | \-ICONST=0 | |-NAME=X | \-NAME=Y |-ASGN | |-NAME=Y | \-FUNREF | |-NAME=SIN | \-NAME=X \-END

Uwe Naumann, Michael Maier RWTH Aachen, Germany[1ex] Bruce Christianson, Dmitrij Gendler, Jan Riehme University of Hertfordshire, Second Order Adjoints with the NAG Fortran 95 Compiler May 22, 2007 9 / 29

slide-21
SLIDE 21

TLM by parsetree manipulation – Example

subroutine F(x,y) double precision :: x, y y=sin(x) end subroutine

|-ROOT \-S SUBP |-SUBROUTINE | |-NAME=F | \-LIST | |-NAME=X | \-NAME=Y |-TYPE | |-DOUBLE P | |-NAME=X | \-NAME=Y |-ASGN | |-NAME=Y | \-FUNREF | |-NAME=SIN | \-NAME=X \-END |-ROOT \-S SUBP |-SUBROUTINE | |-NAME=F | \-LIST | |-NAME=X | \-NAME=Y |-USE | \-NAME=compad module |-TYPE | |-USERTYPE | | |-NAME=COMPAD TYPE | | \-COMMENT | | |-DOUBLE P | | |-ICONST=6 | | \-ICONST=0 | |-NAME=X | \-NAME=Y |-ASGN | |-NAME=Y | \-FUNREF | |-NAME=SIN | \-NAME=X \-END

Uwe Naumann, Michael Maier RWTH Aachen, Germany[1ex] Bruce Christianson, Dmitrij Gendler, Jan Riehme University of Hertfordshire, Second Order Adjoints with the NAG Fortran 95 Compiler May 22, 2007 9 / 29

slide-22
SLIDE 22

TLM by parsetree manipulation – Example

subroutine F(x,y) double precision :: x, y y=sin(x) end subroutine

|-ROOT \-S SUBP |-SUBROUTINE | |-NAME=F | \-LIST | |-NAME=X | \-NAME=Y |-TYPE | |-DOUBLE P | |-NAME=X | \-NAME=Y |-ASGN | |-NAME=Y | \-FUNREF | |-NAME=SIN | \-NAME=X \-END |-ROOT \-S SUBP |-SUBROUTINE | |-NAME=F | \-LIST | |-NAME=X | \-NAME=Y |-USE | \-NAME=compad module |-TYPE | |-USERTYPE | | |-NAME=COMPAD TYPE | | \-COMMENT | | |-DOUBLE P | | |-ICONST=6 | | \-ICONST=0 | |-NAME=X | \-NAME=Y |-ASGN | |-NAME=Y | \-FUNREF | |-NAME=SIN | \-NAME=X \-END

Uwe Naumann, Michael Maier RWTH Aachen, Germany[1ex] Bruce Christianson, Dmitrij Gendler, Jan Riehme University of Hertfordshire, Second Order Adjoints with the NAG Fortran 95 Compiler May 22, 2007 9 / 29

slide-23
SLIDE 23

TLM — Deactivation of parameters

Parameter statement: real(kind=2), DIMENSION(2) :: X, P, Y PARAMETER( P = 3.14 ) When parser hits the parameter statement, P is already activated. So we have to find declaration node of P in parse tree detect original datatype of P create new declaration node with original data type copy all necessary attribute subtrees (dimension, allocatable, ..) remove P from the active declaration subtree keep parsetree consistent Parameter attribute: Much easier than parameter statement since attribute is detected before insertion of nodes into parsetree all variables declared are passive, so just restore original type

Uwe Naumann, Michael Maier RWTH Aachen, Germany[1ex] Bruce Christianson, Dmitrij Gendler, Jan Riehme University of Hertfordshire, Second Order Adjoints with the NAG Fortran 95 Compiler May 22, 2007 10 / 29

slide-24
SLIDE 24

TLM — Deactivation of parameters

Parameter statement: real(kind=2), DIMENSION(2) :: X, P, Y PARAMETER( P = 3.14 ) When parser hits the parameter statement, P is already activated. So we have to find declaration node of P in parse tree detect original datatype of P create new declaration node with original data type copy all necessary attribute subtrees (dimension, allocatable, ..) remove P from the active declaration subtree keep parsetree consistent Parameter attribute: Much easier than parameter statement since attribute is detected before insertion of nodes into parsetree all variables declared are passive, so just restore original type

Uwe Naumann, Michael Maier RWTH Aachen, Germany[1ex] Bruce Christianson, Dmitrij Gendler, Jan Riehme University of Hertfordshire, Second Order Adjoints with the NAG Fortran 95 Compiler May 22, 2007 10 / 29

slide-25
SLIDE 25

TLM — Deactivation of parameters

Parameter statement: real(kind=2), DIMENSION(2) :: X, P, Y PARAMETER( P = 3.14 ) When parser hits the parameter statement, P is already activated. So we have to find declaration node of P in parse tree detect original datatype of P create new declaration node with original data type copy all necessary attribute subtrees (dimension, allocatable, ..) remove P from the active declaration subtree keep parsetree consistent Parameter attribute: Much easier than parameter statement since attribute is detected before insertion of nodes into parsetree all variables declared are passive, so just restore original type

Uwe Naumann, Michael Maier RWTH Aachen, Germany[1ex] Bruce Christianson, Dmitrij Gendler, Jan Riehme University of Hertfordshire, Second Order Adjoints with the NAG Fortran 95 Compiler May 22, 2007 10 / 29

slide-26
SLIDE 26

TLM – Before / After deactivation of P

|-TYPE |-USERTYPE | |-NAME=COMPAD TYPE | \-COMMENT | |-REAL | | \-PAIR | | |-ID=KIND | | \-ICONST=2 | |-ICONST=2 | \-ICONST=2 |-DIMENSION | \-E BOUND | \-ICONST=2 |-NAME=X |-NAME=P \-NAME=Y |-TYPE |-USERTYPE | |-NAME=COMPAD TYPE | \-COMMENT | |-REAL | | \-PAIR | | |-ID=KIND | | \-ICONST=2 | |-ICONST=2 | \-ICONST=2 |-DIMENSION | \-E BOUND | \-ICONST=2 |-NAME=X \-NAME=Y |-TYPE |-REAL | \-PAIR | |-ID=KIND | \-ICONST=2 |-DIMENSION | \-E BOUND | \-ICONST=2 \-NAME=P

Uwe Naumann, Michael Maier RWTH Aachen, Germany[1ex] Bruce Christianson, Dmitrij Gendler, Jan Riehme University of Hertfordshire, Second Order Adjoints with the NAG Fortran 95 Compiler May 22, 2007 11 / 29

slide-27
SLIDE 27

TLM – Before / After deactivation of P

|-TYPE |-USERTYPE | |-NAME=COMPAD TYPE | \-COMMENT | |-REAL | | \-PAIR | | |-ID=KIND | | \-ICONST=2 | |-ICONST=2 | \-ICONST=2 |-DIMENSION | \-E BOUND | \-ICONST=2 |-NAME=X |-NAME=P \-NAME=Y |-TYPE |-USERTYPE | |-NAME=COMPAD TYPE | \-COMMENT | |-REAL | | \-PAIR | | |-ID=KIND | | \-ICONST=2 | |-ICONST=2 | \-ICONST=2 |-DIMENSION | \-E BOUND | \-ICONST=2 |-NAME=X \-NAME=Y |-TYPE |-REAL | \-PAIR | |-ID=KIND | \-ICONST=2 |-DIMENSION | \-E BOUND | \-ICONST=2 \-NAME=P

Uwe Naumann, Michael Maier RWTH Aachen, Germany[1ex] Bruce Christianson, Dmitrij Gendler, Jan Riehme University of Hertfordshire, Second Order Adjoints with the NAG Fortran 95 Compiler May 22, 2007 11 / 29

slide-28
SLIDE 28

TLM – Before / After deactivation of P

|-TYPE |-USERTYPE | |-NAME=COMPAD TYPE | \-COMMENT | |-REAL | | \-PAIR | | |-ID=KIND | | \-ICONST=2 | |-ICONST=2 | \-ICONST=2 |-DIMENSION | \-E BOUND | \-ICONST=2 |-NAME=X |-NAME=P \-NAME=Y |-TYPE |-USERTYPE | |-NAME=COMPAD TYPE | \-COMMENT | |-REAL | | \-PAIR | | |-ID=KIND | | \-ICONST=2 | |-ICONST=2 | \-ICONST=2 |-DIMENSION | \-E BOUND | \-ICONST=2 |-NAME=X \-NAME=Y |-TYPE |-REAL | \-PAIR | |-ID=KIND | \-ICONST=2 |-DIMENSION | \-E BOUND | \-ICONST=2 \-NAME=P

Uwe Naumann, Michael Maier RWTH Aachen, Germany[1ex] Bruce Christianson, Dmitrij Gendler, Jan Riehme University of Hertfordshire, Second Order Adjoints with the NAG Fortran 95 Compiler May 22, 2007 11 / 29

slide-29
SLIDE 29

TLM – Before / After deactivation of P

|-TYPE |-USERTYPE | |-NAME=COMPAD TYPE | \-COMMENT | |-REAL | | \-PAIR | | |-ID=KIND | | \-ICONST=2 | |-ICONST=2 | \-ICONST=2 |-DIMENSION | \-E BOUND | \-ICONST=2 |-NAME=X |-NAME=P \-NAME=Y |-TYPE |-USERTYPE | |-NAME=COMPAD TYPE | \-COMMENT | |-REAL | | \-PAIR | | |-ID=KIND | | \-ICONST=2 | |-ICONST=2 | \-ICONST=2 |-DIMENSION | \-E BOUND | \-ICONST=2 |-NAME=X \-NAME=Y |-TYPE |-REAL | \-PAIR | |-ID=KIND | \-ICONST=2 |-DIMENSION | \-E BOUND | \-ICONST=2 \-NAME=P

Uwe Naumann, Michael Maier RWTH Aachen, Germany[1ex] Bruce Christianson, Dmitrij Gendler, Jan Riehme University of Hertfordshire, Second Order Adjoints with the NAG Fortran 95 Compiler May 22, 2007 11 / 29

slide-30
SLIDE 30

Adjoint code generation

Adjoint code generation

Uwe Naumann, Michael Maier RWTH Aachen, Germany[1ex] Bruce Christianson, Dmitrij Gendler, Jan Riehme University of Hertfordshire, Second Order Adjoints with the NAG Fortran 95 Compiler May 22, 2007 12 / 29

slide-31
SLIDE 31

Adjoint code generation — Source trafo

Input: modified parse tree from overloading mode (forward sweep) Modify forward sweep: Inclusion of module COMPAD SYSTEM providing a value stack and a control stack Replace references to COMPAD TYPE variables by their value component Transform to single assignment code (COMPAD TYPE temporaries) Place a value - push before every overwriting of a value component Place a control - push after of every logical decision Push array indices to control stack Create incremental reverse sweep from forward sweep Process assignments of forward sweep in reverse order Create Adjoint code for every assignment Place value - pops and control - pops where necessary Reset adjoint of left hand side

Uwe Naumann, Michael Maier RWTH Aachen, Germany[1ex] Bruce Christianson, Dmitrij Gendler, Jan Riehme University of Hertfordshire, Second Order Adjoints with the NAG Fortran 95 Compiler May 22, 2007 13 / 29

slide-32
SLIDE 32

Adjoint code generation — Source trafo

Input: modified parse tree from overloading mode (forward sweep) Modify forward sweep: Inclusion of module COMPAD SYSTEM providing a value stack and a control stack Replace references to COMPAD TYPE variables by their value component Transform to single assignment code (COMPAD TYPE temporaries) Place a value - push before every overwriting of a value component Place a control - push after of every logical decision Push array indices to control stack Create incremental reverse sweep from forward sweep Process assignments of forward sweep in reverse order Create Adjoint code for every assignment Place value - pops and control - pops where necessary Reset adjoint of left hand side

Uwe Naumann, Michael Maier RWTH Aachen, Germany[1ex] Bruce Christianson, Dmitrij Gendler, Jan Riehme University of Hertfordshire, Second Order Adjoints with the NAG Fortran 95 Compiler May 22, 2007 13 / 29

slide-33
SLIDE 33

Adjoint code generation — Source trafo

Input: modified parse tree from overloading mode (forward sweep) Modify forward sweep: Inclusion of module COMPAD SYSTEM providing a value stack and a control stack Replace references to COMPAD TYPE variables by their value component Transform to single assignment code (COMPAD TYPE temporaries) Place a value - push before every overwriting of a value component Place a control - push after of every logical decision Push array indices to control stack Create incremental reverse sweep from forward sweep Process assignments of forward sweep in reverse order Create Adjoint code for every assignment Place value - pops and control - pops where necessary Reset adjoint of left hand side

Uwe Naumann, Michael Maier RWTH Aachen, Germany[1ex] Bruce Christianson, Dmitrij Gendler, Jan Riehme University of Hertfordshire, Second Order Adjoints with the NAG Fortran 95 Compiler May 22, 2007 13 / 29

slide-34
SLIDE 34

Adjoint code generation – Further information

Note that Pure mode: operations are intrinsic operators and functions again Derived mode: overloaded operations or intrinsics for VALUE TYPE and DERIV TYPE

More information about adjoint code generation

Listen to Michael Maier in this session

Uwe Naumann, Michael Maier RWTH Aachen, Germany[1ex] Bruce Christianson, Dmitrij Gendler, Jan Riehme University of Hertfordshire, Second Order Adjoints with the NAG Fortran 95 Compiler May 22, 2007 14 / 29

slide-35
SLIDE 35

Adjoint code generation – Further information

Note that Pure mode: operations are intrinsic operators and functions again Derived mode: overloaded operations or intrinsics for VALUE TYPE and DERIV TYPE

More information about adjoint code generation

Listen to Michael Maier in this session

Uwe Naumann, Michael Maier RWTH Aachen, Germany[1ex] Bruce Christianson, Dmitrij Gendler, Jan Riehme University of Hertfordshire, Second Order Adjoints with the NAG Fortran 95 Compiler May 22, 2007 14 / 29

slide-36
SLIDE 36

Second Order Adjoints

Second order adjoints in a Truncated Newton Method

Uwe Naumann, Michael Maier RWTH Aachen, Germany[1ex] Bruce Christianson, Dmitrij Gendler, Jan Riehme University of Hertfordshire, Second Order Adjoints with the NAG Fortran 95 Compiler May 22, 2007 15 / 29

slide-37
SLIDE 37

Truncated Newton — Task and Properties

Task: minx∈Rn f (x) where f : Rn → R f twice continously differentiable f has bounded level sets {x : f (x) ≤ f (x0)} ∀x0 ∈ Rn Some properties of the Truncated Newton approach avoid saddlepoints (detection of directions with negative curvature) requires Hessian vector products only, not the full Hessian automatic exploitation of symmetries in f in theory for exact convergence:

# inner iterations = # distinct eigenvalues of f which is often much smaller than dimension of x

Uwe Naumann, Michael Maier RWTH Aachen, Germany[1ex] Bruce Christianson, Dmitrij Gendler, Jan Riehme University of Hertfordshire, Second Order Adjoints with the NAG Fortran 95 Compiler May 22, 2007 16 / 29

slide-38
SLIDE 38

Truncated Newton — Task and Properties

Task: minx∈Rn f (x) where f : Rn → R f twice continously differentiable f has bounded level sets {x : f (x) ≤ f (x0)} ∀x0 ∈ Rn Some properties of the Truncated Newton approach avoid saddlepoints (detection of directions with negative curvature) requires Hessian vector products only, not the full Hessian automatic exploitation of symmetries in f in theory for exact convergence:

# inner iterations = # distinct eigenvalues of f which is often much smaller than dimension of x

Uwe Naumann, Michael Maier RWTH Aachen, Germany[1ex] Bruce Christianson, Dmitrij Gendler, Jan Riehme University of Hertfordshire, Second Order Adjoints with the NAG Fortran 95 Compiler May 22, 2007 16 / 29

slide-39
SLIDE 39

Truncated Newton — Algorithm outline

Major iteration Input:: start point x0 ∈ Rn Test for convergence, Stop if convergence Minor iteration to find new direction pk

p(0)

k

= −g(xk) requires gradient at old xk Solve H(xk)p(i)

k

= −g(xk) until convergence

Generates a sequence p(i)

k

Repeated Hessian - vector products required at fixed xk

Calculate step length λk requires gradient ∇f (xk + λkpk) xk+1 = xk + λkpk

Uwe Naumann, Michael Maier RWTH Aachen, Germany[1ex] Bruce Christianson, Dmitrij Gendler, Jan Riehme University of Hertfordshire, Second Order Adjoints with the NAG Fortran 95 Compiler May 22, 2007 17 / 29

slide-40
SLIDE 40

Hybrid approach for second order adjoints

Hybrid means here Source code transformation to get the adjoint code Thus we get x and ¯ x But adjoint code deals with variables of derived type TLS TYPE COMPAD TYPE provides overloaded operators for TLS TYPE This gives a forward mode differentiation to get ˙ x and ˙ ¯ x In fact we get a mix of source transformation and Overloading How does the COMPAD TYPE looks like? type compad type type(tls type) :: val ! value x type(tls type) :: drv ! adjoint ¯ x end type compad type type tls type double precision :: val=0 ! tangent ˙ x double precision :: drv=0 ! second order ˙ ¯ x end type tls type

Uwe Naumann, Michael Maier RWTH Aachen, Germany[1ex] Bruce Christianson, Dmitrij Gendler, Jan Riehme University of Hertfordshire, Second Order Adjoints with the NAG Fortran 95 Compiler May 22, 2007 18 / 29

slide-41
SLIDE 41

Hybrid approach for second order adjoints

Hybrid means here Source code transformation to get the adjoint code Thus we get x and ¯ x But adjoint code deals with variables of derived type TLS TYPE COMPAD TYPE provides overloaded operators for TLS TYPE This gives a forward mode differentiation to get ˙ x and ˙ ¯ x In fact we get a mix of source transformation and Overloading How does the COMPAD TYPE looks like? type compad type type(tls type) :: val ! value x type(tls type) :: drv ! adjoint ¯ x end type compad type type tls type double precision :: val=0 ! tangent ˙ x double precision :: drv=0 ! second order ˙ ¯ x end type tls type

Uwe Naumann, Michael Maier RWTH Aachen, Germany[1ex] Bruce Christianson, Dmitrij Gendler, Jan Riehme University of Hertfordshire, Second Order Adjoints with the NAG Fortran 95 Compiler May 22, 2007 18 / 29

slide-42
SLIDE 42

Truncated Newton — Implementation

  • ptnhb.for , a Fortran 77 code from the optima package

C ********************************************************** C * HATFIELD POLYTECHNIC * C * NUMERICAL OPTIMISATION CENTRE * C * THE OPTIMA PACKAGE * C * DOUBLE PRECISION SUBROUTINE: OPTNHP.FOR * C * VERSION 1, JANUARY 1989 * C ********************************************************** C C C PURPOSE: To find a local minimum of a function of several C variables using truncated newton method. C C AUTHOR: Z.A. MAANY, NUMERICAL OPTIMISATION CENTRE C C USAGE:CALL OPTNHP (N,X,FF,G,MAXIT,EPS,IPRINT,W,NW)

Subroutine OPTNHB assume existence of 3 evaluation routines SUBROUTINE CALFUN(X,N,FF,INF) function SUBROUTINE CALGRD(X,N,FF,G,INF) gradient SUBROUTINE CALHP(X,G,N,P,HP,INF) Hessian-vector

Uwe Naumann, Michael Maier RWTH Aachen, Germany[1ex] Bruce Christianson, Dmitrij Gendler, Jan Riehme University of Hertfordshire, Second Order Adjoints with the NAG Fortran 95 Compiler May 22, 2007 19 / 29

slide-43
SLIDE 43

Truncated Newton — Implementation

  • ptnhb.for , a Fortran 77 code from the optima package

C ********************************************************** C * HATFIELD POLYTECHNIC * C * NUMERICAL OPTIMISATION CENTRE * C * THE OPTIMA PACKAGE * C * DOUBLE PRECISION SUBROUTINE: OPTNHP.FOR * C * VERSION 1, JANUARY 1989 * C ********************************************************** C C C PURPOSE: To find a local minimum of a function of several C variables using truncated newton method. C C AUTHOR: Z.A. MAANY, NUMERICAL OPTIMISATION CENTRE C C USAGE:CALL OPTNHP (N,X,FF,G,MAXIT,EPS,IPRINT,W,NW)

Subroutine OPTNHB assume existence of 3 evaluation routines SUBROUTINE CALFUN(X,N,FF,INF) function SUBROUTINE CALGRD(X,N,FF,G,INF) gradient SUBROUTINE CALHP(X,G,N,P,HP,INF) Hessian-vector

Uwe Naumann, Michael Maier RWTH Aachen, Germany[1ex] Bruce Christianson, Dmitrij Gendler, Jan Riehme University of Hertfordshire, Second Order Adjoints with the NAG Fortran 95 Compiler May 22, 2007 19 / 29

slide-44
SLIDE 44

Truncated Newton — Test run on Rosenbrock

Rosenbrock problem: f (x, y) = (1 − x)2 + 100(y − x2)2

2 4 6 8 10 12 0.2 0.4 0.6 0.8 1 1.2 1.4

  • 1
  • 0.5

0.5 1 1.5 2 2.5 3 2 4 6 8 10 12 ((1-x)**2 + 100*(y-x*x)**2)/100 Uwe Naumann, Michael Maier RWTH Aachen, Germany[1ex] Bruce Christianson, Dmitrij Gendler, Jan Riehme University of Hertfordshire, Second Order Adjoints with the NAG Fortran 95 Compiler May 22, 2007 20 / 29

slide-45
SLIDE 45

TN — Results

main x 3.0000000000000000 4.0000000000000000 SUBROUTINE OPTNHP STARTED ************************* ’CONJUGATE GRADIENT TRUNCATED NEWTON’ CONVERGE IF EUCLIDEAN NORM OF GRAD F < 1.0000000D-08 . . . Second order adjoints NEW POINT 1.0000000D-00 1.0000000D-00 FUNCTION = 1.7998949D-18 GRADIENT

  • 1.8523072D-10 -1.2475354D-09

29 FUNCTION CALLS 31 MINOR ITERATIONS 20 MAJOR ITERATIONS CONVERGENCE Divided Differences (ǫ = 10−7) NEW POINT 1.0000089D+00 1.0000178D+00 FUNCTION = 7.9281058D-11 GRADIENT 2.1743654D-05 2.8036572D-05 33884 FUNCTION CALLS 8757 MINOR ITERATIONS 6017 MAJOR ITERATIONS NO FURTHER FUNCTION DECREASE OBTAINED

Uwe Naumann, Michael Maier RWTH Aachen, Germany[1ex] Bruce Christianson, Dmitrij Gendler, Jan Riehme University of Hertfordshire, Second Order Adjoints with the NAG Fortran 95 Compiler May 22, 2007 21 / 29

slide-46
SLIDE 46

TN — Results

main x 3.0000000000000000 4.0000000000000000 SUBROUTINE OPTNHP STARTED ************************* ’CONJUGATE GRADIENT TRUNCATED NEWTON’ CONVERGE IF EUCLIDEAN NORM OF GRAD F < 1.0000000D-08 . . . Second order adjoints NEW POINT 1.0000000D-00 1.0000000D-00 FUNCTION = 1.7998949D-18 GRADIENT

  • 1.8523072D-10 -1.2475354D-09

29 FUNCTION CALLS 31 MINOR ITERATIONS 20 MAJOR ITERATIONS CONVERGENCE Divided Differences (ǫ = 10−7) NEW POINT 1.0000089D+00 1.0000178D+00 FUNCTION = 7.9281058D-11 GRADIENT 2.1743654D-05 2.8036572D-05 33884 FUNCTION CALLS 8757 MINOR ITERATIONS 6017 MAJOR ITERATIONS NO FURTHER FUNCTION DECREASE OBTAINED

Uwe Naumann, Michael Maier RWTH Aachen, Germany[1ex] Bruce Christianson, Dmitrij Gendler, Jan Riehme University of Hertfordshire, Second Order Adjoints with the NAG Fortran 95 Compiler May 22, 2007 21 / 29

slide-47
SLIDE 47

TN — Results

main x 3.0000000000000000 4.0000000000000000 SUBROUTINE OPTNHP STARTED ************************* ’CONJUGATE GRADIENT TRUNCATED NEWTON’ CONVERGE IF EUCLIDEAN NORM OF GRAD F < 1.0000000D-08 . . . Second order adjoints NEW POINT 1.0000000D-00 1.0000000D-00 FUNCTION = 1.7998949D-18 GRADIENT

  • 1.8523072D-10 -1.2475354D-09

29 FUNCTION CALLS 31 MINOR ITERATIONS 20 MAJOR ITERATIONS CONVERGENCE Divided Differences (ǫ = 10−7) NEW POINT 1.0000089D+00 1.0000178D+00 FUNCTION = 7.9281058D-11 GRADIENT 2.1743654D-05 2.8036572D-05 33884 FUNCTION CALLS 8757 MINOR ITERATIONS 6017 MAJOR ITERATIONS NO FURTHER FUNCTION DECREASE OBTAINED

Uwe Naumann, Michael Maier RWTH Aachen, Germany[1ex] Bruce Christianson, Dmitrij Gendler, Jan Riehme University of Hertfordshire, Second Order Adjoints with the NAG Fortran 95 Compiler May 22, 2007 21 / 29

slide-48
SLIDE 48

TN — Main program

PROGRAM tn ros ovl IMPLICIT NONE INTEGER, PARAMETER :: N = 2 DOUBLE PRECISION :: X(N),G(N), FF INTEGER, PARAMETER :: NW = 160 DOUBLE PRECISION :: W(NW), EPS INTEGER :: IOC, IPRINT, IMAX COMMON/DEVICE/IOC IOC=6 WRITE(IOC,101) READ(5,111)IPRINT,IMAX,EPS X(1) = 3 ; X(2) = 4 print *,’main x ’, x CALL OPTNHP(N,X,FF,G,IMAX,EPS,IPRINT,W,NW) STOP

Uwe Naumann, Michael Maier RWTH Aachen, Germany[1ex] Bruce Christianson, Dmitrij Gendler, Jan Riehme University of Hertfordshire, Second Order Adjoints with the NAG Fortran 95 Compiler May 22, 2007 22 / 29

slide-49
SLIDE 49

TN — Subprogram rosenbrock sub.f

SUBROUTINE rosenbrock ( a, res) implicit none DOUBLE PRECISION :: a(2) , res INTEGER :: i = 1 ! Rosenbrock Function. ! f(x,y) = (1-x)**2 + 100(y-x**2)**2 DOUBLE PRECISION :: t1, t2, t3 t1 = 1 - a(i) t2 = t1*t1 t3 = a(2) - a(i)*a(i) res = t2 + 100 * t3*t3 END SUBROUTINE

Uwe Naumann, Michael Maier RWTH Aachen, Germany[1ex] Bruce Christianson, Dmitrij Gendler, Jan Riehme University of Hertfordshire, Second Order Adjoints with the NAG Fortran 95 Compiler May 22, 2007 23 / 29

slide-50
SLIDE 50

TN — Subprogram rosenbrock sub.f II

Since compiler does not create new names for differentiated routines at the moment . . .

#ifdef AD ADJOINT SUBROUTINE rosenbrock adj (a, res) USE COMPAD SYSTEM #elif AD OVL SUBROUTINE rosenbrock ovl (a, res) #else SUBROUTINE rosenbrock ( a, res) #endif

But he will in the future ...

Uwe Naumann, Michael Maier RWTH Aachen, Germany[1ex] Bruce Christianson, Dmitrij Gendler, Jan Riehme University of Hertfordshire, Second Order Adjoints with the NAG Fortran 95 Compiler May 22, 2007 24 / 29

slide-51
SLIDE 51

TN — Subroutine CALFUN

SUBROUTINE CALFUN(X,N,FF,INF) IMPLICIT NONE INTEGER :: N, INF DOUBLE PRECISION :: X(N), FF ! INF != 0 signals infeasibility INF=0 call rosenbrock(x,FF) RETURN END

Uwe Naumann, Michael Maier RWTH Aachen, Germany[1ex] Bruce Christianson, Dmitrij Gendler, Jan Riehme University of Hertfordshire, Second Order Adjoints with the NAG Fortran 95 Compiler May 22, 2007 25 / 29

slide-52
SLIDE 52

TN — Subroutine CALHP

SUBROUTINE CALHP(X,G,N,P,HP,INF) use compad module use compad system IMPLICIT NONE INTEGER :: N, INF, I DOUBLE PRECISION :: X(N),G(N),P(N),HP(N) TYPE(COMPAD TYPE) :: CTX(2), CTF INF=0 DO i=1,2 ctx(i)%val%val = X(i) ! xi value ctx(i)%val%drv = P(i) ! ˙ xi ctx(i)%drv%val = 0.D0 ! ¯ xi ctx(i)%drv%drv = 0.D0 ! ˙ ¯ xi ENDDO ctf%val%val = 0.D0 ! f value ctf%val%drv = 0.D0 ! ˙ f ctf%drv%val = 1.D0 ! ¯ f ctf%drv%drv = 0.D0 ! ˙ ¯ f call rosenbrock adj(ctx,ctf) HP(1) = ctx(1)%drv%drv HP(2) = ctx(2)%drv%drv

Uwe Naumann, Michael Maier RWTH Aachen, Germany[1ex] Bruce Christianson, Dmitrij Gendler, Jan Riehme University of Hertfordshire, Second Order Adjoints with the NAG Fortran 95 Compiler May 22, 2007 26 / 29

slide-53
SLIDE 53

TN — Subroutine CALHP

SUBROUTINE CALHP(X,G,N,P,HP,INF) use compad module use compad system IMPLICIT NONE INTEGER :: N, INF, I DOUBLE PRECISION :: X(N),G(N),P(N),HP(N) TYPE(COMPAD TYPE) :: CTX(2), CTF INF=0 DO i=1,2 ctx(i)%val%val = X(i) ! xi value ctx(i)%val%drv = P(i) ! ˙ xi ctx(i)%drv%val = 0.D0 ! ¯ xi ctx(i)%drv%drv = 0.D0 ! ˙ ¯ xi ENDDO ctf%val%val = 0.D0 ! f value ctf%val%drv = 0.D0 ! ˙ f ctf%drv%val = 1.D0 ! ¯ f ctf%drv%drv = 0.D0 ! ˙ ¯ f call rosenbrock adj(ctx,ctf) HP(1) = ctx(1)%drv%drv HP(2) = ctx(2)%drv%drv

Uwe Naumann, Michael Maier RWTH Aachen, Germany[1ex] Bruce Christianson, Dmitrij Gendler, Jan Riehme University of Hertfordshire, Second Order Adjoints with the NAG Fortran 95 Compiler May 22, 2007 26 / 29

slide-54
SLIDE 54

TN — Subroutine CALHP

SUBROUTINE CALHP(X,G,N,P,HP,INF) use compad module use compad system IMPLICIT NONE INTEGER :: N, INF, I DOUBLE PRECISION :: X(N),G(N),P(N),HP(N) TYPE(COMPAD TYPE) :: CTX(2), CTF INF=0 DO i=1,2 ctx(i)%val%val = X(i) ! xi value ctx(i)%val%drv = P(i) ! ˙ xi ctx(i)%drv%val = 0.D0 ! ¯ xi ctx(i)%drv%drv = 0.D0 ! ˙ ¯ xi ENDDO ctf%val%val = 0.D0 ! f value ctf%val%drv = 0.D0 ! ˙ f ctf%drv%val = 1.D0 ! ¯ f ctf%drv%drv = 0.D0 ! ˙ ¯ f call rosenbrock adj(ctx,ctf) HP(1) = ctx(1)%drv%drv HP(2) = ctx(2)%drv%drv

Uwe Naumann, Michael Maier RWTH Aachen, Germany[1ex] Bruce Christianson, Dmitrij Gendler, Jan Riehme University of Hertfordshire, Second Order Adjoints with the NAG Fortran 95 Compiler May 22, 2007 26 / 29

slide-55
SLIDE 55

TN — Subroutine CALHP

SUBROUTINE CALHP(X,G,N,P,HP,INF) use compad module use compad system IMPLICIT NONE INTEGER :: N, INF, I DOUBLE PRECISION :: X(N),G(N),P(N),HP(N) TYPE(COMPAD TYPE) :: CTX(2), CTF INF=0 DO i=1,2 ctx(i)%val%val = X(i) ! xi value ctx(i)%val%drv = P(i) ! ˙ xi ctx(i)%drv%val = 0.D0 ! ¯ xi ctx(i)%drv%drv = 0.D0 ! ˙ ¯ xi ENDDO ctf%val%val = 0.D0 ! f value ctf%val%drv = 0.D0 ! ˙ f ctf%drv%val = 1.D0 ! ¯ f ctf%drv%drv = 0.D0 ! ˙ ¯ f call rosenbrock adj(ctx,ctf) HP(1) = ctx(1)%drv%drv HP(2) = ctx(2)%drv%drv

Uwe Naumann, Michael Maier RWTH Aachen, Germany[1ex] Bruce Christianson, Dmitrij Gendler, Jan Riehme University of Hertfordshire, Second Order Adjoints with the NAG Fortran 95 Compiler May 22, 2007 26 / 29

slide-56
SLIDE 56

TN — Subroutine CALGRAD

Identical with CALHP() beside different initialisation

DO i=1,2 ctx(i)%val%val = X(i) ! xi value ctx(i)%val%drv = 0.D0 ! ˙ xi ctx(i)%drv%val = 0.D0 ! ¯ xi ctx(i)%drv%drv = 0.D0 ! ˙ ¯ xi ENDDO ctf%val%val = 0.D0 ! f value ctf%val%drv = 0.D0 ! ˙ f ctf%drv%val = 1.D0 ! ¯ f ctf%drv%drv = 0.D0 ! ˙ ¯ f call rosenbrock adj(ctx,ctf) G(1) = ctx(1)%drv%val G(2) = ctx(2)%drv%val

Note that second order adjoints are computed here too. To avoid that a second variant of COMPAD TYPE would be required (major effort).

Uwe Naumann, Michael Maier RWTH Aachen, Germany[1ex] Bruce Christianson, Dmitrij Gendler, Jan Riehme University of Hertfordshire, Second Order Adjoints with the NAG Fortran 95 Compiler May 22, 2007 27 / 29

slide-57
SLIDE 57

TN — Subroutine CALGRAD

Identical with CALHP() beside different initialisation

DO i=1,2 ctx(i)%val%val = X(i) ! xi value ctx(i)%val%drv = 0.D0 ! ˙ xi ctx(i)%drv%val = 0.D0 ! ¯ xi ctx(i)%drv%drv = 0.D0 ! ˙ ¯ xi ENDDO ctf%val%val = 0.D0 ! f value ctf%val%drv = 0.D0 ! ˙ f ctf%drv%val = 1.D0 ! ¯ f ctf%drv%drv = 0.D0 ! ˙ ¯ f call rosenbrock adj(ctx,ctf) G(1) = ctx(1)%drv%val G(2) = ctx(2)%drv%val

Note that second order adjoints are computed here too. To avoid that a second variant of COMPAD TYPE would be required (major effort).

Uwe Naumann, Michael Maier RWTH Aachen, Germany[1ex] Bruce Christianson, Dmitrij Gendler, Jan Riehme University of Hertfordshire, Second Order Adjoints with the NAG Fortran 95 Compiler May 22, 2007 27 / 29

slide-58
SLIDE 58

TN — Subroutine CALGRAD

Identical with CALHP() beside different initialisation

DO i=1,2 ctx(i)%val%val = X(i) ! xi value ctx(i)%val%drv = 0.D0 ! ˙ xi ctx(i)%drv%val = 0.D0 ! ¯ xi ctx(i)%drv%drv = 0.D0 ! ˙ ¯ xi ENDDO ctf%val%val = 0.D0 ! f value ctf%val%drv = 0.D0 ! ˙ f ctf%drv%val = 1.D0 ! ¯ f ctf%drv%drv = 0.D0 ! ˙ ¯ f call rosenbrock adj(ctx,ctf) G(1) = ctx(1)%drv%val G(2) = ctx(2)%drv%val

Note that second order adjoints are computed here too. To avoid that a second variant of COMPAD TYPE would be required (major effort).

Uwe Naumann, Michael Maier RWTH Aachen, Germany[1ex] Bruce Christianson, Dmitrij Gendler, Jan Riehme University of Hertfordshire, Second Order Adjoints with the NAG Fortran 95 Compiler May 22, 2007 27 / 29

slide-59
SLIDE 59

Outlook

Main targets: Robust adjoint compiler Adjoints of MINPACK-II problems Adjoints of CNS – code Checkpointing of loops Secondary targets: Coupling of adjoint compiler with ADAC Coupling of adjoint compiler with ADOLC

Uwe Naumann, Michael Maier RWTH Aachen, Germany[1ex] Bruce Christianson, Dmitrij Gendler, Jan Riehme University of Hertfordshire, Second Order Adjoints with the NAG Fortran 95 Compiler May 22, 2007 28 / 29

slide-60
SLIDE 60

Outlook

Main targets: Robust adjoint compiler Adjoints of MINPACK-II problems Adjoints of CNS – code Checkpointing of loops Secondary targets: Coupling of adjoint compiler with ADAC Coupling of adjoint compiler with ADOLC

Uwe Naumann, Michael Maier RWTH Aachen, Germany[1ex] Bruce Christianson, Dmitrij Gendler, Jan Riehme University of Hertfordshire, Second Order Adjoints with the NAG Fortran 95 Compiler May 22, 2007 28 / 29

slide-61
SLIDE 61

Thank you!

Uwe Naumann, Michael Maier RWTH Aachen, Germany[1ex] Bruce Christianson, Dmitrij Gendler, Jan Riehme University of Hertfordshire, Second Order Adjoints with the NAG Fortran 95 Compiler May 22, 2007 29 / 29