c This is an example driver for running BACOLRI-3 (referred to in this
c file simply as BACOLRI. See comments within bacolri.f for details 
c regarding its use.

c This example hard codes parameters npde, xa, xb, tstart, nintinit, 
c nintmx, kcol, nderiv, est and nxoutmx; edit these as needed.
c
c The user is prompted for an output time, a value to use for relative
c and absolute error tolerances, and the number of points to provide
c output (uniformly spaced in [xa, xb]).

c Output is printed as columns to the standard output unit and
c optionally, interpolant data is written to the file 'Bsplines'.
c
c When prompted to input a tolerance, it is recommended to use
c engineering notation (e.g., '1e-6'.)
c
c This driver should be linked with with bacolri.f and bacolri-aux.f
c (version 3) and a problem definition file, e.g., RCDsys.f
c
c-------------------------------------------------------------------
c Example Python plotting code using serialized B-spline data:
c
c   import matplotlib as mpl
c   mpl.use('AGG')  # for systems not running a GUI
c   import matplotlib.pyplot as plt
c   import numpy as np
c   from scipy.interpolate import splev
c
c   with open('Bsplines') as f: lines = f.readlines()
c   npde = int(lines[0])
c   xbs = np.fromstring(lines[1], sep=' ')
c   y   = np.fromstring(lines[2], sep=' ').reshape((npde,-1), order='F')
c   p = int(lines[3])
c
c   x = np.linspace(0, 1, 1000)
c   for i in range(npde):
c     u = splev(x, (xbs,y[i],p))
c     plt.plot(x, u, label='$u_{%d}$'%(i+1))
c
c   plt.xlabel('$x$')
c   plt.ylabel('$u$')
c   plt.title('$u(t,x)$ at $t=$tout')
c   plt.legend()
c   plt.grid()
c   plt.savefig('curve.png')
c
c scipy.interpolate.splev is an interface to SPLEV of FITPACK.
c-----------------------------------------------------------------------
        implicit none
c constants:
        integer                 npde
        parameter              (npde = 4)
c                               number of PDEs in the system.
        double precision        xa
        parameter              (xa = 0d0)
c                               the left spatial boundary point.
        double precision        xb
        parameter              (xb = 1d0)
c                               the right spatial boundary point.
        double precision        tstart
        parameter              (tstart = 0d0)
c                               the start of the time domain, used to
c                               initialise t0.
        integer                 nintinit
        parameter              (nintinit = 20)
c                               the number of subintervals in the
c                               initial uniform mesh.
        integer                 nintmx
        parameter              (nintmx = 200)
c                               maximum number of subintervals that
c                               bacolri may use
        integer                 nderiv
        parameter              (nderiv = 1)
c                               number of spatial partial derivatives
c                               to output approximations for, along
c                               with the solution value approximations.
c                               note that the solution approximation is
c                               a piecewise polynomial of degree kcol-1,
c                               so nderiv is limited.
        integer                 nxoutmx
        parameter              (nxoutmx = 1001)
c                               this is an arbitrary maximum number of
c                               points at which this driver will output
c                               solution or spatial derivative
c                               approximations.
        integer                 est
        parameter              (est = 1)
c                               set est = 1 to use the LOI/LE scheme,
c                               or set est = 2 to use the SCI/ST scheme.
        integer                 kcol
        parameter              (kcol = 4)
c                               kcol is the number of collocation
c                               points to be used in each subinterval.
c                                  3 <= kcol <= 10
c                               A value of 4 is most efficient in
c                               typical cases, and higher when
c                               a sharp tolerance is enforced.
c                               The sci scheme can have issues finding
c                               an initial mesh at kcol = 3,
c                               and the loi scheme can be extremely
c                               inefficient with sharp tolerances at
c                               kcol = 3.
        integer                 maxvec 
        parameter              (maxvec = npde*(kcol*nintmx+2)) 
c                               the maximum size of the vector of
c                               B-spline coefficients, y.
        integer                 lrp
        parameter              (lrp = 55+8*npde*npde*nintmx*kcol+39*npde
     +                              + 9*kcol+28*nintmx+24*nintmx*kcol
     +                              + 4*npde*npde*nintmx*kcol*kcol
     +                              + 17*npde*nintmx*kcol+21*npde*npde
     +                              + 6*nintmx*kcol*kcol
     +                              + 7*npde*nintmx
     -                              + 2+3*kcol+kcol*kcol-18*nintmx
     -                              - nintmx*kcol*kcol
     -                              - 9*nintmx*kcol-3*npde*nintmx)
c                               length of bacolri's floating point
c                               work array, rpar.
c                               the lines with a `-' continuation
c                               character may be uncommented if only
c                               the loi estimate is used.
        integer                 lip
        parameter              (lip = 102+3*npde*(nintmx*kcol+4))
c                               length of bacolri's integer
c                               work array, ipar.
        integer                 lcp
        parameter              (lcp = npde*(kcol*nintmx+2)
     +                            + npde*npde*(4+nintmx*kcol*(kcol+2)))
c                               lcp is the size of the cpar complex
c                               work array. 
        integer                 lenwrk
        parameter              (lenwrk = (kcol+2)*(nderiv+1)
     +                                 + kcol*(nintmx+1)+4)
c                               length of the work array used by VALUES
c                               after a successful return from bacolri.
        integer                 lmf
        parameter              (lmf = 12)
c                               lmf is the length of the mflag vector.
        integer                 ixbs
        parameter              (ixbs = 33)
c                               the breakpoint sequence used by BACOLRI's
c                               B-splines is stored in rpar(ipar(ixbs)).
c-----------------------------------------------------------------------
c bacolri input:
        double precision        t0
c                               t0 < tout is the initial time for a
c                               given call to bacolri.
        double precision        tout
c                               tout is the desired output time.
        double precision        atol(npde)
c                               atol is the absolute error tolerance
c                               request and is a scalar quantity if
c                               mflag(2) = 0.
        double precision        rtol(npde)
c                               rtol is the relative error tolerance
c                               request and is a scalar quantity if
c                               mflag(2) = 0.
        integer                 nint
c                               nint is the number of subintervals
c                               defined by the spatial mesh x.
        double precision        x(nintmx+1)
c                               x is the spatial mesh which divides the
c                               interval [xa,xb] as:
c                               xa = x(1) < x(2) < .. < x(nint+1) = xb
        integer                 mflag(lmf)
c                               this vector is used for more specific
c                               settings. see documentation in bacolri.f
        external                f
        external                bndxa
        external                bndxb
        external                uinit
        external                dum
c                               PDE system definition subroutines.
c                               see burg1.f for an example, or the
c                               bacolri.f preamble for descriptions.
c-----------------------------------------------------------------------
c bacolri work storage:
        double precision        rpar(lrp)
c                               rpar is a floating point work array
c                               of size lrp.
        integer                 ipar(lip)
c                               ipar is an integer work array
c                               of size lip.
        double complex          cpar(lcp)
c                               cpar is a complex work array of
c                               size lcp.
c-----------------------------------------------------------------------
c bacolri output:
        double precision        y(maxvec)
c                               on successful return from bacolri, y is
c                               the vector of time dependent B-spline
c                               coefficients at the current time t0.
        integer                 idid
c                               idid is the bacolri exit status flag
c-----------------------------------------------------------------------
c miscellany:
        double precision        uout(npde*nxoutmx*(nderiv+1))
c                               contains solution and spatial derivative
c                               approximations, which are obtained from
c                               the time dependent B-spline coefficients
c                               returned by bacolri (in y) through
c                               calling the subroutine values.
        double precision        valwrk(lenwrk)
c                               valwrk is the work array for values.
        double precision        xout(nxoutmx)
c                               xout is a set of spatial points for
c                               which output is desired.
        integer                 nxout
c                               nxout is the number of uniformly spaced
c                               points at which the user desired output
c                               at tout.
        integer                 i, j, k, ij
c                               loop indices
        character*3             eststr
c                               eststr is printed in the input summary.
c                               user input: follow the prompt.
        character               choice
c                               user input: follow the prompt.
        logical                 serialize
c                               should the B-spline information be saved
c                               to disk? the user decides.
c-----------------------------------------------------------------------
c problem specific variables:
        double precision        eps
        common /burger/         eps
c                               eps is used in the burg1.f, burg2.f and
c                               Cahn_Allen.f problem definition files.
c
        double precision        a1,a2,d1,d2,r,c,n,pe1,pe2
        common /rcdstuff/       a1,a2,d1,d2,r,c,n,pe1,pe2
c                               used in RCDsys.f.
c-----------------------------------------------------------------------
c subroutines called:
c                               bacolri
c                               values
c-----------------------------------------------------------------------
c  Common block initializations

      eps = 1d-3

      a1 = 30.d0
      a2 = 30.d0
      d1 = 1.5d0
      d2 = 1.2d0
      r = 1000.d0
      c = 0.96d0
      n = 1.d0
      pe1 = 10000.d0
      pe2 = 10000.d0

c-----------------------------------------------------------------------
c  Get some user input

      write(*,*) 'Please enter an output time'
      read(*,*,err=999) tout

      write(*,*) 'Please choose an error tolerance'
      read(*,*,err=999) atol(1)
      rtol(1) = atol(1)

      write(*,*) 'Output from this driver will be uniformly spaced in x'
      write(*,*) 'Please choose the number of output points'
      read(*,*,err=999) nxout
      if (nxout .gt. nxoutmx) then
         write(*,*) 'Warning: This driver limits nxout to', nxoutmx
         nxout = nxoutmx
      end if

      write(*,*) 'Would you like B-spline information saved to disk?',
     &           ' [y/N]'
      read(*,*,err=999) choice
      if (choice .eq. 'y' .or. choice .eq. 'Y') then
         serialize = .true.
      else
         serialize = .false.
      end if

c-----------------------------------------------------------------------
c  BACOLRI parameter initializations

      nint = nintinit
      t0 = tstart

c     Initialize rpar and ipar with zeros.
      do i = 1, lrp
         rpar(i) = 0d0
      end do
      do i = 1, lip
         ipar(i) = 0
      end do
      do i = 1, lcp
         cpar(i) = 0
      end do

c     Initialize the mflag vector.
      do i = 1, lmf
         mflag(i) = 0
      end do

c     Define the mesh based on a uniform subinterval size.
      x(1) = xa
      do i = 2, nint
         x(i) = xa + (i-1) * (xb-xa) / nint
      end do
      x(nint+1) = xb


c     Set mflag(3) = 1 for Dirichlet boundary conditions, 0 otherwise.
c     Recall that Dirichlet boundary conditions involve only function
c     values, i.e., U_x is not involved in bndxa or bndxb.
      mflag(3) = 0

      if (est .eq. 1) then
         mflag(5) = 0
         eststr = 'LOI'
      elseif (est .eq. 2) then
         mflag(5) = 1
         eststr = 'SCI'
      else
         stop
      endif

c-----------------------------------------------------------------------
c  Computation & output section

      write(*,'(/a)') 'THE INPUT IS'
      write(*,9000) 'kcol = ', kcol, ', nint0 = ', nint,
     &              ', npde = ', npde, ', tout = ', tout
      write(*,9001) 'atol = ', atol(1), ', rtol = ', rtol(1),
     &              ',                 ', eststr

      call bacolri(t0, tout, atol, rtol, npde, kcol, nintmx, nint,
     &            x, mflag, rpar, lrp, ipar, lip, cpar, lcp, y, idid,
     &            f, dum, bndxa, dum, bndxb, dum, uinit)

c     Check for errors running BACOLRI.
      if (idid .le. 0) goto 700
      write(*,9002) 'Number of subintervals in the current mesh:', nint

c     Output values column-wise (npde+1 columns).
      if (nxout .gt. 0) then

c        Define output spacial point sequence, uniformly spaced.
         xout(1) = xa
         do i = 2, nxout-1
            xout(i) = xa + (i-1) * (xb-xa) / (nxout-1)
         end do
         xout(nxout) = xb

c        Obtain approximate solution and spatial derivative values.
         call values(kcol, xout, nint, x, npde, nxout, nderiv, uout, y,
     &               valwrk)

          write(*,9003) 'THE OUTPUT IS'
          write(*,9004) 'XOUT', 'UOUT'
          do i = 1, nxout
             ij = (i-1)*npde
             write(*,*) xout(i), (uout(ij+k), k = 1, npde)
          end do
          if (nderiv .gt. 0) then
             do j = 1, nderiv
                write(*,9005) 'SPATIAL PARTIAL DERIVATIVE', j
                write(*,9004) 'XOUT', 'UOUT'
                do i = 1, nxout
                   ij = j*nxout*npde + (i-1)*npde
                   write(*,*) xout(i), (uout(ij+k), k = 1, npde)
                end do
             end do
          end if
      end if

c     Write B-spline information to file.
      if (serialize) then
         open(unit=20,file='Bsplines')
c        line 1: number of PDEs in the system
         write(20,*) npde
c        line 2: breakpoint/knot sequence
         write(20,*) (rpar(ipar(ixbs)+i-1), i = 1, (nint+1)*kcol+4)
c        line 3: coefficient values at time t0
         write(20,*) (y(i), i = 1, (nint*kcol+2)*npde)
c        line 4: degree of the B-spline interpolant
         write(20,*) kcol+1
         close(20)
      end if

      stop
c-----------------------------------------------------------------------
c  Unsuccessful return section
  700 continue
      write(*,'(a)')        'CANNOT PROCEED DUE TO ERROR FROM BACOLRI.'
      write(*,'(a,i21)')    'IDID         =', idid
      write(*,'(a,g25.18)') 'Failure at T =', t0
      write(*,'(a,i21)')    'NINT         =', nint
      write(*,'(a,i21)')    'Init. remesh.:', ipar(9)
      write(*,'(a,i21)')    'Remeshings   :', ipar(11)
      write(*,'(a,i21)')    'Cold starts  :', ipar(14)

      stop
c-----------------------------------------------------------------------
c  Invalid input error section
  999 continue
      write(*,*) 'Error: Improperly formatted input'

      stop
c-----------------------------------------------------------------------
c  Formats!
 9000 format( a, i2, a, i4, a, i3, a, e7.1 )
 9001 format( 2(a, e7.1), 2a )
 9002 format( a, i8 )
 9003 format( /a )
 9004 format( a13, a27 )
 9005 format( /a, i2 )
c-----------------------------------------------------------------------
      end

      subroutine dum
c     This is a dummy subroutine to be provided to bacolri in place of
c     derivf, difbxa and difbxb when using approximated derivatives.
      end
