! D02AGF Example Program Text ! Mark 24 Release. NAG Copyright 2012. Module d02agfe_mod ! D02AGF Example Program Module: ! Global Parameters ! iprint: set iprint = 1 for output at each Newton iteration. ! nin: the input channel number ! nout: the output channel number ! .. Use Statements .. Use nag_library, Only: nag_wp ! .. Implicit None Statement .. Implicit None ! .. Parameters .. Integer, Parameter :: iprint = 0, nin = 5, nout = 6 End Module d02agfe_mod Module d02agfe_ex1 ! D02AGF Example Program Module: ! Parameters and User-defined Routines For Example 1 ! n : number of differential equations, ! n1: number of parameters. ! .. Use Statements .. Use nag_library, Only: nag_wp Use d02agfe_mod, Only: iprint, nout ! .. Implicit None Statement .. Implicit None ! .. Parameters .. Integer, Parameter :: n = 2, n1 = 2 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 End Module d02agfe_ex1 Module d02agfe_ex2 ! D02AGF Example Program Module: ! Parameters and User-defined Routines For Example 2 ! n : number of differential equations, ! n1: number of parameters. ! .. Use Statements .. Use nag_library, Only: nag_wp Use d02agfe_mod, Only: iprint, nout ! .. Implicit None Statement .. Implicit None ! .. Parameters .. Integer, Parameter :: n = 3, n1 = 3 Contains 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_ex2 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: nin Use d02agfe_ex1, Only: aux1, bcaux1, n, n1, 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),e(n),mat(n,n),param(n),parerr(n),wspac1(n), & wspac2(n),wspace(n,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) Read (nin,*) parerr(1:n1) Read (nin,*) soler e(1:n) = 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,n1,m1,aux1,bcaux1,raaux1,prsol1,mat, & dummy,wspace,wspac1,wspac2,ifail) Write (nout,*) 'Final parameters' Write (nout,99999)(param(i),i=1,n1) 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) 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: nin Use d02agfe_ex2, Only: aux2, bcaux2, n, n1, 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),e(n),mat(n,n),param(n),parerr(n),wspac1(n), & wspac2(n),wspace(n,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) Read (nin,*) parerr(1:n1) Read (nin,*) soler e(1:n) = 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,n1,m1,aux2,bcaux2,raaux2,prsol2,mat, & dummy,wspace,wspac1,wspac2,ifail) Write (nout,*) 'Final parameters' Write (nout,99999)(param(i),i=1,n) 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) xx = xx + h End Do Return 99999 Format (1X,3E16.6) 99998 Format (1X,F7.0,3E13.4) End Subroutine ex2 End Program d02agfe