!   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
