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

    Module d02rafe_mod

!     D02RAF Example Program Module:
!            Parameters and User-defined Routines

!     .. Use Statements ..
      Use nag_library, Only: nag_wp
!     .. Implicit None Statement ..
      Implicit None
!     .. Accessibility Statements ..
      Private
      Public                           :: fcn, g, jaceps, jacgep, jacobf,      &
                                          jacobg
!     .. Parameters ..
      Real (Kind=nag_wp), Parameter    :: one = 1.0_nag_wp
      Real (Kind=nag_wp), Parameter    :: two = 2.0_nag_wp
      Real (Kind=nag_wp), Parameter    :: zero = 0.0_nag_wp
      Integer, Parameter, Public       :: iset = 1, n = 3, nin = 5, nout = 6
    Contains
      Subroutine fcn(x,eps,y,f,n)

!       .. Scalar Arguments ..
        Real (Kind=nag_wp), Intent (In) :: eps, x
        Integer, Intent (In)           :: n
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (Out) :: f(n)
        Real (Kind=nag_wp), Intent (In) :: y(n)
!       .. Executable Statements ..
        f(1) = y(2)
        f(2) = y(3)
        f(3) = -y(1)*y(3) - two*(one-y(2)*y(2))*eps
        Return
      End Subroutine fcn

      Subroutine g(eps,ya,yb,bc,n)

!       .. Scalar Arguments ..
        Real (Kind=nag_wp), Intent (In) :: eps
        Integer, Intent (In)           :: n
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (Out) :: bc(n)
        Real (Kind=nag_wp), Intent (In) :: ya(n), yb(n)
!       .. Executable Statements ..
        bc(1) = ya(1)
        bc(2) = ya(2)
        bc(3) = yb(2) - one
        Return
      End Subroutine g

      Subroutine jaceps(x,eps,y,f,n)

!       .. Scalar Arguments ..
        Real (Kind=nag_wp), Intent (In) :: eps, x
        Integer, Intent (In)           :: n
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (Out) :: f(n)
        Real (Kind=nag_wp), Intent (In) :: y(n)
!       .. Executable Statements ..
        f(1:2) = zero
        f(3) = -two*(one-y(2)*y(2))
        Return
      End Subroutine jaceps

      Subroutine jacgep(eps,ya,yb,bcep,n)

!       .. Scalar Arguments ..
        Real (Kind=nag_wp), Intent (In) :: eps
        Integer, Intent (In)           :: n
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (Out) :: bcep(n)
        Real (Kind=nag_wp), Intent (In) :: ya(n), yb(n)
!       .. Executable Statements ..
        bcep(1:n) = zero
        Return
      End Subroutine jacgep

      Subroutine jacobf(x,eps,y,f,n)

!       .. Scalar Arguments ..
        Real (Kind=nag_wp), Intent (In) :: eps, x
        Integer, Intent (In)           :: n
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (Out) :: f(n,n)
        Real (Kind=nag_wp), Intent (In) :: y(n)
!       .. Executable Statements ..
        f(1:n,1:n) = zero
        f(1,2) = one
        f(2,3) = one
        f(3,1) = -y(3)
        f(3,2) = two*two*y(2)*eps
        f(3,3) = -y(1)
        Return
      End Subroutine jacobf

      Subroutine jacobg(eps,ya,yb,aj,bj,n)

!       .. Scalar Arguments ..
        Real (Kind=nag_wp), Intent (In) :: eps
        Integer, Intent (In)           :: n
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (Out) :: aj(n,n), bj(n,n)
        Real (Kind=nag_wp), Intent (In) :: ya(n), yb(n)
!       .. Executable Statements ..
        aj(1:n,1:n) = zero
        bj(1:n,1:n) = zero
        aj(1,1) = one
        aj(2,2) = one
        bj(3,2) = one
        Return
      End Subroutine jacobg
    End Module d02rafe_mod

    Program d02rafe

!     D02RAF Example Main Program

!     .. Use Statements ..
      Use d02rafe_mod, Only: fcn, g, iset, jaceps, jacgep, jacobf, jacobg, n,  &
                             nin, nout
      Use nag_library, Only: d02raf, nag_wp, x04abf
!     .. Implicit None Statement ..
      Implicit None
!     .. Local Scalars ..
      Real (Kind=nag_wp)               :: deleps, tol
      Integer                          :: ifail, ijac, init, j, ldy, liwork,   &
                                          lwork, mnp, np, numbeg, nummix,      &
                                          outchn
!     .. Local Arrays ..
      Real (Kind=nag_wp), Allocatable  :: abt(:), work(:), x(:), y(:,:)
      Integer, Allocatable             :: iwork(:)
!     .. Executable Statements ..
      Write (nout,*) 'D02RAF Example Program Results'
!     Skip heading in data file
      Read (nin,*)
      Read (nin,*) mnp, np
      ldy = n
      liwork = mnp*(2*n+1) + n
      lwork = mnp*(3*n*n+6*n+2) + 4*n*n + 3*n
      Allocate (abt(n),work(lwork),x(mnp),y(ldy,mnp),iwork(liwork))

      outchn = nout
      Write (nout,*)
      Call x04abf(iset,outchn)
      Read (nin,*) tol, deleps
      Read (nin,*) init, ijac, numbeg, nummix
      Read (nin,*) x(1), x(np)

!     ifail: behaviour on error exit
!             =1 for quiet-soft exit
!     * Set IFAIL to 111 to obtain monitoring information *
      ifail = 1
      Call d02raf(n,mnp,np,numbeg,nummix,tol,init,x,y,ldy,abt,fcn,g,ijac,      &
        jacobf,jacobg,deleps,jaceps,jacgep,work,lwork,iwork,liwork,ifail)

      If (ifail==0 .Or. ifail==4) Then
        Write (nout,*) 'Calculation using analytic Jacobians'
        If (ifail==4) Then
          Write (nout,99996) 'On exit from D02RAF IFAIL = 4'
        End If
        Write (nout,*)
        Write (nout,99999) 'Solution on final mesh of ', np, ' points'
        Write (nout,*) '      X(I)        Y1(I)        Y2(I)        Y3(I)'
        Write (nout,99998)(x(j),y(1:n,j),j=1,np)
        Write (nout,*)
        Write (nout,*) 'Maximum estimated error by components'
        Write (nout,99997) abt(1:n)
      Else
        Write (nout,99996) ' ** D02RAF returned with IFAIL = ', ifail
      End If

99999 Format (1X,A,I2,A)
99998 Format (1X,F10.6,3F13.4)
99997 Format (11X,1P,3E13.2)
99996 Format (1X,A,I5)
    End Program d02rafe