! D02BJ_T1W_F Example Program Text
! Mark 28.7 Release. NAG Copyright 2022.
Module d02bj_t1w_fe_mod
! Data for D02BJ_T1W_F example program
! .. Use Statements ..
Use iso_c_binding, Only: c_ptr
Use nagad_library, Only: cos, nagad_t1w_w_rtype, tan, 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_t1w_w_rtype), Intent (In) :: x
! .. Array Arguments ..
Type (nagad_t1w_w_rtype), Intent (Inout) :: f(*), ruser(*)
Type (nagad_t1w_w_rtype), Intent (In) :: y(*)
Integer, Intent (Inout) :: iuser(*)
! .. Local Scalars ..
Type (nagad_t1w_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_t1w_w_rtype), Intent (Out) :: retval
Type (nagad_t1w_w_rtype), Intent (In) :: x
! .. Array Arguments ..
Type (nagad_t1w_w_rtype), Intent (Inout) :: ruser(*)
Type (nagad_t1w_w_rtype), Intent (In) :: y(*)
Integer, Intent (Inout) :: iuser(*)
! .. Executable Statements ..
retval = y(1)
Return
End Subroutine g
End Module d02bj_t1w_fe_mod
Program d02bj_t1w_fe
! D02BJ_T1W_F Example Main Program
! .. Use Statements ..
Use d02bj_t1w_fe_mod, Only: fcn, g, n, nin, nout
Use iso_c_binding, Only: c_ptr
Use nagad_library, Only: d02bj_t1w_f, d02bj_t1w_x, nagad_t1w_w_rtype, &
x10aa_t1w_f, x10ab_t1w_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_t1w_w_rtype) :: tol, x, xend_ad
Real (Kind=nag_wp) :: da, db, xend, xinit
Integer :: i, ifail, iw, kinit
! .. Local Arrays ..
Type (nagad_t1w_w_rtype) :: ruser(4)
Type (nagad_t1w_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_T1W_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,99996) 'no intermediate output, root-finding'
tol = 1.0E-5_nag_wp
Write (nout,*)
Write (nout,99999) ' Calculation with TOL =', tol%value
! Create AD configuration data object
ifail = 0
Call x10aa_t1w_f(ad_handle,ifail)
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
Do i = 1, 5
x = xinit
xend_ad = xend
y(1:n) = y_in(1:n)
If (i<3) Then
ruser(i)%tangent = 1.0_nag_wp
Else
y(i-2)%tangent = 1.0_nag_wp
End If
ifail = 0
Call d02bj_t1w_f(ad_handle,x,xend_ad,n,y,fcn,tol,'Default', &
d02bj_t1w_x,g,w,iuser,ruser,ifail)
If (i<3) Then
ruser(i)%tangent = 0.0_nag_wp
If (i==1) Then
da = x%tangent
Else
db = x%tangent
End If
Else
y(i-2)%tangent = 0.0_nag_wp
yinit(i-2) = x%tangent
End If
End Do
Write (nout,99998) ' Root of Y(1) = 0.0 at', x%value
Write (nout,99997) ' Solution is', (y(i)%value,i=1,n)
Write (nout,*)
99999 Format (1X,A,E8.1)
99998 Format (1X,A,F7.3)
99997 Format (1X,A,3F13.4)
99996 Format (1X,'Case : ',A)
Write (nout,*)
Write (nout,*) ' Derivatives calculated: First order tangents'
Write (nout,*) ' Computational mode : algorithmic'
Write (nout,*)
Write (nout,*) ' Derivatives: (hit point w.r.t. parameters)'
Write (nout,99995) ' dx/dg =', -da
Write (nout,99995) ' dx/ddrag =', -db
Write (nout,99995) ' dx/dheight =', yinit(1)
Write (nout,99995) ' dx/dvel =', yinit(2)
Write (nout,99995) ' dx/dangle =', yinit(3)
99995 Format (1X,A,1X,E12.5)
! Remove computational data object
ifail = 0
Call x10ab_t1w_f(ad_handle,ifail)
End Program d02bj_t1w_fe