!   D03RBF Example Program Text
!   Mark 25 Release. NAG Copyright 2014.

    Module d03rbfe_mod

!     D03RBF 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                               :: bndary, inidom, monitr, pdedef,  &
                                              pdeiv
!     .. Parameters ..
      Real (Kind=nag_wp), Parameter, Public :: one = 1.0_nag_wp
      Real (Kind=nag_wp), Parameter, Public :: zero = 0.0_nag_wp
      Real (Kind=nag_wp), Parameter, Public :: twant(2) = (/0.25_nag_wp,one/)
      Integer, Parameter, Public           :: itrace = 0, nin = 5, nout = 6,   &
                                              npde = 2
!     .. Local Scalars ..
      Integer, Public, Save                :: iout
    Contains
      Subroutine pdeiv(npts,npde,t,x,y,u)

!       .. Parameters ..
        Real (Kind=nag_wp), Parameter        :: eps = 0.001_nag_wp
!       .. Scalar Arguments ..
        Real (Kind=nag_wp), Intent (In)      :: t
        Integer, Intent (In)                 :: npde, npts
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (Out)     :: u(npts,npde)
        Real (Kind=nag_wp), Intent (In)      :: x(npts), y(npts)
!       .. Local Scalars ..
        Real (Kind=nag_wp)                   :: a
        Integer                              :: i
!       .. Intrinsic Procedures ..
        Intrinsic                            :: exp
!       .. Executable Statements ..
        Do i = 1, npts
          a = (4.0_nag_wp*(y(i)-x(i))-t)/(32.0_nag_wp*eps)
          If (a<=zero) Then
            u(i,1) = 0.75_nag_wp - 0.25_nag_wp/(one+exp(a))
            u(i,2) = 0.75_nag_wp + 0.25_nag_wp/(one+exp(a))
          Else
            a = -a
            u(i,1) = 0.75_nag_wp - 0.25_nag_wp*exp(a)/(one+exp(a))
            u(i,2) = 0.75_nag_wp + 0.25_nag_wp*exp(a)/(one+exp(a))
          End If
        End Do

        Return
      End Subroutine pdeiv
      Subroutine inidom(maxpts,xmin,xmax,ymin,ymax,nx,ny,npts,nrows,nbnds, &
        nbpts,lrow,irow,icol,llbnd,ilbnd,lbnd,ierr)

!       .. Use Statements ..
        Use nag_library, Only: d03ryf
!       .. Scalar Arguments ..
        Real (Kind=nag_wp), Intent (Out)     :: xmax, xmin, ymax, ymin
        Integer, Intent (Inout)              :: ierr
        Integer, Intent (In)                 :: maxpts
        Integer, Intent (Out)                :: nbnds, nbpts, npts, nrows, nx, &
                                                ny
!       .. Array Arguments ..
        Integer, Intent (Inout)              :: icol(*), ilbnd(*), irow(*),    &
                                                lbnd(*), llbnd(*), lrow(*)
!       .. Local Scalars ..
        Integer                              :: i, ifail, j, leniwk
!       .. Local Arrays ..
        Integer                              :: icold(105), ilbndd(28),        &
                                                irowd(11), iwk(122),           &
                                                lbndd(72), llbndd(28), lrowd(11)
        Character (33)                       :: pgrid(11)
!       .. Executable Statements ..
        icold(1:105) = (/0,1,2,0,1,2,3,4,5,6,7,8,9,10,0,1,2,3,4,5,6,7,8,9,10, &
          0,1,2,3,4,5,6,7,8,9,10,0,1,2,3,4,5,8,9,10,0,1,2,3,4,5,6,7,8,9,10,0, &
          1,2,3,4,5,6,7,8,9,10,0,1,2,3,4,5,6,7,8,9,10,0,1,2,3,4,5,6,7,8,0,1,2, &
          3,4,5,6,7,8,0,1,2,3,4,5,6,7,8/)
        ilbndd(1:28) = (/1,2,3,4,1,4,1,2,3,4,3,4,1,2,12,23,34,41,14,41,12,23, &
          34,41,43,14,21,32/)
        irowd(1:11) = (/0,1,2,3,4,5,6,7,8,9,10/)

        lbndd(1:72) = (/2,4,15,26,37,46,57,68,79,88,98,99,100,101,102,103,104, &
          96,86,85,84,83,82,70,59,48,39,28,17,6,8,9,10,11,12,13,18,29,40,49, &
          60,72,73,74,75,76,77,67,56,45,36,25,33,32,42,52,53,43,1,97,105,87, &
          81,3,7,71,78,14,31,51,54,34/)

        llbndd(1:28) = (/1,2,11,18,19,24,31,37,42,48,53,55,56,58,59,60,61,62, &
          63,64,65,66,67,68,69,70,71,72/)
        lrowd(1:11) = (/1,4,15,26,37,46,57,68,79,88,97/)

        nx = 11
        ny = 11

!       Check MAXPTS against rough estimate of NPTS

        npts = nx*ny
        If (maxpts<npts) Then
          ierr = -1
          Return
        End If

        xmin = zero
        ymin = zero
        xmax = one
        ymax = one

        nrows = 11
        npts = 105
        nbnds = 28
        nbpts = 72

        Do i = 1, nrows
          lrow(i) = lrowd(i)
          irow(i) = irowd(i)
        End Do

        Do i = 1, nbnds
          llbnd(i) = llbndd(i)
          ilbnd(i) = ilbndd(i)
        End Do

        Do i = 1, nbpts
          lbnd(i) = lbndd(i)
        End Do

        Do i = 1, npts
          icol(i) = icold(i)
        End Do

        Write (nout,*) 'Base grid:'
        Write (nout,*)
        leniwk = 122
        ifail = -1

        Call d03ryf(nx,ny,npts,nrows,nbnds,nbpts,lrow,irow,icol,llbnd,ilbnd, &
          lbnd,iwk,leniwk,pgrid,ifail)

        If (ifail==0) Then
          Write (nout,*) ' '
          Do j = 1, ny
            Write (nout,*) pgrid(j)
            Write (nout,*) ' '
          End Do
          Write (nout,*) ' '
        End If

        Return
      End Subroutine inidom
      Subroutine pdedef(npts,npde,t,x,y,u,ut,ux,uy,uxx,uxy,uyy,res)

!       .. Parameters ..
        Real (Kind=nag_wp), Parameter        :: eps = 1E-3_nag_wp
!       .. Scalar Arguments ..
        Real (Kind=nag_wp), Intent (In)      :: t
        Integer, Intent (In)                 :: npde, npts
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (Out)     :: res(npts,npde)
        Real (Kind=nag_wp), Intent (In)      :: u(npts,npde), ut(npts,npde),   &
                                                ux(npts,npde), uxx(npts,npde), &
                                                uxy(npts,npde), uy(npts,npde), &
                                                uyy(npts,npde), x(npts), y(npts)
!       .. Local Scalars ..
        Integer                              :: i
!       .. Executable Statements ..
        !$Omp Do
        Do i = 1, npts
          res(i,1) = -u(i,1)*ux(i,1) - u(i,2)*uy(i,1) + &
            eps*(uxx(i,1)+uyy(i,1))
          res(i,2) = -u(i,1)*ux(i,2) - u(i,2)*uy(i,2) + &
            eps*(uxx(i,2)+uyy(i,2))
          res(i,1) = ut(i,1) - res(i,1)
          res(i,2) = ut(i,2) - res(i,2)
        End Do
        !$Omp End Do

        Return
      End Subroutine pdedef
      Subroutine bndary(npts,npde,t,x,y,u,ut,ux,uy,nbnds,nbpts,llbnd,ilbnd, &
        lbnd,res)

!       .. Parameters ..
        Real (Kind=nag_wp), Parameter        :: eps = 1E-3_nag_wp
!       .. Scalar Arguments ..
        Real (Kind=nag_wp), Intent (In)      :: t
        Integer, Intent (In)                 :: nbnds, nbpts, npde, npts
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (Inout)   :: res(npts,npde)
        Real (Kind=nag_wp), Intent (In)      :: u(npts,npde), ut(npts,npde),   &
                                                ux(npts,npde), uy(npts,npde),  &
                                                x(npts), y(npts)
        Integer, Intent (In)                 :: ilbnd(nbnds), lbnd(nbpts),     &
                                                llbnd(nbnds)
!       .. Local Scalars ..
        Real (Kind=nag_wp)                   :: a
        Integer                              :: i, k
!       .. Intrinsic Procedures ..
        Intrinsic                            :: exp
!       .. Executable Statements ..
        Do k = llbnd(1), nbpts
          i = lbnd(k)
          a = (-4.0_nag_wp*x(i)+4.0_nag_wp*y(i)-t)/(32.0_nag_wp*eps)
          If (a<=zero) Then
            res(i,1) = 0.75_nag_wp - 0.25_nag_wp/(one+exp(a))
            res(i,2) = 0.75_nag_wp + 0.25_nag_wp/(one+exp(a))
          Else
            a = -a
            res(i,1) = 0.75_nag_wp - 0.25_nag_wp*exp(a)/(one+exp(a))
            res(i,2) = 0.75_nag_wp + 0.25_nag_wp*exp(a)/(one+exp(a))
          End If
          res(i,1:2) = u(i,1:2) - res(i,1:2)
        End Do

        Return
      End Subroutine bndary
      Subroutine monitr(npde,t,dt,dtnew,tlast,nlev,xmin,ymin,dxb,dyb,lgrid, &
        istruc,lsol,sol,ierr)

!       .. Use Statements ..
        Use nag_library, Only: d03rzf
!       .. Parameters ..
        Integer, Parameter                   :: maxpts = 2500, nout = 6
!       .. Scalar Arguments ..
        Real (Kind=nag_wp), Intent (In)      :: dt, dtnew, dxb, dyb, t, xmin,  &
                                                ymin
        Integer, Intent (Inout)              :: ierr
        Integer, Intent (In)                 :: nlev, npde
        Logical, Intent (In)                 :: tlast
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (In)      :: sol(*)
        Integer, Intent (In)                 :: istruc(*), lgrid(*), lsol(nlev)
!       .. Local Scalars ..
        Integer                              :: i, ifail, ipsol, level, npts
!       .. Local Arrays ..
        Real (Kind=nag_wp)                   :: uex(105,2), x(maxpts), y(maxpts)
!       .. Executable Statements ..
        ifail = -1
levels: Do level = 1, nlev
          If (.Not. tlast) Exit levels
          ipsol = lsol(level)

!         Get grid information

          Call d03rzf(level,nlev,xmin,ymin,dxb,dyb,lgrid,istruc,npts,x,y, &
            maxpts,ifail)
          If (ifail/=0) Then
            ierr = 1
            Exit levels
          End If

!         Skip printing if iout<2 or level>1.
          If (iout/=2 .Or. level/=1) Cycle levels

!         Get exact solution
          Call pdeiv(npts,npde,t,x,y,uex)
          Write (nout,*)
          Write (nout,99999) t
          Write (nout,*)
          Write (nout,99998) 'x', 'y', 'approx u', 'exact u', 'approx v', &
            'exact v'
          Write (nout,*)
          ipsol = lsol(level)
          Do i = 1, npts, 2
            Write (nout,99997) x(i), y(i), sol(ipsol+i), uex(i,1), &
              sol(ipsol+npts+i), uex(i,2)
          End Do
          Write (nout,*)

        End Do levels

        Return
99999   Format (' Solution at every 2nd grid point in level 1 at time ',F8.4, &
          ':')
99998   Format (7X,A,10X,A,8X,A,5X,A,4X,A,4X,A)
99997   Format (6(1X,E11.4))
      End Subroutine monitr
    End Module d03rbfe_mod
    Program d03rbfe

!     D03RBF Example Main Program

!     .. Use Statements ..
      Use nag_library, Only: d03rbf, nag_wp
      Use d03rbfe_mod, Only: bndary, inidom, iout, itrace, monitr, nin, nout,  &
                             npde, one, pdedef, pdeiv, twant, zero
!     .. Implicit None Statement ..
      Implicit None
!     .. Local Scalars ..
      Real (Kind=nag_wp)                   :: tols, tolt, tout, ts
      Integer                              :: i, ifail, ind, j, leniwk,        &
                                              lenlwk, lenrwk, maxlev, mxlev,   &
                                              npts
!     .. Local Arrays ..
      Real (Kind=nag_wp)                   :: dt(3)
      Real (Kind=nag_wp), Allocatable      :: optr(:,:), rwk(:)
      Integer, Allocatable                 :: iwk(:)
      Integer                              :: opti(4)
      Logical, Allocatable                 :: lwk(:)
!     .. Executable Statements ..
      Write (nout,*) 'D03RBF Example Program Results'
!     Skip heading in data file
      Read (nin,*)
      Read (nin,*) npts, mxlev

      leniwk = npts*(5*mxlev+14) + 2 + 7*mxlev
      lenlwk = 2*npts
      lenrwk = npts*npde*(5*mxlev+9+18*npde) + 2*npts
      Allocate (rwk(lenrwk),iwk(leniwk),lwk(lenlwk),optr(3,npde))


!     Specify that we are starting the integration in time (ind = 0 normally).
!     Note: we have parallelised the loop in the function pdedef using OpenMP
!     so set the alternative value of ind to indicate that this can be run in
!     parallel if we are using a multithreaded implementation. Either option
!     is OK for serial NAG Library implementations from Mark 25 onwards.
      ind = 10

      ts = zero
      dt(1) = 0.001_nag_wp
      dt(2) = 1.0E-7_nag_wp
      dt(3) = zero
      tols = 0.1_nag_wp
      tolt = 0.05_nag_wp
      opti(1) = 5
      maxlev = opti(1)
      opti(2:4) = 0
      optr(1:3,1:npde) = one

!     Call main routine
      Do iout = 1, 2
        tout = twant(iout)

!       ifail: behaviour on error exit   
!              =0 for hard exit, =1 for quiet-soft, =-1 for noisy-soft  
        ifail = 0
        Call d03rbf(npde,ts,tout,dt,tols,tolt,inidom,pdedef,bndary,pdeiv, &
          monitr,opti,optr,rwk,lenrwk,iwk,leniwk,lwk,lenlwk,itrace,ind,ifail)

!       Print statistics

        Write (nout,99999) 'Statistics:'
        Write (nout,99998) 'Time = ', ts
        Write (nout,99997) 'Total number of accepted timesteps =', iwk(1)
        Write (nout,99997) 'Total number of rejected timesteps =', iwk(2)
        Write (nout,99996)
        maxlev = opti(1)
        Do j = 1, maxlev
          If (iwk(j+2)/=0) Then
            Write (nout,'(I6,4I10)') j, (iwk(j+2+i*maxlev),i=0,3)
          End If
        End Do
        Write (nout,99995)
        Do j = 1, maxlev
          If (iwk(j+2)/=0) Then
            Write (nout,'(I6,2I14)') j, (iwk(j+2+i*maxlev),i=4,5)
          End If
        End Do
        Write (nout,*)

      End Do

99999 Format (1X,A)
99998 Format (1X,A,F8.4)
99997 Format (1X,A,I5)
99996 Format (1X/'             T o t a l  n u m b e r  o f    '/'   ', &
        '       Residual  Jacobian    Newton   Lin sys'/'       ', &
        '      evals     evals     iters     iters'/' At level ')
99995 Format (1X/'             M a x i m u m  n u m b e r  o f'/'   ', &
        '           Newton iters    Lin sys iters '/' At level ')
    End Program d03rbfe