!********************************************************************************
!* Automatic Differentiation (AD) example program, Version 1 (2012),
!*
!* written by Dr. Katate Masatsuka (info[at]cfdbooks.com),
!*
!* the author of useful CFD books, "I do like CFD" (http://www.cfdbooks.com).
!*
!* This file is written and made available for an educational purpose.
!* It is not complete and may not be the best implementation.
!* It only attempts to let beginners understand how derivatives can be
!* computed by the automatic differentiation.
!*
!* If it did help beginners understand the idea of automatic differentiation,
!* it would make the author feel happy.
!*
!* Suppose you have a function that computes g(x) for a given value of x,
!* where g and x are assumed to be scalars. For example, g(x) can be
!*
!* function g(x)
!* real :: x, g
!* g = x**2 + 5.0*x + sqrt(x)
!* end function
!*
!* To compute the derivative, dg/dx, we rewrite this function in terms of
!* "derivative data type", instead of "real":
!*
!* function g(x)
!* type(derivative_data_type) :: x, g
!* g = x**2 + 5.0*x + ddt_sqrt(x)
!* end function
!*
!* Then, this computes the function value as well as the derivative!
!* Trick is that the derivative_data_type carries the function value and derivative
!* as its components (x%f, x%df, and g%f, g%df) and the derivative components are
!* computed simultaneously at every step of the computations because of the newly
!* defined operators (**, *, ddt_sqrt, etc, are overloaded for derivative_data_type).
!* For example, multiplication operator * is re-defined such that it applies to
!* both the function value and the derivative, e.g., 5.0*x is applied to the
!* function value (%f) as well as the derivative component (%df):
!*
!* 5.0 * x%f
!* 5.0 * x%df
!*
!* This way, the function looks very much like just a function, but it is
!* actually a function that computes the function value and the derivative
!* simulntanesously.
!*
!* In CFD, AD can be useful for computing flux Jacobians for implicit schemes.
!*
!* This file may be updated in future.
!*
!* Katate Masatsuka, November 2012. http://www.cfdbooks.com
!********************************************************************************
! This module contains the definition of the derivative data type (ddt), and the
! associated operators for ddt. This one is specifically prepared for vairables
! carrying 3 derivatives. To extend it to accomodate more derivatives, generate
! module_ddt4.f90, module_ddt5.f90, module_ddt6.f90, and so on, by modifying the
! size of the derivative array (i.e., %df(:) ). See the file for details.
include 'module_ddt3.f90'
!********************************************************************************
!* This program runs a test problem for automatic differentiation.
!*
!* --- Compute dg/dx, where g = (g1(x),g2(x),g3(x)) and x=(x1,x2,x3) ---
!*
!* In the test problem, the AD result is compared with the analytical formula.
!* The difference will be printed, which must be zero.
!*
!* Compile and run this program (e.g., gfortran -O2 ad_driver.f90).
!* See screen output for details.
!*
!* You may want to make a loop over, say, 10000 evaluations, and measure and
!* compare the CPU time, AD versus exact derivative formula.
!*
!* Katate Masatsuka, November 2012. http://www.cfdbooks.com
!********************************************************************************
program test_automatic_differentiation
use derivative_data_df3 ! <- To use the derivative data type with 3 derivatives.
integer , parameter :: p2 = selected_real_kind(15)
real(p2), parameter :: zero = 0.0_p2
real(p2) , dimension(3) :: x_real
real(p2) , dimension(3,3) :: exact_derivatives
type(derivative_data_type_df3), dimension(3) :: g, x
integer :: j
! Note: g and x are derivative_data containing a function and its derivatives. E.g.,
! g%f is the function value, and g%df(1:3) is the array that store derivatives.
! We will store the 3x3 matrix dg/dx as g(i)%df(j), i=1,3, j=1,3.
write(*,*) "*******************************************************************"
write(*,*) " Test problem: Derivatives of a simple vector function"
write(*,*) "*******************************************************************"
write(*,*)
write(*,*) " In this program, we compute the gradient of a vector function, dg/dx,"
write(*,*) " where g=(g1,g2,g3), x=(x1,x2,x3) - Note: dg/dx is a matrix."
write(*,*) " We compute dg/dx at x=(2, 5, 7)."
write(*,*) " "
write(*,*) " x is defined as an array of derivative_data_type variables, x(1:3)"
write(*,*) " g is defined as an array of derivative_data_type variables, g(1:3)"
write(*,*) " "
! Some random function values at which dg/dx will be computed.
! Set any number you like.
x_real(1) = 2.0_p2
x_real(2) = 5.0_p2
x_real(3) = 7.0_p2
! Initialization of x
write(*,*)
write(*,*) "Initializing x..."
! Set function values by x_real; derivatives are zero at this point.
x = x_real ! This will set the function values, i.e., x%f = x_real (and x%df = 0.0)
write(*,*)
write(*,'(a,es25.15,a,3es25.15)') " x(1)%f=",x(1)%f, " x(1)%df(1), x(1)%df(2), x(1)%df(3) =",x(1)%df
write(*,'(a,es25.15,a,3es25.15)') " x(2)%f=",x(2)%f, " x(2)%df(1), x(2)%df(2), x(2)%df(3) =",x(2)%df
write(*,'(a,es25.15,a,3es25.15)') " x(3)%f=",x(3)%f, " x(3)%df(1), x(3)%df(2), x(3)%df(3) =",x(3)%df
write(*,*)
write(*,*) "Initializing g..."
g = zero ! This will set g%f = 0.0 and g%df = 0.0
write(*,*)
write(*,'(a,es25.15,a,3es25.15)') " g(1)%f=",g(1)%f, " g(1)%df(1), g(1)%df(2), g(1)%df(3) =",g(1)%df
write(*,'(a,es25.15,a,3es25.15)') " g(2)%f=",g(2)%f, " g(2)%df(1), g(2)%df(2), g(2)%df(3) =",g(2)%df
write(*,'(a,es25.15,a,3es25.15)') " g(3)%f=",g(3)%f, " g(3)%df(1), g(3)%df(2), g(3)%df(3) =",g(3)%df
write(*,*)
write(*,*) "Initialization done."
! Set correct derivatives for x, which is the identity matrix since the derivative is dx/dx.
! That is, we set x(i)%df(i) = 1.0 for i=1,2,3.
! (Note: We don't know correct values for g%df. That's exactly what we are going to compute.)
write(*,*)
write(*,*) "Seeding x..."
call ddt_seed(x) ! This will set x(i)%df(i) = 1.0 for i=1,2,3.
write(*,*)
write(*,'(a,es25.15,a,3es25.15)') " x(1)%f=",x(1)%f, " x(1)%df(1), x(1)%df(2), x(1)%df(3) =",x(1)%df
write(*,'(a,es25.15,a,3es25.15)') " x(2)%f=",x(2)%f, " x(2)%df(1), x(2)%df(2), x(2)%df(3) =",x(2)%df
write(*,'(a,es25.15,a,3es25.15)') " x(3)%f=",x(3)%f, " x(3)%df(1), x(3)%df(2), x(3)%df(3) =",x(3)%df
write(*,*)
write(*,*) "Seeding done."
! Compute the derivative of function g: dg/dx by calling the function g.
write(*,*)
write(*,*) "Calling a function subroutine that computes g... 'g = function_g(x)'"
g = function_g(x)
! Display the derivatives computed by AD.
write(*,*) "Returned from the subroutine."
write(*,*)
write(*,*) "Now, you don't see it, but dg/dx has been computed and stored in g%df."
write(*,*) "See below."
write(*,*)
write(*,*) "Computed dg/dx:"
write(*,'(a,3es25.15)') " g(1)%df(1), g(1)%df(2), g(1)%df(3) = ", g(1)%df
write(*,'(a,3es25.15)') " g(2)%df(1), g(2)%df(2), g(2)%df(3) = ", g(2)%df
write(*,'(a,3es25.15)') " g(3)%df(1), g(3)%df(2), g(3)%df(3) = ", g(3)%df
! Below are the exact derivatives of the funciton defined in function_g(x), worked out by hand.
! Note: these are correct for x=(2,5,7), but may not be correct for other values because
! the function g contains min, max, sign functions. See function_g(x).
exact_derivatives(1,1) = 2.0_p2*x_real(1)
exact_derivatives(1,2) = 5.0_p2+0.5_p2/sqrt(x_real(2))
exact_derivatives(1,3) = 2.0_p2**x_real(3)*log(2.0_p2)
exact_derivatives(2,1) = 0.0_p2
exact_derivatives(2,2) = 1.0_p2+2.0_p2*x_real(2)
exact_derivatives(2,3) = 6.0_p2
exact_derivatives(3,1) = 7.0_p2 + 1.0_p2/x_real(3)
exact_derivatives(3,2) = 0.0_p2
exact_derivatives(3,3) = 2.0_p2*x_real(3)-x_real(1)/x_real(3)**2
! Display the exact derivatives.
write(*,*)
write(*,*) "Exact dg/dx:"
write(*,'(a,3es25.15)') " dg1/dx1, dg1/dx2, dg1/dx3 = ", &
exact_derivatives(1,1),exact_derivatives(1,2),exact_derivatives(1,3)
write(*,'(a,3es25.15)') " dg2/dx1, dg2/dx2, dg2/dx3 = ", &
exact_derivatives(2,1),exact_derivatives(2,2),exact_derivatives(2,3)
write(*,'(a,3es25.15)') " dg3/dx1, dg3/dx2, dg3/dx3 = ", &
exact_derivatives(3,1),exact_derivatives(3,2),exact_derivatives(3,3)
write(*,*)
write(*,*) "Difference between the exact and computed derivatives:"
write(*,*)
write(*,*) (g(1)%df(j)-exact_derivatives(1,j), j=1,3)
write(*,*) (g(2)%df(j)-exact_derivatives(2,j), j=1,3)
write(*,*) (g(3)%df(j)-exact_derivatives(3,j), j=1,3)
write(*,*)
write(*,*) "Zero, zero, zero! Congratulations! AD is working!"
write(*,*)
stop
contains
!********************************************************************************
!* This is a function that computes the value of g(x) for a given x.
!*
!* Note: g and x are vectors of dimension 3.
!*
!* Note: The derivative comptued, dg/dx, is a matrix.
!*
!* Note: This is just a function, but written in terms of derivative data type,
!* so that the derivatives are computed simultaneously.
!* At the time of return, the derivatives are computed (behind the scene)
!* and stored in function_g%df (which corresponds exactly to dg/dx).
!* So, it is automatic! You don't have to do anything but compute the
!* function. It's nice, isn't it?
!*
!* Note: This function needs to be called with x having x(i)%df(i)=1.0 and
!* x(i)%df(j)=0.0 for i not equal to j, i.e., derivatives of x need to be
!* set correctly.
!*
!* This function is used in the example program above.
!********************************************************************************
function function_g(x)
use derivative_data_df3
implicit none
integer , parameter :: p2 = selected_real_kind(15)
type(derivative_data_type_df3), dimension(3), intent(in) :: x
type(derivative_data_type_df3), dimension(3) :: function_g
! Compute the function g.
function_g(1) = x(1)**2 + x(2)*5.0_p2 + 2.0_p2**x(3) + ddt_sqrt(x(2))
function_g(2) = x(2)**2 + x(3)*5.0_p2 + ddt_max(x(1),x(2)) + ddt_sign(1.0_p2,x(2)) + ddt_abs(x(3))
function_g(3) = x(3)**2 + x(1)*5.0_p2 + ddt_min(x(1),x(3)) + x(1)/x(3) + ddt_min(x(1),x(3))
end function function_g
end program test_automatic_differentiation