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

    MODULE d02tkfe_mod

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

!      .. Use Statements ..
       USE nag_library, ONLY : nag_wp
!      .. Implicit None Statement ..
       IMPLICIT NONE
!      .. Parameters ..
       INTEGER, PARAMETER                  :: mmax = 3, neq = 3, nin = 5,      &
                                              nlbc = 3, nout = 6, nrbc = 3
!      .. Local Scalars ..
       REAL (KIND=nag_wp)                  :: omega
       REAL (KIND=nag_wp)                  :: one = 1.0_nag_wp
       REAL (KIND=nag_wp)                  :: sqrofr
!      .. Local Arrays ..
       INTEGER                             :: m(neq) = (/ 1, 3, 2/)
    CONTAINS
!dec$ attributes offload : mic :: ffun
       SUBROUTINE ffun(x,y,neq,m,f)

!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. Scalar Arguments ..
          REAL (KIND=nag_wp), INTENT (IN)     :: x
          INTEGER, INTENT (IN)                :: neq
!         .. Array Arguments ..
          REAL (KIND=nag_wp), INTENT (OUT)    :: f(neq)
          REAL (KIND=nag_wp), INTENT (IN)     :: y(neq,0:*)
          INTEGER, INTENT (IN)                :: m(neq)
!         .. Executable Statements ..
          f(1) = y(2,0)
          f(2) = -(y(1,0)*y(2,2)+y(3,0)*y(3,1))*sqrofr
          f(3) = (y(2,0)*y(3,0)-y(1,0)*y(3,1))*sqrofr
          RETURN
!dec$ attributes offload : mic :: ffun
       END SUBROUTINE ffun
!dec$ attributes offload : mic :: fjac
       SUBROUTINE fjac(x,y,neq,m,dfdy)

!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. Scalar Arguments ..
          REAL (KIND=nag_wp), INTENT (IN)     :: x
          INTEGER, INTENT (IN)                :: neq
!         .. Array Arguments ..
          REAL (KIND=nag_wp), INTENT (INOUT)  :: dfdy(neq,neq,0:*)
          REAL (KIND=nag_wp), INTENT (IN)     :: y(neq,0:*)
          INTEGER, INTENT (IN)                :: m(neq)
!         .. Executable Statements ..
          dfdy(1,2,0) = one
          dfdy(2,1,0) = -y(2,2)*sqrofr
          dfdy(2,2,2) = -y(1,0)*sqrofr
          dfdy(2,3,0) = -y(3,1)*sqrofr
          dfdy(2,3,1) = -y(3,0)*sqrofr
          dfdy(3,1,0) = -y(3,1)*sqrofr
          dfdy(3,2,0) = y(3,0)*sqrofr
          dfdy(3,3,0) = y(2,0)*sqrofr
          dfdy(3,3,1) = -y(1,0)*sqrofr
          RETURN
!dec$ attributes offload : mic :: fjac
       END SUBROUTINE fjac
!dec$ attributes offload : mic :: gafun
       SUBROUTINE gafun(ya,neq,m,nlbc,ga)

!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. Scalar Arguments ..
          INTEGER, INTENT (IN)                :: neq, nlbc
!         .. Array Arguments ..
          REAL (KIND=nag_wp), INTENT (OUT)    :: ga(nlbc)
          REAL (KIND=nag_wp), INTENT (IN)     :: ya(neq,0:*)
          INTEGER, INTENT (IN)                :: m(neq)
!         .. Executable Statements ..
          ga(1) = ya(1,0)
          ga(2) = ya(2,0)
          ga(3) = ya(3,0) - omega
          RETURN
!dec$ attributes offload : mic :: gafun
       END SUBROUTINE gafun
!dec$ attributes offload : mic :: gbfun
       SUBROUTINE gbfun(yb,neq,m,nrbc,gb)

!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. Scalar Arguments ..
          INTEGER, INTENT (IN)                :: neq, nrbc
!         .. Array Arguments ..
          REAL (KIND=nag_wp), INTENT (OUT)    :: gb(nrbc)
          REAL (KIND=nag_wp), INTENT (IN)     :: yb(neq,0:*)
          INTEGER, INTENT (IN)                :: m(neq)
!         .. Executable Statements ..
          gb(1) = yb(1,0)
          gb(2) = yb(2,0)
          gb(3) = yb(3,0) + omega
          RETURN
!dec$ attributes offload : mic :: gbfun
       END SUBROUTINE gbfun
!dec$ attributes offload : mic :: gajac
       SUBROUTINE gajac(ya,neq,m,nlbc,dgady)

!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. Scalar Arguments ..
          INTEGER, INTENT (IN)                :: neq, nlbc
!         .. Array Arguments ..
          REAL (KIND=nag_wp), INTENT (INOUT)  :: dgady(nlbc,neq,0:*)
          REAL (KIND=nag_wp), INTENT (IN)     :: ya(neq,0:*)
          INTEGER, INTENT (IN)                :: m(neq)
!         .. Executable Statements ..
          dgady(1,1,0) = one
          dgady(2,2,0) = one
          dgady(3,3,0) = one
          RETURN
!dec$ attributes offload : mic :: gajac
       END SUBROUTINE gajac
!dec$ attributes offload : mic :: gbjac
       SUBROUTINE gbjac(yb,neq,m,nrbc,dgbdy)

!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. Scalar Arguments ..
          INTEGER, INTENT (IN)                :: neq, nrbc
!         .. Array Arguments ..
          REAL (KIND=nag_wp), INTENT (INOUT)  :: dgbdy(nrbc,neq,0:*)
          REAL (KIND=nag_wp), INTENT (IN)     :: yb(neq,0:*)
          INTEGER, INTENT (IN)                :: m(neq)
!         .. Executable Statements ..
          dgbdy(1,1,0) = one
          dgbdy(2,2,0) = one
          dgbdy(3,3,0) = one
          RETURN
!dec$ attributes offload : mic :: gbjac
       END SUBROUTINE gbjac
!dec$ attributes offload : mic :: guess
       SUBROUTINE guess(x,neq,m,y,dym)

!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. Scalar Arguments ..
          REAL (KIND=nag_wp), INTENT (IN)     :: x
          INTEGER, INTENT (IN)                :: neq
!         .. Array Arguments ..
          REAL (KIND=nag_wp), INTENT (OUT)    :: dym(neq)
          REAL (KIND=nag_wp), INTENT (INOUT)  :: y(neq,0:*)
          INTEGER, INTENT (IN)                :: m(neq)
!         .. Local Scalars ..
          REAL (KIND=nag_wp)                  :: five = 5.0_nag_wp
          REAL (KIND=nag_wp)                  :: half = 0.5_nag_wp
          REAL (KIND=nag_wp)                  :: ten = 10.0_nag_wp
          REAL (KIND=nag_wp)                  :: two = 2.0_nag_wp
!         .. Executable Statements ..
          y(1,0) = -(x-half)*(x*(x-one))**2
          y(2,0) = -x*(x-one)*(five*x*(x-one)+one)
          y(2,1) = -(two*x-one)*(ten*x*(x-one)+one)
          y(2,2) = -12.0_nag_wp*(five*x*(x-one)+x)
          y(3,0) = -8.0_nag_wp*omega*(x-half)**3
          y(3,1) = -24.0_nag_wp*omega*(x-half)**2
          dym(1) = y(2,0)
          dym(2) = -120.0_nag_wp*(x-half)
          dym(3) = -56.0_nag_wp*omega*(x-half)
          RETURN
!dec$ attributes offload : mic :: guess
       END SUBROUTINE guess
    END MODULE d02tkfe_mod
    PROGRAM d02tkfe

!      D02TKF Example Main Program

!      .. Use Statements ..
       USE nag_library, ONLY : d02tkf, d02tvf, d02txf, d02tyf, d02tzf
       USE d02tkfe_mod, ONLY : ffun, fjac, gafun, gajac, gbfun, gbjac, guess,  &
                               m, mmax, nag_wp, neq, nin, nlbc, nout, nrbc,    &
                               omega, one, sqrofr
!      .. Implicit None Statement ..
       IMPLICIT NONE
!      .. Local Scalars ..
       REAL (KIND=nag_wp)                  :: dx, ermx, r
       INTEGER                             :: i, iermx, ifail, ijermx, j,      &
                                              liwork, lrwork, mxmesh, ncol,    &
                                              ncont, nmesh
!      .. Local Arrays ..
       REAL (KIND=nag_wp), ALLOCATABLE     :: mesh(:), tol(:), work(:), y(:,:)
       INTEGER, ALLOCATABLE                :: ipmesh(:), iwork(:)
!      .. Intrinsic Functions ..
       INTRINSIC                              real, sqrt
!      .. Executable Statements ..
       WRITE (nout,*) 'D02TKF Example Program Results'
       WRITE (nout,*)
!      Skip heading in data file
       READ (nin,*)
       READ (nin,*) ncol, nmesh, mxmesh

       liwork = mxmesh*(11*neq+6)
       lrwork = mxmesh*(109*neq**2+78*neq+7)
       ALLOCATE (mesh(mxmesh),tol(neq),work(lrwork),y(neq,0:mmax-1), &
          ipmesh(mxmesh),iwork(liwork))

       READ (nin,*) omega
       READ (nin,*) tol(1:neq)

       dx = one/real(nmesh-1,kind=nag_wp)
       mesh(1) = 0.0_nag_wp
       DO i = 2, nmesh - 1
          mesh(i) = mesh(i-1) + dx
       END DO
       mesh(nmesh) = one

       ipmesh(1) = 1
       ipmesh(2:nmesh-1) = 2
       ipmesh(nmesh) = 1

!      Initial integrator for given problem.
       ifail = 0
       CALL d02tvf(neq,m,nlbc,nrbc,ncol,tol,mxmesh,nmesh,mesh,ipmesh,work, &
          lrwork,iwork,liwork,ifail)

!      Number of continuation steps (last r=100**ncont, sqrofr=10**ncont)
       READ (nin,*) ncont
!      Initialize problem continuation parameter.
       READ (nin,*) r
       sqrofr = sqrt(r)

CONTN: DO j = 1, ncont
          WRITE (nout,99999) tol(1), r

!         Solve problem.
          ifail = -1
          CALL d02tkf(ffun,fjac,gafun,gbfun,gajac,gbjac,guess,work,iwork, &
             ifail)

!         Extract mesh
          ifail = -1
          CALL d02tzf(mxmesh,nmesh,mesh,ipmesh,ermx,iermx,ijermx,work,iwork, &
             ifail)
          IF (ifail==1) THEN
             EXIT CONTN
          END IF

!         Print mesh and error statistics.
          WRITE (nout,99998) nmesh, ermx, iermx, ijermx
          WRITE (nout,99997) (i,ipmesh(i),mesh(i),i=1,nmesh)
!         Print solution components on mesh.
          WRITE (nout,99996)
          DO i = 1, nmesh
             ifail = 0
             CALL d02tyf(mesh(i),y,neq,mmax,work,iwork,ifail)
             WRITE (nout,99995) mesh(i), y(1:neq,0)
          END DO

          IF (j==ncont) EXIT CONTN

!         Modify continuation parameter.
          r = 100.0_nag_wp*r
          sqrofr = sqrt(r)
!         Select mesh for continuation and call continuation primer routine.
          ipmesh(2:nmesh-1) = 2
          ifail = 0
          CALL d02txf(mxmesh,nmesh,mesh,ipmesh,work,iwork,ifail)

       END DO CONTN

99999  FORMAT (/' Tolerance = ',1P,E8.1,'  R = ',E10.3)
99998  FORMAT (/' Used a mesh of ',I4,' points'/' Maximum error = ',E10.2, &
          '  in interval ',I4,' for component ',I4/)
99997  FORMAT (/' Mesh points:'/4(I4,'(',I1,')',E11.4))
99996  FORMAT (/'      x        f        f''       g')
99995  FORMAT (' ',F8.3,1X,3F9.4)
    END PROGRAM d02tkfe
