c-----------------------------------------------------------------------
c driver-L2.f
c With the true solution defined in TRUU, this driver uses DQAGE from
c quadpack to approximate the L2-norm of the error in the approximate
c solution after BACOLI has returned.
c Some bits need to be modified to suit the current test being done,
c such as the parameters npde, nintmx, xa and xb,
c or common block declarations and initializations for problem-specific
c parameters.
c Also the output given here probably isn't what you want, so do
c something about that.
c
c When prompted to input the tolerance or epsilon, it is recommended to
c use engineering notation (like `1e-7'.)
c
c This driver should be compiled with:
c   some version of BACOLI, eg, bacoli.f
c   a problem definition file, eg, burg1.f
c   quadpack routines
c
c-----------------------------------------------------------------------
c Use preprocessor directives for certain aspects of this driver.
#ifndef NPDE
#  define NPDE 1
#endif
#ifndef NINTMX
#  define NINTMX 1000
#endif
#ifndef XA
#   define XA 0d0
#endif
#ifndef XB
#   define XB 1d0
#endif
#ifndef TOUT
#   define TOUT 1d0
#endif
c Define BACOL to use this driver with the old BACOL code, otherwise
c BACOLI will be used.
c Define DIRICHLET if the boundary conditions are such.
c burg1.f and burg2.f have Dirichlet boundary conditions, for example.
c-----------------------------------------------------------------------
        implicit none
c constants:
        integer                 kcolmx
c                               kcol is the number of collocation points
c                               to be used in each subinterval, which is
c                               equal to the degree of the piecewise
c                               polynomials minus one.
c                                  2 < kcol < 11
        parameter              (kcolmx = 11)
        integer                 npde
c                               number of pdes
        parameter              (npde = NPDE)
        integer                 nintmx
c                               maximum number of intervals required
        parameter              (nintmx = NINTMX)
        integer                 maxvec
c                               the dimension of the vector of
c                               B-spline coefficients
        parameter              (maxvec = npde*(nintmx*kcolmx+2))
        integer                 lrp
c                               uncomment the lines with a `-'
c                               continuation character
c                               and fix the parentheses if you
c                               only want to use the loi estimate.
#ifdef BACOL
        parameter              (lrp=134+nintmx*(35+35*kcolmx+31*npde
     +                               +38*npde*kcolmx+8*kcolmx*kcolmx)
     +                               +14*kcolmx+79*npde+npde*npde*(21
     +                               +4*nintmx*kcolmx*kcolmx+12*nintmx
     +                               *kcolmx+6*nintmx))
#else
        parameter              (lrp=113+57*npde+27*nintmx
     +                               +13*npde*npde+9*kcolmx+24*kcolmx
     +                               *nintmx+6*nintmx*kcolmx*kcolmx
     +                               +7*nintmx*npde+27*npde*nintmx
     +                               *kcolmx+2*npde*npde*nintmx*kcolmx
     +                               *kcolmx+4*npde*npde*kcolmx*nintmx)
c    -                               -15*nintmx+3*kcolmx+kcolmx*kcolmx
c    -                               -8*kcolmx*nintmx-3*nintmx*npde
c    -                               -nintmx*kcolmx*kcolmx+2
#endif
        integer                 lip
c                               see the comment for ipar

#ifdef BACOL
        parameter              (lip=115+npde*(nintmx*(2*kcolmx+1)+4))
#else
        parameter              (lip=101+npde*(nintmx*kcolmx)+2)

        external                f
        external                uinit
        external                derivf
        external                bndxa
        external                bndxb
        external                difbxa
        external                difbxb
#endif

        integer                 lenwrk
c                               the dimension of array work when we
c                               call values
        parameter              (lenwrk =(kcolmx+2)+kcolmx*(nintmx+1)+4)
        double precision        xa
c                               the left boundary point
        parameter              (xa = XA)
        double precision        xb
c                               the right boundary point
        parameter              (xb = XB)
c-----------------------------------------------------------------------
c bacoli input:
        integer                 kcol
c                               the number of collocation pointer per
c                               subinterval
        integer                 refkcol
        double precision        t0
c                               t0 < tout is the initial time.
        double precision        tout
c                               tout is the desired final 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                 nintinit
        integer                 nint
        integer                 refnint
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 [x_a,x_b] as: x_a = x(1) <
c                               x(2) < x(3) < ... < x(nint+1) = x_b.
        double precision        refx(nintmx+1)
        integer                 mflag(12)
c                               this vector of user input determines
c                               the interaction of bacoli with dassl.
c                               PDE system definition subroutines.
c-----------------------------------------------------------------------
c bacoli 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.
c-----------------------------------------------------------------------
c bacoli output:
        double precision        y(maxvec)
c                               on successful return from bacoli, y is
c                               the vector of bspline
c                               coefficients at the current time t0.
        double precision        refy(maxvec)
        integer                 idid
c                               idid is the bacoli exit status flag
c                               which is based on the exit status from
c                               dassl on error checking performed by
c                               bacoli on initialization.
c-----------------------------------------------------------------------
        double precision        yy(maxvec)
c                               put this in the common block, rather
c                               than y directly. the BAC/LE scheme is
c                               implemented by simply copying the second
c                               set of B-spline coefficients directly
c                               from rpar.
#ifdef BACOL
        integer                 iy
c                               used to index the B-spline coefficients
c                               inside of rpar.
c                               rpar(ipar(24)) has the lower order y
c                               rpar(ipar(43)) has the higher order y
#endif

        double precision        valwrk(lenwrk)
c                               valwrk is a work array in values
        double precision        tstart
c                               stores initial variables that the
c                               subroutine bacoli will change, so they
c                               can be used in second call to bacoli
        integer                 est
c                               1 to use LOI, 2 to use SCI
        integer                 i, j, k, m
c                               loop counters
        character*5             eststr
        integer                 khigh
        integer                 comp
c-----------------------------------------------------------------------
c L2-norm error calculation:
        integer                 neval, ier, last
        double precision        result, maxresult, abserr
        common /L21/       x, yy, valwrk, t0
        common /L22/       atol, rtol, nint, kcol, khigh
        common /L23/       refx, refy, comp, refkcol, refnint

c-----------------------------------------------------------------------
c problem specific variables:
        double precision        eps
        common /burger/         eps
c
        double precision        a1,a2,d1,d2,r,c,n,pe1,pe2
        common /rcdstuff/       a1,a2,d1,d2,r,c,n,pe1,pe2
c-----------------------------------------------------------------------
c timing variables:
        integer                 nruns
        real                    dtvals(2), time

c-----------------------------------------------------------------------
c subroutines called:
c                               bacol
c                               bacoli
c                               values
c                               scaledL2
c                               dtime
c
c-----------------------------------------------------------------------
c Fortran functions used:
c                               dble
c                               dabs
c                               max
c
c-----------------------------------------------------------------------
c  Initializations
      tstart = 0d0
      nintinit = 100
      tout = TOUT

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

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

      write(*,9000) 'Pick a kcol'
      read(*,*) kcol

      write(*,9000) 'Pick the error control mode'
      write(*,8999) '1. LOI/LE'
      write(*,8999) '2. SCI/ST'

      read(*,*) est

      write(*,9000) 'Please input the tolerance to use'
      read(*,*) atol(1)
      rtol(1) = atol(1)

      write(*,9000) 'Please input epsilon, too'
      read(*,*) eps

c-----------------------------------------------------------------------
c  Even more initializations

c  vectors.
      do i = 1, 8
         mflag(i) = 0
      end do
c  set mflag(5) = 1 for dirichlet boundary conditions, 0 otherwise
#ifdef DIRICHLET
      mflag(5) = 1
#endif

c     Set error estimation scheme for BACOLI. For BACOL, set different
c     kcol value for use in call to values in order to extract higher
c     or lower-order interpolant.
#ifndef BACOL
      khigh = kcol
      if (est .eq. 1) then
         mflag(8) = 0
         eststr = 'LOI'
      elseif (est .eq. 2) then
         mflag(8) = 1
         eststr = 'SCI'
      endif
      mflag(9) = 3
#else
      if (est .eq. 1) then
         mflag(8) = 0
         eststr = 'BAC'
         khigh = kcol
      elseif (est .eq. 2) then
         mflag(8) = 1
         eststr = 'BLE'
         khigh = kcol + 1
      endif
#endif

c       This is to restrict the order used so that a solution is 
c        actually reached.      
#ifdef BACOLISCHROD
      mflag(7) = 1
      ipar(15) = 2
#endif

c-----------------------------------------------------------------------
c  Computation & output section
      write(*,'(/a)') 'THE INPUT IS'

      write(*,9001) 'kcol = ', kcol, ', nint0 = ', nintinit,
     &              ', npde = ', npde, ', tout = ', tout
      write(*,9002) 'atol = rtol = ', rtol(1),
     &              ',    ', eststr, ',    epsilon = ', eps

c     choose nruns appropriate for expected runtime
      if      (atol(1) .gt. 1d-2) then
         nruns = 80
      else if (atol(1) .gt. 1d-4) then
         nruns = 40
      else if (atol(1) .gt. 1d-6) then
         nruns = 20
      else if (atol(1) .gt. 1d-8) then
         nruns = 10
      else
         nruns = 5
      end if

c     CALCULATION
      call dtime(dtvals, time)
      do j = 1, nruns
c        some per-run initialisations
         t0 = tstart
c         rpar(1) = tout part of a no longer existing flag
         nint = nintinit
         x(1) = xa
         do i = 2, nint
            x(i) = xa + ((i-1) * (xb - xa)) / nint
         end do
         x(nint+1) = xb
         mflag(1) = 0


#ifdef BACOL
      call bacol(t0, tout, atol, rtol, npde, kcol, nintmx, nint,
     &            x, mflag, rpar, lrp, ipar, lip, y, idid)
#else
      call bacoli(t0, tout, atol, rtol, npde, kcol, nintmx, nint,
     &            x, mflag, rpar, lrp, ipar, lip, y, idid, f, derivf,
     &            bndxa, difbxa, bndxb, difbxb, uinit)
#endif
c        Check for errors running BACOLI.
         if (idid .le. 0) goto 700
      end do
      call dtime(dtvals, time)
      time = time / nruns

#ifdef BACOL
c     Take either the higher or lower order collocation solution out
c     of the work array to effect either the BAC/ST or BAC/LE scheme.
      if (est .eq. 1) iy = 24
      if (est .eq. 2) iy = 43
      call dcopy(maxvec, rpar(ipar(iy)), 1, yy, 1)
#else
      call dcopy(maxvec, y, 1, yy, 1)
#endif

c-----------------------------------------------------------------------
c     ACCESS REFERENCE SOLUTION
c     Reference solution is computed by BACOL with atol = rtol = 1d-12,
c     kcol = 3
      open(102, file='parameters.txt', status='UNKNOWN')
      read(102, *) refkcol
      read(102, *) refnint
c      write(0,*) '----------------------'
c      write(0,*) 'Checking kcol, nint'
c      write(0,*) refkcol, ', ', refnint
      close(102)

      open(101, file='reference_x.txt', status='UNKNOWN')
      do i = 1, refnint+1
        read(101, *) refx(i)
      end do
c      write(0,*) '----------------------'
c      write(0,*) 'Checking x'

c      do i = 1, refnint+1
c        write(0,*) refx(i)
c      end do
      close(101)

      open(100, file='reference_coefs.txt', status='UNKNOWN')
      do i = 1, npde*(refnint*refkcol+2)
        read(100, *) refy(i)
      end do
c      write(0,*) '----------------------'
c      write(0,*) 'Checking x'
c      do i = 1, npde*(nintmx*refkcol+2)
c        write(0,*) refy(i)
c      end do
c      write(0,*) '----------------------'
      close(100)



c-----------------------------------------------------------------------
c     Check the L2-norm of the error at tout
c     (to see how much success bacol had in controlling it.)
c   See how this works without this section, otherwise look into it

c     Scaled L2-norm.
      do comp = 1,1
        call scaledL2(xa, xb, atol(1), rtol(1),
     &              result, abserr, neval, ier, last)

        if (comp .eq. 1) then
          maxresult = result
        else if (result .gt. maxresult) then
          maxresult = result
        end if

c     Check for errors in approximating the scaled L2-norm.
        if (ier .eq. 6) then
           write(*,*) 'Input error for scaled L2 norm approximation.'
        elseif ((ier .ne. 0 .and. ier .ne. 2)
     &   .or.(abserr*1d2 .gt. result)) then
           write(*,*) 'Warning: scaled L2 norm likely unreliable.'
           write(*,*) 'ier    = ', ier
           write(*,*) 'abserr = ', abserr
        endif
      end do

c-----------------------------------------------------------------------
c  Sucessful return section

      write(*,9006) 'IDID         = ', idid
      write(*,9007) 'DTIME        = ', time
      write(*,9004) 'Final time   = ', t0
      write(*,9005) 'Final nint   = ', nint
      write(*,9005) 'ipar(irshin) = ', ipar(9)
      write(*,9005) 'ipar(irmesh) = ', ipar(11)
      write(*,9009) 'L2 error     = ', result
      write(*,9005) 'quad nevals  = ', neval
      write(*,9005) 'quad nints   = ', last
      write(0,*) 'SUCCESS', maxresult, time, atol(1)

      stop

c-----------------------------------------------------------------------
c  Unsucessful return section
  700 continue
#ifndef BACOL
      if (est .eq. 1) then
          write(*,'(a)') 'CANNOT PROCEED DUE TO ERROR FROM BACOLI/LE.'
      else
          write(*,'(a)') 'CANNOT PROCEED DUE TO ERROR FROM BACOLI/ST.'
      end if
#else
      if (est .eq. 1) then
          write(*,'(a)') 'CANNOT PROCEED DUE TO ERROR FROM BACOL/ST.'
      else
          write(*,'(a)') 'CANNOT PROCEED DUE TO ERROR FROM BACOL/LE.'
      end if
#endif
      write(*,'(a,g25.18)') 'Failure at T = ', t0
      write(*,'(a,i3)')     'IDID = ', idid
      write(*,'(a,i10)')    'NINT = ', nint

      write(0,*) 'FAILURE', idid, kcol, atol(1), ipar(11)+ipar(9)

      stop
c-----------------------------------------------------------------------
c  Formats!
 8999 format( a )
 9000 format( /a )
 9001 format( a, i2, a, i4, a, i3, a, es7.1 )
 9002 format( a, es7.1, 3a, es7.1 )
 9004 format( a, es20.8 )
 9005 format( a, i10 )
 9006 format( a, i4 )
 9007 format( a, f15.3 )
 9009 format( a, es25.18 )
c-----------------------------------------------------------------------
      end
      subroutine scaledL2(a, b, atol, rtol, result, abserr, neval,
     &                    ier, last)
cliffsnotes
c     This subroutine approximates the L2 norm of |GE| / (1 + |U|)
c     after a call to bacol using dqage from quadpack.
c     On sucessful return, result should contain the norm.
c     If ier is (0, 0), it should be reliable.
c     See dqage documentation for more details.
      implicit none
      double precision  a, b, atol, rtol
      double precision  result, abserr
      integer           neval, ier, last
      double precision  scaledErr
      external          scaledErr
constants
      integer           key
      parameter        (key = 6)
      integer           limit
      parameter        (limit = 500)
crucial variables
      integer           iord(limit)
      double precision  alist(limit), blist(limit), elist(limit)
      double precision  rlist(limit)
code
c     approximate integral of (|GE| / (1 + |U|)) ** 2
      call dqage(scaledErr, a, b, 0d0, max(1d-13, rtol*1d-8),
     &           key, limit, result, abserr, neval, ier,
     &           alist, blist, rlist, elist, iord, last)
      result = sqrt(result)
      return
      end
c-----------------------------------------------------------------------
      double precision function scaledErr(x)
      implicit none
      double precision  x
constants
      integer           kcolmx
      parameter        (kcolmx = 11)
      integer           npde
      parameter        (npde = NPDE)
      integer           nintmx
      parameter        (nintmx = NINTMX)
      integer           maxvec
      parameter        (maxvec = npde*(nintmx*kcolmx+2))
      integer           lenwrk
      parameter        (lenwrk = (kcolmx+2)+kcolmx*(nintmx+1)+4)
      external          truu
commons
      integer           i
      double precision  mesh(nintmx+1), y(maxvec), valwrk(lenwrk)
      double precision  tout
      double precision  atol(npde), rtol(npde)
      integer           nint, kcol, khigh, comp
      double precision  refmesh(nintmx+1), refy(maxvec)
      integer           refnint, refkcol
      common /L21/       mesh, y, valwrk, tout
      common /L22/       atol, rtol, nint, kcol, khigh
      common /L23/       refmesh, refy, comp, refkcol, refnint
crucial variables
      double precision  U(npde,1,1)
      double precision  refU(npde,1,1)
code
      call values(khigh, x, nint, mesh, npde, 1, 0, U, y, valwrk)
      call values(refkcol, x, refnint, refmesh, npde, 1, 0, refU, refy,
     &            valwrk)

      scaledErr = (U(comp,1,1) - refU(comp,1,1)) /
     &                (1d0 + abs(U(comp,1,1)))
      scaledErr = scaledErr * scaledErr
      end