! D02BJ_A1T1W_F Example Program Text
! Mark 28.3 Release. NAG Copyright 2022.
Module d02bj_a1t1w_fe_mod
! Data for D02BJ_A1T1W_F example program
! .. Use Statements ..
Use iso_c_binding, Only: c_ptr
Use nagad_library, Only: cos, nagad_a1t1w_w_rtype, tan, Assignment (=), &
Operator (+), Operator (*), Operator (/), &
Operator (**)
! .. Implicit None Statement ..
Implicit None
! .. Accessibility Statements ..
Private
Public :: fcn, g
! .. Parameters ..
Integer, Parameter, Public :: n = 3, nin = 5, nout = 6
! n: number of differential equations
Contains
Subroutine fcn(ad_handle,x,y,f,iuser,ruser)
! .. Scalar Arguments ..
Type (c_ptr), Intent (Inout) :: ad_handle
Type (nagad_a1t1w_w_rtype), Intent (In) :: x
! .. Array Arguments ..
Type (nagad_a1t1w_w_rtype), Intent (Inout) :: f(*), ruser(*)
Type (nagad_a1t1w_w_rtype), Intent (In) :: y(*)
Integer, Intent (Inout) :: iuser(*)
! .. Local Scalars ..
Type (nagad_a1t1w_w_rtype) :: alpha, beta
! .. Executable Statements ..
alpha = ruser(1)
beta = ruser(2)
f(1) = tan(y(3))
f(2) = alpha*tan(y(3))/y(2) + beta*y(2)/cos(y(3))
f(3) = alpha/y(2)**2
Return
End Subroutine fcn
Subroutine g(ad_handle,x,y,retval,iuser,ruser)
! .. Scalar Arguments ..
Type (c_ptr), Intent (Inout) :: ad_handle
Type (nagad_a1t1w_w_rtype), Intent (Out) :: retval
Type (nagad_a1t1w_w_rtype), Intent (In) :: x
! .. Array Arguments ..
Type (nagad_a1t1w_w_rtype), Intent (Inout) :: ruser(*)
Type (nagad_a1t1w_w_rtype), Intent (In) :: y(*)
Integer, Intent (Inout) :: iuser(*)
! .. Executable Statements ..
retval = y(1)
Return
End Subroutine g
End Module d02bj_a1t1w_fe_mod
Program d02bj_a1t1w_fe
! D02BJ_A1T1W_F Example Main Program
! .. Use Statements ..
Use d02bj_a1t1w_fe_mod, Only: fcn, g, n, nin, nout
Use iso_c_binding, Only: c_ptr
Use nagad_library, Only: d02bj_a1t1w_f, d02bj_a1t1w_x, &
nagad_a1t1w_get_derivative, &
nagad_a1t1w_inc_derivative, &
nagad_a1t1w_ir_interpret_adjoint, &
nagad_a1t1w_ir_register_variable, &
nagad_a1t1w_ir_remove, nagad_a1t1w_w_rtype, &
nagad_t1w_w_rtype, x10aa_a1t1w_f, &
x10ab_a1t1w_f, x10za_a1t1w_f, Assignment (=)
Use nag_library, Only: nag_wp
! .. Implicit None Statement ..
Implicit None
! .. Parameters ..
Real (Kind=nag_wp), Parameter :: alpha = -0.032E0_nag_wp
Real (Kind=nag_wp), Parameter :: beta = -0.02E0_nag_wp
! .. Local Scalars ..
Type (c_ptr) :: ad_handle
Type (nagad_a1t1w_w_rtype) :: tol, x, xend_ad
Type (nagad_t1w_w_rtype) :: t_t
Real (Kind=nag_wp) :: da, xend, xinit
Integer :: i, ifail, iw, kinit
! .. Local Arrays ..
Type (nagad_a1t1w_w_rtype) :: ruser(4)
Type (nagad_a1t1w_w_rtype), Allocatable :: w(:), y(:), y_in(:)
Real (Kind=nag_wp), Allocatable :: yinit(:)
Integer :: iuser(1)
! .. Intrinsic Procedures ..
Intrinsic :: real
! .. Executable Statements ..
Write (nout,*) 'D02BJ_A1T1W_F Example Program Results'
iw = 20*n
Allocate (w(iw),y(n),yinit(n),y_in(n))
! Skip heading in data file
Read (nin,*)
! xinit: initial x value, xend: final x value.
! yinit: initial solution values
Read (nin,*) xinit, xend
Read (nin,*) yinit(1:n)
Read (nin,*) kinit
Write (nout,99999) 'no intermediate output, root-finding'
99999 Format (1X,'Case : ',A)
tol = 1.0E-5_nag_wp
Write (nout,*)
Write (nout,99998) ' Calculation with TOL =', tol%value%value
99998 Format (1X,A,E8.1)
! Create AD tape
Call x10za_a1t1w_f
! Create AD configuration data object
ifail = 0
Call x10aa_a1t1w_f(ad_handle,ifail)
x = xinit
xend_ad = xend
y_in(1:n) = yinit(1:n)
ruser(1) = alpha
ruser(2) = beta
ruser(3) = (xend-xinit)/real(kinit+1,kind=nag_wp)
ruser(4) = xend
ruser(1:2)%value%tangent = 1.0_nag_wp
y_in(1:n)%value%tangent = 1.0_nag_wp
! Register variables to differentiate w.r.t.
Call nagad_a1t1w_ir_register_variable(ruser(1:2))
Call nagad_a1t1w_ir_register_variable(y_in)
y(1:n) = y_in(1:n)
ifail = 0
Call d02bj_a1t1w_f(ad_handle,x,xend_ad,n,y,fcn,tol,'Default', &
d02bj_a1t1w_x,g,w,iuser,ruser,ifail)
Write (nout,99997) ' Root of Y(1) = 0.0 at', x%value%value
Write (nout,99996) ' Solution is', (y(i)%value%value,i=1,n)
Write (nout,*)
99997 Format (1X,A,F7.3)
99996 Format (1X,A,3F13.4)
Write (nout,*)
Write (nout,*) &
' Derivatives calculated: Second order, adjoints of tangents'
Write (nout,*) ' Computational mode : algorithmic'
Write (nout,*)
Write (nout,*) ' Derivative: hit point w.r.t. all parameters summed'
! Setup evaluation of derivatives via adjoints
t_t%value = 1.0_nag_wp
t_t%tangent = 0.0_nag_wp
Call nagad_a1t1w_inc_derivative(x,t_t)
ifail = 0
Call nagad_a1t1w_ir_interpret_adjoint(ifail)
t_t = nagad_a1t1w_get_derivative(ruser(1))
da = t_t%tangent
t_t = nagad_a1t1w_get_derivative(ruser(2))
da = da + t_t%tangent
Do i = 1, n
t_t = nagad_a1t1w_get_derivative(y_in(i))
da = da + t_t%tangent
End Do
Write (nout,99996) ' Sum of Hessian terms for x = ', da
! Remove computational data object and tape
Call x10ab_a1t1w_f(ad_handle,ifail)
Call nagad_a1t1w_ir_remove
End Program d02bj_a1t1w_fe