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

    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
!dec$ attributes offload : mic :: myjac1
       SUBROUTINE myjac1(neq,ml,mu,t,y,ydot,pd,cj)

!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. 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
!dec$ attributes offload : mic :: myjac1
       END SUBROUTINE myjac1
!dec$ attributes offload : mic :: myjac2
       SUBROUTINE myjac2(neq,t,y,ydot,pd,cj)

!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. 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 Functions ..
          INTRINSIC                              exp
!         .. Executable Statements ..
          pd(1) = -two*y(1) + 0.1E0_nag_wp*t*y(1)*exp(y(1))

          RETURN
!dec$ attributes offload : mic :: myjac2
       END SUBROUTINE myjac2

!dec$ attributes offload : mic :: res1
       SUBROUTINE res1(neq,t,y,ydot,r,ires,iuser,ruser)

!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. 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
!dec$ attributes offload : mic :: res1
       END SUBROUTINE res1
!dec$ attributes offload : mic :: jac1
       SUBROUTINE jac1(neq,t,y,ydot,pd,cj,iuser,ruser)

!         .. Use Statements ..
          USE nag_library, ONLY : d02nez
!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. 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
!dec$ attributes offload : mic :: jac1
       END SUBROUTINE jac1
!dec$ attributes offload : mic :: res2
       SUBROUTINE res2(neq,t,y,ydot,r,ires,iuser,ruser)

!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. 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 Functions ..
          INTRINSIC                              exp
!         .. Executable Statements ..
          r(1) = 4.0_nag_wp - y(1)**2 + t*0.1E0_nag_wp*exp(y(1))
          RETURN
!dec$ attributes offload : mic :: res2
       END SUBROUTINE res2
!dec$ attributes offload : mic :: jac2
       SUBROUTINE jac2(neq,t,y,ydot,pd,cj,iuser,ruser)

!         .. Use Statements ..
          USE nag_library, ONLY : d02nez
!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. 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
!dec$ attributes offload : mic :: jac2
       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
!dec$ attributes offload : mic :: ex1
       SUBROUTINE ex1

!         .. Use Statements ..
          USE nag_library, ONLY : d02mcf, d02mwf, d02nef, d02npf, nag_wp
          USE d02nefe_mod, ONLY : jac1, ml, mu, neq1, nin, res1
!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. 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)
!dec$ attributes offload : mic :: ex1
       END SUBROUTINE ex1
!dec$ attributes offload : mic :: ex2
       SUBROUTINE ex2

!         .. Use Statements ..
          USE nag_library, ONLY : d02mcf, d02mwf, d02nef, nag_wp
          USE d02nefe_mod, ONLY : jac2, neq2, nin, res2
!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. 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)
!dec$ attributes offload : mic :: ex2
       END SUBROUTINE ex2
    END PROGRAM d02nefe
