Getting Started (A Toy Problem)

What do we mean when we say that a tangent-linear and an adjoint model are generated for a given numerical simulation program?

Consider the following toy Fortran 90 code (also contained in the downloadable test setup)

        subroutine head(x,y) 
          double precision,intent(in) :: x
          double precision,intent(out) :: y
          y=sin(x*x)
        end subroutine

We need to indicate to the system what the independent and the dependent variables are with pragmas:

        subroutine head(x,y) 
          double precision,intent(in) :: x
          double precision,intent(out) :: y
c$openad INDEPENDENT(x)
          y=sin(x*x)
c$openad DEPENDENT(y)
        end subroutine

Tangent Linear Code

The tangent linear code computes

• all intermediate and the final function values;
  (A single intermediate value is computed when evaluating head.)
• partial derivatives of the arithmetic operations and intrinsic function calls;
  (We compute the two partials of the product and the single partial of the sine intrinsic. The original code augmented with statements for the computation of all local partials is referred to as linearized.)
• the Jacobian by preaccumulation using elimination techniques;
  (Since the program consists of a single assignment the Jacobian of the assignment that can be obtained by preaccumulation is equivalent to the Jacobian of the program.)
• the product of the Jacobian with a direction in the input space;
  (This is what a tangent-linear model is all about.)

      SUBROUTINE head(X, Y)
      use w2f__types
      use active_module
      IMPLICIT NONE
      REAL(w2f__8) OpenAD_Symbol_0
      REAL(w2f__8) OpenAD_Symbol_1
      REAL(w2f__8) OpenAD_Symbol_2
      REAL(w2f__8) OpenAD_Symbol_3
      REAL(w2f__8) OpenAD_Symbol_5
      type(active) :: X
      INTENT(IN)  X
      type(active) :: Y
      INTENT(OUT)  Y
      OpenAD_Symbol_0 = (X%v*X%v)
      Y%v = SIN(OpenAD_Symbol_0)
      OpenAD_Symbol_2 = X%v
      OpenAD_Symbol_3 = X%v
      OpenAD_Symbol_1 = COS(OpenAD_Symbol_0)
      OpenAD_Symbol_5 = ((OpenAD_Symbol_3 + OpenAD_Symbol_2) * OpenAD_Symbol_1)
      CALL sax(OpenAD_Symbol_5,X,Y)
      RETURN
      END SUBROUTINE

The code has been modified for the purpose of presentation. You may want to look at the complete tangent-linear code that is generated by OpenAD. A new user-defined data type (active) is introduced to hold for a given scalar intermediate value both its value (%v) and the value of its directional derivative (%d) with respect to the independent variables (X) at the current point. The extra variables generated by xaifBooster all have the OpenAD_Symbol_ prefix. The runtime support library implements the subroutine sax(a,x,y) as: y%d=a*x%d. Note, that this example shows some extraneous assignments of the X%v to generated variables. They can be removed by a subsequent optimizing compilation of this code. The XAIF representation of head can be found here.

Adjoint Code

The main emphasis of our work is on the generation of robust and efficient adjoint code. In many applications the code will consist of multiple subroutines and we support several reversal schemes at the subroutine level. While not really relevant for a single subroutine we need to complete the setup by indicating the reversal scheme as pragma such as c$openad XXX Template ad_template.f  to complete the input code. All relevant local partial derivatives of the arithmetic operators and intrinsic functions are stored on a tape (double_tape) during a single forward evaluation of the linearized code.

      OpenAD_Symbol_0 = (X%v*X%v)
      Y%v = SIN(OpenAD_Symbol_0)
      OpenAD_Symbol_2 = X%v
      OpenAD_Symbol_3 = X%v
      OpenAD_Symbol_1 = COS(OpenAD_Symbol_0)
      OpenAD_Symbol_5 = ((OpenAD_Symbol_3 + OpenAD_Symbol_2) * OpenAD_Symbol_1)
      double_tape(double_tape_pointer) = OpenAD_Symbol_5
      double_tape_pointer = double_tape_pointer+1

The reverse sweep restores the partial derivatives in reverse order and uses them to compute local transposed Jacobian-vector products.

      double_tape_pointer = double_tape_pointer-1
      OpenAD_Symbol_7 = double_tape(double_tape_pointer)
      X%d = X%d+Y%d*OpenAD_Symbol_7
      Y%d = 0.0d0

The adjoint of y%d=a*x%d (sax(a,x,y)) is x%d=x%d+a*y%d followed by nullifying y%d. Adjoints of all intermediate variable with respect to the dependent variables (y) are stored in the %d component. The complete adjoint code contains all sections for a variety of reversal schemes. A report that describes the various approaches to the reversal of the call graph is under development.