!   D02AGF Example Program Text
!   Mark 26.1 Release. NAG Copyright 2016.

    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 d02agfe_mod, Only: aux1, bcaux1, n1_ex1, nin, n_ex1, prsol1,       &
                               raaux1
        Use nag_library, Only: d02agf, nag_wp
!       .. 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 d02agfe_mod, Only: aux2, bcaux2, n1_ex2, nin, n_ex2, prsol2,       &
                               raaux2
        Use nag_library, Only: d02agf, nag_wp
!       .. 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