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

    Module d02hbfe_mod

!     Data for D02HBF example programs

!     .. Use Statements ..
      Use nag_library, Only: nag_wp
!     .. Implicit None Statement ..
      Implicit None
!     .. Accessibility Statements ..
      Private
      Public                           :: bc1, bc2, fcn1, fcn2, range1, range2
!     .. Parameters ..
      Integer, Parameter, Public       :: iset = 1, nin = 5, nout = 6
    Contains
      Subroutine fcn1(x,y,f,p)

!       .. Scalar Arguments ..
        Real (Kind=nag_wp), Intent (In) :: x
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (Out) :: f(*)
        Real (Kind=nag_wp), Intent (In) :: p(*), y(*)
!       .. Executable Statements ..
        f(1) = y(2)
        f(2) = (y(1)**3-y(2))/(2.0E0_nag_wp*x)
        Return
      End Subroutine fcn1
      Subroutine range1(a,b,p)

!       .. Scalar Arguments ..
        Real (Kind=nag_wp), Intent (Out) :: a, b
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (In) :: p(*)
!       .. Executable Statements ..
        a = 0.1E0_nag_wp
        b = 16.0E0_nag_wp
        Return
      End Subroutine range1
      Subroutine bc1(g1,g2,p)

!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (Out) :: g1(*), g2(*)
        Real (Kind=nag_wp), Intent (In) :: p(*)
!       .. Local Scalars ..
        Real (Kind=nag_wp)             :: z
!       .. Intrinsic Procedures ..
        Intrinsic                      :: sqrt
!       .. Executable Statements ..
        z = 0.1E0_nag_wp
        g1(1) = 0.1E0_nag_wp + p(1)*sqrt(z)*0.1E0_nag_wp + 0.01E0_nag_wp*z
        g1(2) = p(1)*0.05E0_nag_wp/sqrt(z) + 0.01E0_nag_wp
        g2(1) = 1.0E0_nag_wp/6.0E0_nag_wp
        g2(2) = p(2)
        Return
      End Subroutine bc1
      Subroutine fcn2(x,y,f,p)

!       .. Scalar Arguments ..
        Real (Kind=nag_wp), Intent (In) :: x
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (Out) :: f(*)
        Real (Kind=nag_wp), Intent (In) :: p(*), y(*)
!       .. Intrinsic Procedures ..
        Intrinsic                      :: cos, tan
!       .. Executable Statements ..
        f(1) = tan(y(3))
        f(2) = -p(1)*tan(y(3))/y(2) - 0.00002E0_nag_wp*y(2)/cos(y(3))
        f(3) = -p(1)/y(2)**2
        Return
      End Subroutine fcn2
      Subroutine range2(a,b,p)

!       .. Scalar Arguments ..
        Real (Kind=nag_wp), Intent (Out) :: a, b
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (In) :: p(*)
!       .. Executable Statements ..
        a = 0.0E0_nag_wp
        b = p(2)
        Return
      End Subroutine range2
      Subroutine bc2(g1,g2,p)

!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (Out) :: g1(*), g2(*)
        Real (Kind=nag_wp), Intent (In) :: p(*)
!       .. Executable Statements ..
        g1(1) = 0.0E0_nag_wp
        g1(2) = 500.0E0_nag_wp
        g1(3) = 0.5E0_nag_wp
        g2(1) = 0.0E0_nag_wp
        g2(2) = 450.0E0_nag_wp
        g2(3) = p(3)
        Return
      End Subroutine bc2
    End Module d02hbfe_mod
    Program d02hbfe

!     D02HBF Example Main Program

!     .. Use Statements ..
      Use d02hbfe_mod, Only: nout
!     .. Implicit None Statement ..
      Implicit None
!     .. Executable Statements ..
      Write (nout,*) 'D02HBF Example Program Results'

      Call ex1

      Call ex2

    Contains
      Subroutine ex1

!       .. Use Statements ..
        Use d02hbfe_mod, Only: bc1, fcn1, iset, nin, range1
        Use nag_library, Only: d02hbf, nag_wp, x04abf
!       .. Local Scalars ..
        Real (Kind=nag_wp)             :: h, x, x1, xh
        Integer                        :: i, ifail, m1, n, n1, outchn, sdw
!       .. Local Arrays ..
        Real (Kind=nag_wp), Allocatable :: e(:), p(:), pe(:), soln(:,:),       &
                                          w(:,:)
!       .. Intrinsic Procedures ..
        Intrinsic                      :: real
!       .. Executable Statements ..
!       Skip heading in data file
        Read (nin,*)
!       m1: controls exit values,  n: number of differential equations,
!       n1: number of parameters.
        Read (nin,*) m1, n, n1
        sdw = 3*n + 14 + 11
        Allocate (e(n),p(n1),pe(n1),soln(n,m1),w(n,sdw))
        Write (nout,*)
        outchn = nout
        Write (nout,*)
        Call x04abf(iset,outchn)
!       p: estimates for the parameters p,  e: bound on the local error.
        Read (nin,*) p(1:n1)
        Read (nin,*) pe(1:n1)
        Read (nin,*) e(1:n)

        Write (nout,*) 'Case 1'
        Write (nout,*)
!       ifail: behaviour on error exit
!              =1 for quiet-soft exit
!       * Set ifail to 111 to obtain monitoring information *
        ifail = 1
        Call d02hbf(p,n1,pe,e,n,soln,m1,fcn1,bc1,range1,w,sdw,ifail)

        If (ifail==0) Then
          Write (nout,*) 'Final parameters'
          Write (nout,99999)(p(i),i=1,n1)
          Write (nout,*)
          Write (nout,*) 'Final solution'
          Write (nout,*) 'X-value     Components of solution'
          Call range1(x,x1,p)
          h = (x1-x)/real(m1-1,kind=nag_wp)
          xh = x
          Do i = 1, m1
            Write (nout,99998) xh, soln(1:n,i)
            xh = xh + h
          End Do
        Else
          Write (nout,99996) ifail
          If (ifail>1 .And. ifail<=5) Then
            Write (nout,99997) w(1,2), (w(i,1),i=1,n)
          End If
        End If

        Return

99999   Format (1X,1P,3E15.3)
99998   Format (1X,F7.2,2F13.4)
99997   Format (/,1X,'W(1,2) = ',F9.4,' W(.,1) = ',10E10.3)
99996   Format (1X,/,1X,' ** D02HBF returned with IFAIL = ',I5)
      End Subroutine ex1
      Subroutine ex2

!       .. Use Statements ..
        Use d02hbfe_mod, Only: bc2, fcn2, iset, nin, range2
        Use nag_library, Only: d02hbf, nag_wp, x04abf
!       .. Local Scalars ..
        Real (Kind=nag_wp)             :: h, x, x1, xh
        Integer                        :: i, ifail, m1, n, n1, outchn, sdw
!       .. Local Arrays ..
        Real (Kind=nag_wp), Allocatable :: e(:), p(:), pe(:), soln(:,:),       &
                                          w(:,:)
!       .. Intrinsic Procedures ..
        Intrinsic                      :: real
!       .. Executable Statements ..
        Read (nin,*)
!       m1: controls exit values,  n: number of differential equations,
!       n1: number of parameters.
        Read (nin,*) m1, n, n1
        sdw = 3*n + 14 + 11
        Allocate (e(n),p(n1),pe(n1),soln(n,m1),w(n,sdw))
        outchn = nout
        Call x04abf(iset,outchn)
!       p: estimates for the parameters p,  e: bound on the local error.
        Read (nin,*) p(1:n1)
        Read (nin,*) pe(1:n1)
        Read (nin,*) e(1:n)

        Write (nout,*)
        Write (nout,*)
        Write (nout,*) 'Case 2'
        Write (nout,*)
!       ifail: behaviour on error exit
!              =1 for quiet-soft exit
!       * Set ifail to 111 to obtain monitoring information *
        ifail = 1
        Call d02hbf(p,n1,pe,e,n,soln,m1,fcn2,bc2,range2,w,sdw,ifail)

        If (ifail==0) Then
          Write (nout,*) 'Final parameters'
          Write (nout,99999)(p(i),i=1,n1)
          Write (nout,*)
          Write (nout,*) 'Final solution'
          Write (nout,*) 'X-value     Components of solution'
          Call range2(x,x1,p)
          h = (x1-x)/real(m1-1,kind=nag_wp)
          xh = x
          Do i = 1, m1
            Write (nout,99998) xh, soln(1:n,i)
            xh = xh + h
          End Do
        Else
          Write (nout,99996) ifail
          If (ifail>1 .And. ifail<=5) Then
            Write (nout,99997) w(1,2), (w(i,1),i=1,n)
          End If
        End If

        Return

99999   Format (1X,1P,3E15.3)
99998   Format (1X,F7.0,2F13.1,F13.3)
99997   Format (/,1X,'W(1,2) = ',F9.4,' W(.,1) = ',10E10.3)
99996   Format (1X,/,1X,' ** D02HBF returned with IFAIL = ',I5)
      End Subroutine ex2
    End Program d02hbfe