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

    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
!      .. Parameters ..
       REAL (KIND=nag_wp), PARAMETER       :: one = 1.0_nag_wp
       REAL (KIND=nag_wp), PARAMETER       :: zero = 0.0_nag_wp
       REAL (KIND=nag_wp), PARAMETER       :: twant(2) = (/ 0.25_nag_wp, one/)
       INTEGER, PARAMETER                  :: itrace = 0, nin = 5, nout = 6,   &
                                              npde = 2
!      .. Local Scalars ..
       INTEGER                             :: iout
    CONTAINS
       SUBROUTINE pdeiv(npts,npde,t,x,y,u)

!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. 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 Functions ..
          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
!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. 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)

!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. 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 ..
          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))
          END DO
          res(1:npts,1:2) = ut(1:npts,1:2) - res(1:npts,1:2)

          RETURN
       END SUBROUTINE pdedef
       SUBROUTINE bndary(npts,npde,t,x,y,u,ut,ux,uy,nbnds,nbpts,llbnd,ilbnd, &
             lbnd,res)

!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. 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 Functions ..
          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
!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. 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
       USE d03rbfe_mod, ONLY : bndary, inidom, iout, itrace, monitr, nag_wp,   &
                               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))

       ind = 0
       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