! D02AGF Example Program Text ! Mark 23 Release. NAG Copyright 2011. 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) ! .. Implicit None Statement .. IMPLICIT NONE ! .. 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) ! .. Implicit None Statement .. IMPLICIT NONE ! .. 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) ! .. Implicit None Statement .. IMPLICIT NONE ! .. 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 Functions .. 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) ! .. Implicit None Statement .. IMPLICIT NONE ! .. 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) ! .. Implicit None Statement .. IMPLICIT NONE ! .. 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 Functions .. 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) ! .. Implicit None Statement .. IMPLICIT NONE ! .. 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) ! .. Implicit None Statement .. IMPLICIT NONE ! .. 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) ! .. Implicit None Statement .. IMPLICIT NONE ! .. 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 USE d02agfe_mod, ONLY : nin USE d02agfe_ex1, ONLY : aux1, bcaux1, n, n1, nag_wp, prsol1, raaux1 ! .. Implicit None Statement .. IMPLICIT NONE ! .. 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 Functions .. 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 USE d02agfe_mod, ONLY : nin USE d02agfe_ex2, ONLY : aux2, bcaux2, n, n1, nag_wp, prsol2, raaux2 ! .. Implicit None Statement .. IMPLICIT NONE ! .. 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