!   D02TKF Example Program Text
!   Mark 26 Release. NAG Copyright 2016.

    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
!     .. Accessibility Statements ..
      Private
      Public                           :: ffun, fjac, gafun, gajac, gbfun,     &
                                          gbjac, guess
!     .. Parameters ..
      Integer, Parameter, Public       :: mmax = 3, neq = 3, nin = 5,          &
                                          nlbc = 3, nout = 6, nrbc = 3
!     .. Local Scalars ..
      Real (Kind=nag_wp), Public, Save :: omega
      Real (Kind=nag_wp), Public, Save :: one = 1.0_nag_wp
      Real (Kind=nag_wp), Public, Save :: sqrofr
!     .. Local Arrays ..
      Integer, Public, Save            :: m(neq) = (/1,3,2/)
    Contains
      Subroutine ffun(x,y,neq,m,f)

!       .. 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
      End Subroutine ffun
      Subroutine fjac(x,y,neq,m,dfdy)

!       .. 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
      End Subroutine fjac
      Subroutine gafun(ya,neq,m,nlbc,ga)

!       .. 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
      End Subroutine gafun
      Subroutine gbfun(yb,neq,m,nrbc,gb)

!       .. 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
      End Subroutine gbfun
      Subroutine gajac(ya,neq,m,nlbc,dgady)

!       .. 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
      End Subroutine gajac
      Subroutine gbjac(yb,neq,m,nrbc,dgbdy)

!       .. 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
      End Subroutine gbjac
      Subroutine guess(x,neq,m,y,dym)

!       .. 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)
!       .. Executable Statements ..
        y(1,0) = -(x-0.5_nag_wp)*(x*(x-one))**2
        y(2,0) = -x*(x-one)*(5._nag_wp*x*(x-one)+one)
        y(2,1) = -(2._nag_wp*x-one)*(10._nag_wp*x*(x-one)+one)
        y(2,2) = -12.0_nag_wp*(5._nag_wp*x*(x-one)+x)
        y(3,0) = -8.0_nag_wp*omega*(x-0.5_nag_wp)**3
        y(3,1) = -24.0_nag_wp*omega*(x-0.5_nag_wp)**2
        dym(1) = y(2,0)
        dym(2) = -120.0_nag_wp*(x-0.5_nag_wp)
        dym(3) = -56.0_nag_wp*omega*(x-0.5_nag_wp)
        Return
      End Subroutine guess
    End Module d02tkfe_mod
    Program d02tkfe

!     D02TKF Example Main Program

!     .. Use Statements ..
      Use nag_library, Only: d02tkf, d02tvf, d02txf, d02tyf, d02tzf, nag_wp
      Use d02tkfe_mod, Only: ffun, fjac, gafun, gajac, gbfun, gbjac, guess, m, &
                             mmax, 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, licomm,  &
                                          lrcomm, mxmesh, ncol, ncont, nmesh
!     .. Local Arrays ..
      Real (Kind=nag_wp), Allocatable  :: mesh(:), tol(:), work(:), y(:,:)
      Integer, Allocatable             :: icomm(:), ipmesh(:)
!     .. Intrinsic Procedures ..
      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

      licomm = 23 + neq + mxmesh
      lrcomm = mxmesh*(109*neq**2+78*neq+7)
      Allocate (mesh(mxmesh),tol(neq),work(lrcomm),y(neq,0:mmax-1),            &
        ipmesh(mxmesh),icomm(licomm))

      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,      &
        lrcomm,icomm,licomm,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,icomm,ifail)

!       Extract mesh
        ifail = -1
        Call d02tzf(mxmesh,nmesh,mesh,ipmesh,ermx,iermx,ijermx,work,icomm,     &
          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,icomm,ifail)
          Write (nout,99995) mesh(i), y(1:neq,0)
        End Do

        If (j==ncont) Then
          Exit contn
        End If

!       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,icomm,ifail)

      End Do contn

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