!   D02NEF Example Program Text
!   Mark 24 Release. NAG Copyright 2012.

    Module d02nefe_mod

!     D02NEF 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        :: alpha = 0.04_nag_wp
      Real (Kind=nag_wp), Parameter        :: beta = 1.0E4_nag_wp
      Real (Kind=nag_wp), Parameter        :: gamma = 3.0E7_nag_wp
      Real (Kind=nag_wp), Parameter        :: one = 1.0_nag_wp
      Real (Kind=nag_wp), Parameter        :: two = 2.0_nag_wp
      Integer, Parameter                   :: ml = 1, mu = 2, neq1 = 3,        &
                                              neq2 = 1, nin = 5, nout = 6
    Contains
      Subroutine myjac1(neq,ml,mu,t,y,ydot,pd,cj)

!       .. Scalar Arguments ..
        Real (Kind=nag_wp), Intent (In)      :: cj, t
        Integer, Intent (In)                 :: ml, mu, neq
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (Out)     :: pd(2*ml+mu+1,neq)
        Real (Kind=nag_wp), Intent (In)      :: y(neq), ydot(neq)
!       .. Local Scalars ..
        Integer                              :: md, ms
!       .. Executable Statements ..
!       Main diagonal pdfull(i,i), i=1,neq
        md = mu + ml + 1
        pd(md,1) = -alpha - cj
        pd(md,2) = -beta*y(3) - two*gamma*y(2) - cj
        pd(md,3) = -cj
!       1 Sub-diagonal pdfull(i+1:i), i=1,neq-1
        ms = md + 1
        pd(ms,1) = alpha
        pd(ms,2) = two*gamma*y(2)
!       First super-diagonal pdfull(i-1,i), i=2, neq
        ms = md - 1
        pd(ms,2) = beta*y(3)
        pd(ms,3) = -beta*y(2)
!       Second super-diagonal pdfull(i-2,i), i=3, neq
        ms = md - 2
        pd(ms,3) = beta*y(2)

        Return
      End Subroutine myjac1
      Subroutine myjac2(neq,t,y,ydot,pd,cj)

!       .. Scalar Arguments ..
        Real (Kind=nag_wp), Intent (In)      :: cj, t
        Integer, Intent (In)                 :: neq
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (Out)     :: pd(neq*neq)
        Real (Kind=nag_wp), Intent (In)      :: y(neq), ydot(neq)
!       .. Intrinsic Procedures ..
        Intrinsic                            :: exp
!       .. Executable Statements ..
        pd(1) = -two*y(1) + 0.1E0_nag_wp*t*y(1)*exp(y(1))

        Return
      End Subroutine myjac2

      Subroutine res1(neq,t,y,ydot,r,ires,iuser,ruser)

!       .. Scalar Arguments ..
        Real (Kind=nag_wp), Intent (In)      :: t
        Integer, Intent (Inout)              :: ires
        Integer, Intent (In)                 :: neq
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (Out)     :: r(neq)
        Real (Kind=nag_wp), Intent (Inout)   :: ruser(*)
        Real (Kind=nag_wp), Intent (In)      :: y(neq), ydot(neq)
        Integer, Intent (Inout)              :: iuser(*)
!       .. Executable Statements ..
        r(1) = -alpha*y(1) + beta*y(2)*y(3) - ydot(1)
        r(2) = alpha*y(1) - beta*y(2)*y(3) - gamma*y(2)*y(2) - ydot(2)
        r(3) = gamma*y(2)*y(2) - ydot(3)
        Return
      End Subroutine res1
      Subroutine jac1(neq,t,y,ydot,pd,cj,iuser,ruser)

!       .. Use Statements ..
        Use nag_library, Only: d02nez
!       .. Scalar Arguments ..
        Real (Kind=nag_wp), Intent (In)      :: cj, t
        Integer, Intent (In)                 :: neq
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (Inout)   :: pd(*), ruser(*)
        Real (Kind=nag_wp), Intent (In)      :: y(neq), ydot(neq)
        Integer, Intent (Inout)              :: iuser(*)
!       .. Local Scalars ..
        Integer                              :: ijac, ml, mu
!       .. Executable Statements ..
        ml = iuser(1)
        mu = iuser(2)
        ijac = iuser(3)

        If (ijac==1) Then
          Call myjac1(neq,ml,mu,t,y,ydot,pd,cj)
        Else
          Call d02nez(neq,t,y,ydot,pd,cj,iuser,ruser)
        End If

        Return
      End Subroutine jac1
      Subroutine res2(neq,t,y,ydot,r,ires,iuser,ruser)

!       .. Scalar Arguments ..
        Real (Kind=nag_wp), Intent (In)      :: t
        Integer, Intent (Inout)              :: ires
        Integer, Intent (In)                 :: neq
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (Out)     :: r(neq)
        Real (Kind=nag_wp), Intent (Inout)   :: ruser(*)
        Real (Kind=nag_wp), Intent (In)      :: y(neq), ydot(neq)
        Integer, Intent (Inout)              :: iuser(*)
!       .. Intrinsic Procedures ..
        Intrinsic                            :: exp
!       .. Executable Statements ..
        r(1) = 4.0_nag_wp - y(1)**2 + t*0.1E0_nag_wp*exp(y(1))
        Return
      End Subroutine res2
      Subroutine jac2(neq,t,y,ydot,pd,cj,iuser,ruser)

!       .. Use Statements ..
        Use nag_library, Only: d02nez
!       .. Scalar Arguments ..
        Real (Kind=nag_wp), Intent (In)      :: cj, t
        Integer, Intent (In)                 :: neq
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (Inout)   :: pd(*), ruser(*)
        Real (Kind=nag_wp), Intent (In)      :: y(neq), ydot(neq)
        Integer, Intent (Inout)              :: iuser(*)
!       .. Local Scalars ..
        Integer                              :: ijac
!       .. Executable Statements ..
        ijac = iuser(1)

        If (ijac==1) Then
          Call myjac2(neq,t,y,ydot,pd,cj)
        Else
          Call d02nez(neq,t,y,ydot,pd,cj,iuser,ruser)
        End If

        Return
      End Subroutine jac2
    End Module d02nefe_mod
    Program d02nefe

!     D02NEF Example Main Program

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

      Call ex1

      Call ex2

    Contains
      Subroutine ex1

!       .. Use Statements ..
        Use nag_library, Only: d02mcf, d02mwf, d02nef, d02npf, nag_wp
        Use d02nefe_mod, Only: jac1, ml, mu, neq1, nin, res1
!       .. Local Scalars ..
        Real (Kind=nag_wp)                   :: h0, hmax, t, tout
        Integer                              :: i, ifail, ijac, itask, itol,   &
                                                j, lcom, licom, maxord, neq
        Character (8)                        :: jceval
!       .. Local Arrays ..
        Real (Kind=nag_wp), Allocatable      :: atol(:), com(:), rtol(:),      &
                                                y(:), ydot(:)
        Real (Kind=nag_wp)                   :: ruser(1)
        Integer, Allocatable                 :: icom(:)
        Integer                              :: iuser(3)
!       .. Executable Statements ..
        Write (nout,*)
        Write (nout,*) 'D02NEF Example 1'
!       Skip heading in data file
        Read (nin,*)
        Read (nin,*) maxord
        neq = neq1
        lcom = 40 + (maxord+4)*neq + (2*ml+mu+1)*neq + 2*(neq/(ml+mu+1)+1)
        licom = 50 + neq
        Allocate (atol(neq),com(lcom),rtol(neq),y(neq),ydot(neq),icom(licom))
        Read (nin,*) ijac, itol
        Read (nin,*) rtol(1:neq)
        Read (nin,*) atol(1:neq)
        Read (nin,*) ydot(1:neq)
        If (ijac==1) Then
          jceval = 'Analytic'
        Else
          jceval = 'Numeric'
        End If
!       Set initial values
        Read (nin,*) y(1:neq)

!       Initialize the problem, specifying that the Jacobian is to be
!       evaluated analytically using the provided routine jac.

        Read (nin,*) hmax, h0
        Read (nin,*) t, tout

!       ifail: behaviour on error exit   
!              =0 for hard exit, =1 for quiet-soft, =-1 for noisy-soft
        ifail = 0
        Call d02mwf(neq,maxord,jceval,hmax,h0,itol,icom,licom,com,lcom,ifail)

!       Specify that the Jacobian is banded.

        ifail = 0
        Call d02npf(neq,ml,mu,icom,licom,ifail)

!       Use the iuser array to pass the band dimensions through to jac.
!       An alternative would be to hard code values for ml and mu in jac.

        iuser(1) = ml
        iuser(2) = mu
        iuser(3) = ijac

        Write (nout,99999)(i,i=1,neq)
        Write (nout,99998) t, (y(i),i=1,neq)
        itask = 0

!       Obtain the solution at 5 equally spaced values of T.

loop:   Do j = 1, 5
          ifail = -1
          Call d02nef(neq,t,tout,y,ydot,rtol,atol,itask,res1,jac1,icom,com, &
            lcom,iuser,ruser,ifail)
          Write (nout,99998) t, (y(i),i=1,neq)
          If (ifail/=0) Then
            Write (nout,99997) ifail
            Exit loop
          End If
          tout = tout + 0.02_nag_wp
          Call d02mcf(icom)
        End Do loop

        Write (nout,*)
        Write (nout,99996) itask

99999   Format (/1X,'    t ',5X,3('      Y(',I1,')  '))
99998   Format (1X,F8.4,3X,3(F12.6))
99997   Format (1X,' ** D02NEF returned with IFAIL = ',I5)
99996   Format (1X,'The integrator completed task, ITASK = ',I4)
      End Subroutine ex1
      Subroutine ex2

!       .. Use Statements ..
        Use nag_library, Only: d02mcf, d02mwf, d02nef, nag_wp
        Use d02nefe_mod, Only: jac2, neq2, nin, res2
!       .. Local Scalars ..
        Real (Kind=nag_wp)                   :: h0, hmax, t, tout
        Integer                              :: i, ifail, ijac, itask, itol,   &
                                                j, lcom, licom, maxord, neq
        Character (8)                        :: jceval
!       .. Local Arrays ..
        Real (Kind=nag_wp), Allocatable      :: atol(:), com(:), rtol(:),      &
                                                y(:), ydot(:)
        Real (Kind=nag_wp)                   :: ruser(1)
        Integer, Allocatable                 :: icom(:)
        Integer                              :: iuser(1)
!       .. Executable Statements ..
        Write (nout,*)
        Write (nout,*) 'D02NEF Example 2'
        Write (nout,*)
        Read (nin,*)
        Read (nin,*) maxord
        neq = neq2
        lcom = 40 + (maxord+4)*neq + neq*neq
        licom = 50 + neq
        Allocate (atol(neq),com(lcom),rtol(neq),y(neq),ydot(neq),icom(licom))
        Read (nin,*) ijac, itol
        Read (nin,*) rtol(1:neq)
        Read (nin,*) atol(1:neq)
        Read (nin,*) ydot(1:neq)
        If (ijac==1) Then
          jceval = 'Analytic'
        Else
          jceval = 'Numeric'
        End If

!       Initialize the problem, specifying that the Jacobian is to be
!       evaluated analytically using the provided routine jac.

        Read (nin,*) y(1:neq)
        Read (nin,*) hmax, h0
        Read (nin,*) t, tout

        ifail = 0
        Call d02mwf(neq,maxord,jceval,hmax,h0,itol,icom,licom,com,lcom,ifail)

!       Use the iuser array to pass whether numerical or analytic Jacobian
!       is to be used.

        iuser(1) = ijac

        Write (nout,99999)(i,i=1,neq)
        Write (nout,99998) t, y(1:neq)
        itask = 0

!       Obtain the solution at 5 equally spaced values of t.

loop:   Do j = 1, 5

          ifail = -1
          Call d02nef(neq,t,tout,y,ydot,rtol,atol,itask,res2,jac2,icom,com, &
            lcom,iuser,ruser,ifail)

          Write (nout,99998) t, y(1:neq)
          If (ifail/=0) Then
            Write (nout,99997) ifail
            Exit loop
          End If
          tout = tout + 0.2_nag_wp
          Call d02mcf(icom)
        End Do loop

        Write (nout,*)
        Write (nout,99996) itask

99999   Format (/1X,'    t            y(',I1,')')
99998   Format (1X,F8.4,3X,3(F12.6))
99997   Format (1X,' ** D02NEF returned with IFAIL = ',I5)
99996   Format (1X,'The integrator completed task, ITASK = ',I4)
      End Subroutine ex2
    End Program d02nefe
