! D02AGF Example Program Text
! Mark 25 Release. NAG Copyright 2014.
Module d02agfe_mod
! D02AGF Example Program Module:
! Parameters and User-defined Routines
! iprint: set iprint = 1 for output at each Newton iteration.
! nin: the input channel number
! nout: the output channel number
! For Example 1:
! n_ex1 : number of differential equations,
! n1_ex1: number of parameters.
! For Example 2:
! n_ex2 : number of differential equations,
! n1_ex2: number of parameters.
! .. Use Statements ..
Use nag_library, Only: nag_wp
! .. Implicit None Statement ..
Implicit None
! .. Accessibility Statements ..
Private
Public :: aux1, aux2, bcaux1, bcaux2, &
prsol1, prsol2, raaux1, raaux2
! .. Parameters ..
Integer, Parameter :: iprint = 0
Integer, Parameter, Public :: n1_ex1 = 2, n1_ex2 = 3, nin = 5, &
nout = 6, n_ex1 = 2, n_ex2 = 3
Contains
Subroutine aux1(f,y,x,param)
! .. Scalar Arguments ..
Real (Kind=nag_wp), Intent (In) :: x
! .. Array Arguments ..
Real (Kind=nag_wp), Intent (Out) :: f(*)
Real (Kind=nag_wp), Intent (In) :: param(*), y(*)
! .. Executable Statements ..
f(1) = y(2)
f(2) = (y(1)**3-y(2))/(2.0_nag_wp*x)
Return
End Subroutine aux1
Subroutine raaux1(x0,x1,r,param)
! .. Scalar Arguments ..
Real (Kind=nag_wp), Intent (Out) :: r, x0, x1
! .. Array Arguments ..
Real (Kind=nag_wp), Intent (In) :: param(*)
! .. Executable Statements ..
x0 = 0.1_nag_wp
x1 = 16.0_nag_wp
r = 16.0_nag_wp
Return
End Subroutine raaux1
Subroutine bcaux1(g0,g1,param)
! .. Array Arguments ..
Real (Kind=nag_wp), Intent (Out) :: g0(*), g1(*)
Real (Kind=nag_wp), Intent (In) :: param(*)
! .. Local Scalars ..
Real (Kind=nag_wp) :: z
! .. Intrinsic Procedures ..
Intrinsic :: sqrt
! .. Executable Statements ..
z = 0.1_nag_wp
g0(1) = 0.1_nag_wp + param(1)*sqrt(z)*0.1_nag_wp + 0.01_nag_wp*z
g0(2) = param(1)*0.05_nag_wp/sqrt(z) + 0.01_nag_wp
g1(1) = 1.0_nag_wp/6.0_nag_wp
g1(2) = param(2)
Return
End Subroutine bcaux1
Subroutine prsol1(param,res,n1,err)
! .. Scalar Arguments ..
Real (Kind=nag_wp), Intent (In) :: res
Integer, Intent (In) :: n1
! .. Array Arguments ..
Real (Kind=nag_wp), Intent (In) :: err(n1), param(n1)
! .. Local Scalars ..
Integer :: i
! .. Executable Statements ..
If (iprint/=0) Then
Write (nout,99999) 'Current parameters ', (param(i),i=1,n1)
Write (nout,99998) 'Residuals ', (err(i),i=1,n1)
Write (nout,99998) 'Sum of residuals squared ', res
Write (nout,*)
End If
Return
99999 Format (1X,A,6(E14.6,2X))
99998 Format (1X,A,6(E12.4,1X))
End Subroutine prsol1
Subroutine aux2(f,y,x,param)
! .. Parameters ..
Real (Kind=nag_wp), Parameter :: eps = 2.0E-5_nag_wp
! .. Scalar Arguments ..
Real (Kind=nag_wp), Intent (In) :: x
! .. Array Arguments ..
Real (Kind=nag_wp), Intent (Out) :: f(*)
Real (Kind=nag_wp), Intent (In) :: param(*), y(*)
! .. Local Scalars ..
Real (Kind=nag_wp) :: c, s
! .. Intrinsic Procedures ..
Intrinsic :: cos, sin
! .. Executable Statements ..
c = cos(y(3))
s = sin(y(3))
f(1) = s/c
f(2) = -(param(1)*s+eps*y(2)*y(2))/(y(2)*c)
f(3) = -param(1)/(y(2)*y(2))
Return
End Subroutine aux2
Subroutine raaux2(x0,x1,r,param)
! .. Scalar Arguments ..
Real (Kind=nag_wp), Intent (Out) :: r, x0, x1
! .. Array Arguments ..
Real (Kind=nag_wp), Intent (In) :: param(*)
! .. Executable Statements ..
x0 = 0.0_nag_wp
x1 = param(2)
r = param(2)
Return
End Subroutine raaux2
Subroutine bcaux2(g0,g1,param)
! .. Array Arguments ..
Real (Kind=nag_wp), Intent (Out) :: g0(*), g1(*)
Real (Kind=nag_wp), Intent (In) :: param(*)
! .. Executable Statements ..
g0(1) = 0.0E0_nag_wp
g0(2) = 500.0E0_nag_wp
g0(3) = 0.5E0_nag_wp
g1(1) = 0.0E0_nag_wp
g1(2) = 450.0E0_nag_wp
g1(3) = param(3)
Return
End Subroutine bcaux2
Subroutine prsol2(param,res,n1,err)
! .. Scalar Arguments ..
Real (Kind=nag_wp), Intent (In) :: res
Integer, Intent (In) :: n1
! .. Array Arguments ..
Real (Kind=nag_wp), Intent (In) :: err(n1), param(n1)
! .. Local Scalars ..
Integer :: i
! .. Executable Statements ..
If (iprint/=0) Then
Write (nout,99999) 'Current parameters ', (param(i),i=1,n1)
Write (nout,99998) 'Residuals ', (err(i),i=1,n1)
Write (nout,99998) 'Sum of residuals squared ', res
Write (nout,*)
End If
Return
99999 Format (1X,A,6(E14.6,2X))
99998 Format (1X,A,6(E12.4,1X))
End Subroutine prsol2
End Module d02agfe_mod
Program d02agfe
! D02AGF Example Main Program
! .. Use Statements ..
Use d02agfe_mod, Only: nout
! .. Implicit None Statement ..
Implicit None
! .. Executable Statements ..
Write (nout,*) 'D02AGF Example Program Results'
Call ex1
Call ex2
Contains
Subroutine ex1
! .. Use Statements ..
Use nag_library, Only: d02agf, nag_wp
Use d02agfe_mod, Only: aux1, bcaux1, n1_ex1, nin, n_ex1, prsol1, raaux1
! .. Local Scalars ..
Real (Kind=nag_wp) :: h, r, soler, x, x1, xx
Integer :: i, ifail, j, m1
! .. Local Arrays ..
Real (Kind=nag_wp), Allocatable :: c(:,:), e(:), mat(:,:), &
param(:), parerr(:), &
wspac1(:), wspac2(:), &
wspace(:,:)
Real (Kind=nag_wp) :: dummy(1,1)
! .. Intrinsic Procedures ..
Intrinsic :: real
! .. Executable Statements ..
! Skip heading in data file
Read (nin,*)
! m1: final solution calculated at m1 points in range including
! end points.
Read (nin,*) m1
Allocate (c(m1,n_ex1),e(n_ex1),mat(n_ex1,n_ex1),param(n_ex1), &
parerr(n_ex1),wspac1(n_ex1),wspac2(n_ex1),wspace(n_ex1,9))
! h: Step size estimate, param: initial parameter estimates
! parerr: Newton iteration tolerances, soler: bound on the local error.
Read (nin,*) h
Read (nin,*) param(1:n1_ex1)
Read (nin,*) parerr(1:n1_ex1)
Read (nin,*) soler
e(1:n_ex1) = soler
Write (nout,*)
Write (nout,*)
Write (nout,*) 'Case 1'
Write (nout,*)
! ifail: behaviour on error exit
! =0 for hard exit, =1 for quiet-soft, =-1 for noisy-soft
ifail = 0
Call d02agf(h,e,parerr,param,c,n_ex1,n1_ex1,m1,aux1,bcaux1,raaux1, &
prsol1,mat,dummy,wspace,wspac1,wspac2,ifail)
Write (nout,*) 'Final parameters'
Write (nout,99999)(param(i),i=1,n1_ex1)
Write (nout,*)
Write (nout,*) 'Final solution'
Write (nout,*) 'X-value Components of solution'
Call raaux1(x,x1,r,param)
h = (x1-x)/real(m1-1,kind=nag_wp)
xx = x
Do i = 1, m1
Write (nout,99998) xx, (c(i,j),j=1,n_ex1)
xx = xx + h
End Do
Return
99999 Format (1X,3E16.6)
99998 Format (1X,F7.2,3E13.4)
End Subroutine ex1
Subroutine ex2
! .. Use Statements ..
Use nag_library, Only: d02agf, nag_wp
Use d02agfe_mod, Only: aux2, bcaux2, n1_ex2, nin, n_ex2, prsol2, raaux2
! .. Local Scalars ..
Real (Kind=nag_wp) :: h, r, soler, x, x1, xx
Integer :: i, ifail, j, m1
! .. Local Arrays ..
Real (Kind=nag_wp), Allocatable :: c(:,:), e(:), mat(:,:), &
param(:), parerr(:), &
wspac1(:), wspac2(:), &
wspace(:,:)
Real (Kind=nag_wp) :: dummy(1,1)
! .. Executable Statements ..
Read (nin,*)
! Read in problem parameters
! m1: final solution calculated at m1 points in range including
! end points.
Read (nin,*) m1
Allocate (c(m1,n_ex2),e(n_ex2),mat(n_ex2,n_ex2),param(n_ex2), &
parerr(n_ex2),wspac1(n_ex2),wspac2(n_ex2),wspace(n_ex2,9))
Write (nout,*)
Write (nout,*)
Write (nout,*) 'Case 2'
Write (nout,*)
! h: Step size estimate, param: initial parameter estimates
! parerr: Newton iteration tolerances, soler: bound on the local error.
Read (nin,*) h
Read (nin,*) param(1:n1_ex2)
Read (nin,*) parerr(1:n1_ex2)
Read (nin,*) soler
e(1:n_ex2) = soler
! ifail: behaviour on error exit
! =0 for hard exit, =1 for quiet-soft, =-1 for noisy-soft
ifail = 0
Call d02agf(h,e,parerr,param,c,n_ex2,n1_ex2,m1,aux2,bcaux2,raaux2, &
prsol2,mat,dummy,wspace,wspac1,wspac2,ifail)
Write (nout,*) 'Final parameters'
Write (nout,99999)(param(i),i=1,n_ex2)
Write (nout,*)
Write (nout,*) 'Final solution'
Write (nout,*) 'X-value Components of solution'
Call raaux2(x,x1,r,param)
h = (x1-x)/5.0E0_nag_wp
xx = x
Do i = 1, 6
Write (nout,99998) xx, (c(i,j),j=1,n_ex2)
xx = xx + h
End Do
Return
99999 Format (1X,3E16.6)
99998 Format (1X,F7.0,3E13.4)
End Subroutine ex2
End Program d02agfe