!   D02RAF Example Program Text
!   Mark 23 Release. NAG Copyright 2011.

    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
!      .. 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                  :: iset = 1, n = 3, nin = 5, nout = 6
    CONTAINS
       SUBROUTINE fcn(x,eps,y,f,n)

!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. 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)

!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. 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)

!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. 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)

!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. 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)

!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. 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)

!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. 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 nag_library, ONLY : d02raf, x04abf
       USE d02rafe_mod, ONLY : fcn, g, iset, jaceps, jacgep, jacobf, jacobg,   &
                               n, nag_wp, nin, nout
!      .. 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) WRITE (nout,99996) 'On exit from D02RAF IFAIL = 4'
          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.3,3F13.4)
99997  FORMAT (11X,1P,3E13.2)
99996  FORMAT (1X,A,I5)
    END PROGRAM d02rafe