program adapt
!
!*******************************************************************************
!
!! ADAPT solves a 1D problem using an adaptive finite element method.
!
!
!  Discussion:
!
!    The equation to be treated is:
!
!    - d/dX (P dU/dX) + Q U  =  F
!
!    by the finite-element method using piecewise linear basis
!    functions.
!
!    An adaptive method is used to try to reduce the maximum
!    error by refining the mesh in certain places.
!
!    Here U is an unknown scalar function of X defined on the
!    interval [XL,XR], and P, Q and F are given functions of X.
!
!    The values of U at XL and XR are also specified.
!
!
!  The interval [XL,XR] is "meshed" with N+1 points,
!
!  XN(0) = XL, XN(1)=XL+H, XN(2)=XL+2*H, ..., XN(N)=XR.
!
!  This creates N subintervals, with interval number 1
!  having endpoints XN(0) and XN(1), and so on up to interval
!  N, which has endpoints XN(N-1) and XN(N).
!
!
!  The algorithm employed by ADAPT tries to guarantee a certain amount
!  of accuracy by examining the current solution, estimating the error
!  in each subinterval, and, if necessary, subdividing one or more
!  subintervals and repeating the calculation.
!
!  We can think of the adaptive part of the algorithm as a refined
!  problem.  The program re-solves the problem on the pair of
!  intervals J and J+1, which extend from node J-1 to node J+1.
!  The values of U that were just computed at nodes J-1 and J+1
!  will be used as the boundary values for this refined problem.
!  The intervals J and J+1 will each be evenly divided into NY
!  smaller subintervals.  This boundary value problem is solved,
!  and the derivatives of the original and refined solutions are
!  then compared to get an estimate of the error.
!
!  Sample problems:
!
!    #1: u = x, p=1, q=0, f=0, ibc=3, ul=0, ur=1.
!    The program should find the solution exactly, and the
!    adaptive code should find that there is no reason to
!    subdivide any interval.
!
!    #2: u = x*x, p=1, q=0, f=-2, ibc=3, ul=0, ur=1.
!    This problem should find the solution exactly, and
!    the adaptive code should again find there is nothing
!    to do.
!
!    #3: u = sin(pi*x/2), p=1, q=0, f=0.25*pi*pi*sin(pi*x/2), ul=0, ur=1.
!
!    #4: u = cos(pi*x/2), p=1, q=0, f=0.25*pi*pi*cos(pi*x/2), ul=1, ur=0.
!
!    #5: u = x**(beta+2)/((beta+2)*(beta+1)), p=1, q=1,
!    f = -x**beta + (x**(beta+2))/((beta+2)*(beta+1)),
!    ul = 0, ur=1/((beta+2)*(beta+1))
!
!  (beta must be greater than -2, and not equal to -1)
!
!    #6: u = atan((x-0.5)/eps), p=1, q=0,
!    f =  2*eps*(x-0.5) / (eps**2 + (x-0.5)**2) **2,
!    ul = u(0), ur=u(1)
!
!  Parameters:
!
!    REAL ADIAG(NU).
!    ADIAG(I) is the "diagonal" coefficient of the I-th
!    equation in the linear system.  That is, ADIAG(I) is
!    the coefficient of the I-th unknown in the I-th equation.
!
!    REAL ALEFT(NU).
!    ALEFT(I) is the "left hand" coefficient of the I-th
!    equation in the linear system.  That is, ALEFT(I) is the
!    coefficient of the (I-1)-th unknown in the I-th equation.
!    There is no value in ALEFT(1), since the first equation
!    does not refer to a "0-th" unknown.
!
!    REAL ARITE(NU).
!    ARITE(I) is the "right hand" coefficient of the I-th
!    equation in the linear system.  ARITE(I) is the coefficient
!    of the (I+1)-th unknown in the I-th equation.  There is
!    no value in ARITE(NU) because the NU-th equation does not
!    refer to an "NU+1"-th unknown.
!
!    REAL ETA(N).
!    ETA(I) is the error estimate for interval I.  It is computed
!    as the sum of two quantities, one associated with the left
!    and one with the right node of the interval.
!
!    REAL F(N+1) or F(NU).
!    ASSEMB stores into F the right hand side of the linear
!    equations.
!    SOLVE replaces those values of F by the solution of the
!    linear equations.
!
!    REAL FY(M).
!    FY is the right hand side of the linear system of the refined
!    problem.
!
!    REAL H(N)
!    H(I) is the length of subinterval I.  This code uses
!    equal spacing for all the subintervals.
!
!    REAL HY(M).
!    HY(I) is the length of subinterval I in the refined problem.
!
!    INTEGER IBC.
!    IBC declares what the boundary conditions are.
!    1, at the left endpoint, U has the value UL,
!       at the right endpoint, U' has the value UR.
!    2, at the left endpoint, U' has the value UL,
!       at the right endpoint, U has the value UR.
!    3, at the left endpoint, U has the value UL,
!       and at the right endpoint, U has the value UR.
!    4, at the left endpoint, U' has the value UL,
!       at the right endpoint U' has the value UR.
!
!    INTEGER IBCY.
!    IBCY declares the boundary conditions for the refined problem
!    which should always be that the value of U is specified at
!    both the left and right endpoints.  This corresponds to a
!    value of IBCY = 3.
!
!    INTEGER INDX(0:N).
!    For a node I, INDX(I) is the index of the unknown
!    associated with node I.
!    If INDX(I) is equal to -1, then no unknown is associated
!    with the node, because a boundary condition fixing the
!    value of U has been applied at the node instead.
!    Unknowns are numbered beginning with 1.
!    If IBC is 2 or 4, then there is an unknown value of U
!    at node 0, which will be unknown number 1.  Otherwise,
!    unknown number 1 will be associated with node 1.
!    If IBC is 1 or 4, then there is an unknown value of U
!    at node N, which will be unknown N or N+1,
!    depending on whether there was an unknown at node 0.
!
!    INTEGER INDY(0:M).
!    INDY(I) records the index of the unknown associated with
!    node I for the refined problem.
!
!    INTEGER JADD(N).
!    JADD(I) is 1 if the error estimates show that interval I
!    should be subdivided.
!
!    INTEGER M.
!    M is the number of subintervals used in the refined problem.
!    M is equal to NY for computations centered at node 0 or node N,
!    and otherwise, M is equal to 2*NY.
!
!    Integer N
!    The number of subintervals into which the interval
!    [XL,XR] is broken.
!
!    INTEGER NL.
!    The number of basis functions used in a single
!    subinterval.  (NL-1) is the degree of the polynomials
!    used.  For this code, NL is fixed at 2, meaning that
!    piecewise linear functions are used as the basis.
!
!    INTEGER NODE(N,NL).
!    For each subinterval I:
!    NODE(I,1) is the number of the left node, and
!    NODE(I,2) is the number of the right node.
!
!    INTEGER NODE(M,NL).
!    NODEY performs the same function for the refined problem that
!    NODE performs for the full problem, recording the node numbers
!    associated with a particular subinterval.
!
!    INTEGER NQUAD
!    The number of quadrature points used in a subinterval.
!
!    INTEGER NU.
!    NU is the number of unknowns in the linear system.
!    Depending on the value of IBC, there will be N-1,
!    N, or N+1 unknown values, which are the coefficients
!    of basis functions.
!
!    INTEGER NUY.
!    The number of unknowns in the refined problem.
!
!    INTEGER NY.
!    NY is the number of subintervals into which a given interval
!    will be subdivided, before solving the refined probelm.
!
!    REAL TOL.
!    A tolerance that is used to determine whether the estimated
!    error in an interval is so large that it should be subdivided
!    and the problem solved again.
!
!    REAL UL.
!    If IBC is 1 or 3, UL is the value that U is required
!    to have at X = XL.
!    If IBC is 2 or 4, UL is the value that U' is required
!    to have at X = XL.
!
!    REAL UR.
!    If IBC is 2 or 3, UR is the value that U is required
!    to have at X = XR.
!    If IBC is 1 or 4, UR is the value that U' is required
!    to have at X = XR.
!
!    REAL WQUAD(NQUAD).
!    WQUAD(I) is the weight associated with the I-th point
!    of an NQUAD point Gaussian quadrature rule.
!
!    REAL XL.
!    XL is the left endpoint of the interval over which the
!    differential equation is being solved.
!
!    REAL XN(0:N).
!    XN(I) is the location of the I-th node.  XN(0) is XL,
!    and XN(N) is XR.
!
!    REAL XQUAD(NMAX,NQUAD)
!    XQ(I,J) is the location of the J-th quadrature point
!    in interval I.
!
!    REAL XQUADY(NMAY,NQUAD).
!    XQUADY(I,J) is the location of the J-th quadrature point
!    in subinterval I of the refined problem.
!
!    REAL XR.
!    XR is the right endpoint of the interval over which the
!    differential equation is being solved.
!
!    REAL YN(0:M).
!    YN(I) is the location of the I-th node in the refined
!    problem.
!
  implicit none
!
  integer, parameter :: nl = 2
  integer, parameter :: nmax = 30
  integer, parameter :: nquad = 2
!
  real adiag(nmax)
  real aleft(nmax)
  real arite(nmax)
  real eta(nmax)
  real f(nmax)
  real h(nmax)
  integer ibc
  integer indx(0:nmax)
  integer jadd(nmax)
  integer kount
  integer n
  integer node(nmax,nl)
  integer nu
  real tol
  real ul
  real ur
  real wquad(nquad)
  real xl
  real xn(0:nmax)
  real xquad(nmax,nquad)
  real xr
  real xt(0:nmax)
!
  call timestamp ( )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'ADAPT'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'Solve the two-point boundary value problem:'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  -d/dx P du/dx + Q U  =  F'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'on the interval [0,1], specifying the value'
  write ( *, '(a)' ) 'of U at each endpoint.'
  write ( *, '(a)' ) ' '
  write ( *, '(a,i6)' ) 'The number of basis functions per element is ', nl
  write ( *, '(a)' ) ' '
  write ( *, '(a,i6)' ) 'The number of quadrature points per element is ', nquad
!
!  Start out with just 4 subintervals.
!
  n = 4
!
!  Initialize values that define the problem.
!
  call init ( ibc, n, tol, ul, ur, xl, xn, xr )
!
!  Start the iteration counter off at 0.
!
  kount = 0
!
!  Begin the next iteration.
!
  do

    kount = kount + 1

    write ( *, '(a)' ) ' '
    write ( *, '(a,i6,a)' ) 'Begin new iteration with ',n,' nodes.'
    write ( *, '(a)' ) ' '
!
!  Solve the regular problem.
!
    call solvex ( adiag, aleft, arite, f, h, ibc, indx, kount, n, nl, &
      nmax, node, nquad, nu, ul, ur, wquad, xn, xquad )
!
!  Solve N subproblems to get the error estimators.
!
    call solvey ( eta, f, h, n, nu, ul, ur, xn )
!
!  Examine the error estimators, and see how many intervals should
!  be subdivided.
!
    call subdiv ( eta, jadd, kount, n, nmax, tol, xn, xt )
!
!  Solve the problem again, with the new nodes.
!
  end do

  stop
end
subroutine assemb ( adiag, aleft, arite, f, h, n, indx, node, nu, nl, &
  nquad, nmax, ul, ur, wquad, xn, xquad )
!
!*******************************************************************************
!
!! ASSEMB assembles the global matrix.
!
  implicit none
!
  integer n
  integer nl
  integer nmax
  integer nquad
  integer nu
!
  real adiag(nu)
  real aleft(nu)
  real arite(nu)
  real aij
  real f(nu)
  real ff
  real h(n)
  real he
  integer i
  integer ie
  integer ig
  integer il
  integer indx(0:n)
  integer iq
  integer iu
  integer jg
  integer jl
  integer ju
  integer node(nmax,nl)
  real phii
  real phiix
  real phij
  real phijx
  real pp
  real qq
  real ul
  real ur
  real wquad(nquad)
  real wquade
  real xleft
  real xn(0:n)
  real xquad(nmax,nquad)
  real xquade
  real xrite
!
  external ff
  external pp
  external qq
!
!  Zero out the entries.
!
  f(1:nu) = 0.0
  aleft(1:nu) = 0.0
  arite(1:nu) = 0.0
  adiag(1:nu) = 0.0
!
!  For each interval,
!
  do ie = 1,n

    he = h(ie)
    xleft = xn(node(ie,1))
    xrite = xn(node(ie,2))
!
!  For each quadrature point in the interval,
!
    do iq = 1,nquad

      xquade = xquad(ie,iq)
      wquade = wquad(iq)
!
!  Pick a basis function which defines the equation,
!
      do il = 1,nl

        ig = node(ie,il)
        iu = indx(ig)

        if ( iu > 0 ) then

          call phi ( il, xquade, phii, phiix, xleft, xrite )
          f(iu) = f(iu)+he*wquade*ff(xquade)*phii
!
!  Take care of boundary conditions specifying the value of U'.
!
          if ( ig == 0 ) then
            f(iu) = f(iu)-pp(0.0)*ul
          else if ( ig == n ) then
            f(iu) = f(iu)+pp(1.0)*ur
          end if
!
!  Pick a basis function which defines the coefficient
!  being computed.
!
          do jl = 1, nl

            jg = node(ie,jl)
            ju = indx(jg)
            call phi ( jl, xquade, phij, phijx, xleft, xrite )

            aij = he * wquade * &
               ( pp(xquade) * phiix * phijx &
               + qq(xquade) * phii * phij )
!
!  Decide where the coefficient is to be added.
!
            if ( ju <= 0 ) then
 
              if ( jg == 0 ) then
                f(iu) = f(iu) - aij * ul
              else if ( jg == n ) then
                f(iu) = f(iu) - aij * ur
              end if

            else if ( iu == ju ) then
              adiag(iu) = adiag(iu) + aij
            else if ( iu > ju ) then
              aleft(iu) = aleft(iu) + aij
            else
              arite(iu) = arite(iu) + aij
            end if

          end do

        end if

      end do

    end do

  end do

  return
end
function ff(x)
!
!*******************************************************************************
!
!! FF evaluates the function F in the differential equation.
!
!    - d/dx (p du/dx) + q u  =  f
!
!  at the point X.
!
  implicit none
!
  real beta
  real eps
  real ff
  integer iprob
  real, parameter :: pi = 3.14159265E+00
  real x
!
!  Find out which problem we're working on.
!
  call getprb(iprob)
!
  if ( iprob == 1 ) then

    ff = 0.0E+00

  else if ( iprob == 2 ) then

    ff = -2.0E+00 * x

  else if ( iprob == 3 ) then

    ff = 0.25E+00 * pi**2 * sin ( 0.5 * pi * x )

  else if ( iprob == 4 ) then

    ff = 0.25 * pi**2 * cos ( 0.5 * pi * x )

  else if ( iprob == 5 ) then

    call getbet ( beta )

    ff = - ( x**beta ) + ( x**(beta+2.0) ) &
      / ( ( beta + 2.0 ) * ( beta + 1.0 ) )

  else if ( iprob == 6 ) then

    call geteps ( eps )
    ff = 2.0 * eps * ( x - 0.5 ) / ( eps**2 + ( x - 0.5 )**2 )**2

  end if

  return
end
subroutine geomet ( h, ibc, indx, n, nl, nmax, node, nquad, nu, wquad, xn, &
  xquad )
!
!*******************************************************************************
!
!! GEOMET sets up some of the geometric information for the problem.  
!
!
!  Note, however, that the location of the nodes
!  is done outside of this routine, and, in fact, before this
!  routine is called.
!
  implicit none
!
  integer n
  integer nl
  integer nmax
  integer nquad
!
  real alfa
  real h(n)
  integer i
  integer ibc
  integer igl
  integer igr
  integer indx(0:n)
  integer node(nmax,nl)
  integer nu
  real wquad(nquad)
  real xl
  real xn(0:n)
  real xquad(nmax,nquad)
  real xr
!
!  Store in NODE the fact that interval I has node I-1
!  as its left endpoint, and node I as its right endpoint.
!
  do i = 1,n
    node(i,1) = i-1
    node(i,2) = i
  end do
!
!  For every node that is associated with an unknown, we
!  record the number of the unknown in INDX.
!
  nu = 0
  do i = 0, n

    if ( i == 0 .and. ( ibc == 1 .or. ibc == 3 ) ) then
      indx(i) = -1
    else if ( i == n .and. ( ibc == 2 .or. ibc == 3 ) ) then
      indx(i) = -1
    else
      nu = nu+1
      indx(i) = nu
    end if

  end do
!
!  We compute the width of each interval.
!
  do i = 1, n
    igl = node(i,1)
    igr = node(i,2)
    h(i) = xn(igr) - xn(igl)
  end do
!
!  We compute the location of the quadrature points in each
!  interval.
!
  do i = 1, n

    xl = xn(node(i,1))
    xr = xn(node(i,2))

    if ( nquad == 1 ) then
      xquad(i,1) = 0.5 * ( xl + xr )
    else if ( nquad == 2 ) then
      alfa = -0.577350
      xquad(i,1) = ( (1.0-alfa)*xl + (1.0+alfa)*xr ) / 2.0
      alfa = +0.577350
      xquad(i,2) = ( (1.0-alfa)*xl + (1.0+alfa)*xr ) / 2.0
    else if ( nquad == 3 ) then
      alfa = -0.774597
      xquad(i,1) = ( (1.0-alfa)*xl + (1.0+alfa)*xr ) / 2.0
      xquad(i,2) = 0.5*(xl+xr)
      alfa = +0.774597
      xquad(i,3) = ( (1.0-alfa)*xl + (1.0+alfa)*xr ) / 2.0
    end if

  end do
!
!  Store the weights for the quadrature rule.
!
  if ( nquad == 1 ) then
    wquad(1) = 1.0
  else if ( nquad == 2 ) then
    wquad(1) = 0.5
    wquad(2) = 0.5
  else if ( nquad == 3 ) then
    wquad(1) = 4.0 / 9.0
    wquad(2) = 5.0 / 18.0
    wquad(3) = 4.0 / 9.0
  end if

  return
end
subroutine getbet ( beta )
!
!*******************************************************************************
!
!! GETBET returns the value of BETA, for use by problem 5.
!
  implicit none
!
  real beta
!
  beta = -0.9E+00

  return
end
subroutine geteps ( eps )
!
!*******************************************************************************
!
!! GETEPS returns the value of EPS, for use by problem 6.
!
  implicit none
!
  real eps
!
  eps = 0.01E+00

  return
end
subroutine getprb ( iprob )
!
!*******************************************************************************
!
!! GETPRB returns the value of the current problem number.
!
  implicit none
!
  integer iprob
!
  iprob = 6

  return
end
subroutine init ( ibc, n, tol, ul, ur, xl, xn, xr )
!
!*******************************************************************************
!
!! INIT initializes some parameters that define the problem.
!
  implicit none
!
  integer n
!
  real beta
  real eps
  integer i
  integer ibc
  integer iprob
  real tol
  real uexact
  real ul
  real ur
  real xl
  real xn(0:n)
  real xr
!
  tol = 0.01E+00
!
!  Find out which problem we're working on.
!
  call getprb ( iprob )
!
!  Set the boundary conditions for the problem, and
!  print out its title.
!
  if ( iprob == 1 ) then

    ibc = 3
    ul = 0.0
    ur = 1.0
    xl = 0.0
    xr = 1.0
    print *, ' '
    print *, 'Exact solution is U = X'

  else if ( iprob == 2 ) then

    ibc = 3
    ul = 0.0
    ur = 1.0
    xl = 0.0
    xr = 1.0
    print *, ' '
    print *, 'Exact solution is U = X*X'

  else if ( iprob == 3 ) then

    ibc = 3
    ul = 0.0
    ur = 1.0
    xl = 0.0
    xr = 1.0
    print *, ' '
    print *, 'Exact solution is U = SIN(PI*X/2)'

  else if ( iprob == 4 ) then

    ibc = 3
    ul = 1.0
    ur = 0.0
    xl = 0.0
    xr = 1.0
    print *, ' '
    print *, 'Exact solution is U = COS(PI*X/2)'

  else if ( iprob == 5 ) then

    ibc = 3
    beta = -0.9
    call getbet(beta)
    ul = 0.0
    ur = 1.0 / ( (beta+2.0) * (beta+1.0) )
    xl = 0.0
    xr = 1.0
    print *, ' '
    print *, 'Rheinboldt problem'

  else if ( iprob == 6 ) then

    ibc = 3
    call geteps(eps)
    xl = 0.0
    xr = 1.0
    ul = uexact(xl)
    ur = uexact(xr)
    print *, ' '
    print *, 'Arctangent problem'

  end if
!
!  The nodes are defined here, and not in the GEOMET routine.
!  This is because each new iteration chooses the location
!  of the new nodes in a special way.
!
  do i = 0, n
    xn(i) = ( real ( n - i ) * xl + real ( i ) * xr ) / real ( n )
  end do

  print *, 'The equation is to be solved for '
  print *, 'X greater than ',xl,' and less than ',xr
  print *, ' '
  print *, 'The boundary conditions are:'
  print *, ' '

  if ( ibc == 1 .or. ibc == 3 ) then
    print *, '  At X = XL, U=',ul
  else
    print *, '  At X = XL, U''=',ul
  end if

  if ( ibc == 2 .or. ibc == 3 ) then
    print *, '  At X = XR, U=',ur
  else
    print *, '  At X = XR, U''=',ur
  end if

  return
end
subroutine output ( f, ibc, indx, n, nu, ul, ur, xn )
!
!*******************************************************************************
!
!! OUTPUT prints out the computed solution.
!
  implicit none
!
  integer n
  integer nu
!
  real error
  real f(nu)
  integer i
  integer ibc
  integer indx(0:n)
  real u
  real uex
  real uexact
  real ul
  real ur
  real xn(0:n)
!
  print *, ' '
  print *, 'Node    X(I)        U(X(I))        Uexact    ' // '    Error'
  print *, ' '

  do i = 0, n

    if ( i == 0 ) then

      if ( ibc == 1 .or. ibc == 3 ) then
        u = ul
      else
        u = f(indx(i))
      end if

    else if ( i == n ) then

      if ( ibc == 2 .or. ibc == 3 ) then
        u = ur
      else
        u = f(indx(i))
      end if

    else

      u = f(indx(i))

    end if

    uex = uexact(xn(i))
    error = u-uex

    write(*,'(i4,4g14.6)') i, xn(i), u, uex, error

  end do

  return
end
subroutine phi ( il, x, phii, phiix, xleft, xrite )
!
!*******************************************************************************
!
!! PHI evaluates a linear basis function and its derivative.
!
!
!  The functions are evaluated at a point X in an interval.  In any
!  interval, there are just two basis functions.  The first
!  basis function is a line which is 1 at the left endpoint
!  and 0 at the right.  The second basis function is 0 at
!  the left endpoint and 1 at the right.
!
  implicit none
!
  integer il
  real phii
  real phiix
  real x
  real xleft
  real xrite
!
  if ( xleft <= x .and. x <= xrite ) then

    if ( il == 1 ) then
      phii = (xrite-x) / (xrite-xleft)
      phiix = -1.0 / (xrite-xleft)
    else
      phii = (x-xleft) / (xrite-xleft)
      phiix = 1.0 / (xrite-xleft)
    end if
!
!  If X is outside of the interval, then the basis function
!  is always zero.
!
  else
    phii = 0.0
    phiix = 0.0
  end if

  return
end
function pp ( x )
!
!*******************************************************************************
!
!! PP evaluates the function P in the differential equation:
!
!    - d/dx (p du/dx) + q u  =  f
!
!  at the point X.
!
  implicit none
!
  integer iprob
  real pp
  real x
!
!  Find out which problem we're working on.
!
  call getprb ( iprob )

  if ( iprob == 1 ) then
    pp = 1.0E+00
  else if ( iprob == 2 ) then
    pp = 1.0E+00
  else if ( iprob == 3 ) then
    pp = 1.0E+00
  else if ( iprob == 4 ) then
    pp = 1.0E+00
  else if ( iprob == 5 ) then
    pp = 1.0E+00
  else if ( iprob == 6 ) then
    pp = 1.0E+00
  end if

  return
end
subroutine prsys ( adiag, aleft, arite, f, nu )
!
!*******************************************************************************
!
!! PRSYS prints out the tridiagonal linear system.
!
  implicit none
!
  integer nu
!
  real adiag(nu)
  real aleft(nu)
  real arite(nu-1)
  real f(nu)
  integer i
!
  print *, ' '
  print *, 'Printout of tridiagonal linear system:'
  print *, ' '
  print *, 'Equation  A-Left  A-Diag  A-Rite  RHS'
  print *, ' '

  do i = 1,nu
    if ( i == 1 ) then
      write(*,'(i3,14x,3g14.6)')i,adiag(i),arite(i),f(i)
    else if ( i < nu ) then
      write(*,'(i3,4g14.6)')i,aleft(i),adiag(i),arite(i),f(i)
    else
      write(*,'(i3,2g14.6,14x,g14.6)')i,aleft(i),adiag(i),f(i)
    end if
  end do

  return
end
function qq ( x )
!
!*******************************************************************************
!
!! QQ evaluates the function Q in the differential equation.
!
!    - d/dx (p du/dx) + q u  =  f
!
!  at the point X.
!
  implicit none
!
  integer iprob
  real qq
  real x
!
!  Find out which problem we're working on.
!
  call getprb(iprob)
!
  if ( iprob == 1 ) then
    qq = 0.0
  else if ( iprob == 2 ) then
    qq = 0.0
  else if ( iprob == 3 ) then
    qq = 0.0
  else if ( iprob == 4 ) then
    qq = 0.0
  else if ( iprob == 5 ) then
    qq = 1.0
  else if ( iprob == 6 ) then
    qq = 0.0
  end if

  return
end
subroutine solve ( adiag, aleft, arite, f, nu )
!
!*******************************************************************************
!
!! SOLVE solves a tridiagonal matrix system of the form A*x = b.
!
!
!  Parameters:
!
!    Input/output, REAL ADIAG(NU), ALEFT(NU), ARITE(NU).
!
!    On input, ADIAG, ALEFT, and ARITE contain the diagonal,
!    left and right entries of the equations.
!
!    On output, ADIAG and ARITE have been changed in order
!    to compute the solution.
!
!    Note that for the first equation, there is no ALEFT
!    coefficient, and for the last, there is no ARITE.
!    So there is no need to store a value in ALEFT(1), nor
!    in ARITE(NU).
!
!    Input/output, REAL F(NU).
!    On input, F contains the right hand side of the linear
!    system to be solve
!    On output, F contains the solution of the linear system.
!
!    Input, INTEGER NU, the number of equations to be solved.
!
  implicit none
!
  integer nu
!
  real adiag(nu)
  real aleft(nu)
  real arite(*)
  real f(nu)
  integer i
!
!  Handle the special case of a single equation.
!
  if ( nu == 1 ) then
    f(1) = f(1) / adiag(1)
!
!  The general case, when NU is greater than 1.
!
  else
    arite(1) = arite(1) / adiag(1)
    do i = 2, nu-1
      adiag(i) = adiag(i) - aleft(i) * arite(i-1)
      arite(i) = arite(i) / adiag(i)
    end do
    adiag(nu) = adiag(nu)-aleft(nu)*arite(nu-1)

    f(1) = f(1)/adiag(1)
    do i = 2,nu
      f(i) = (f(i)-aleft(i)*f(i-1))/adiag(i)
    end do

    do i = nu-1,1,-1
      f(i) = f(i)-arite(i)*f(i+1)
    end do

  end if

  return
end
subroutine solvex ( adiag, aleft, arite, f, h, ibc, indx, kount, n, nl, &
  nmax, node, nquad, nu, ul, ur, wquad, xn, xquad )
!
!*******************************************************************************
!
!! SOLVEX discretizes and solves a differential equation given the nodes.
!
!
  implicit none
!
  integer n
  integer nl
  integer nmax
  integer nquad
!
  real adiag(nmax)
  real aleft(nmax)
  real arite(nmax)
  real f(nmax)
  real h(n)
  integer ibc
  integer indx(0:n)
  integer kount
  integer node(nmax,nl)
  integer nu
  real ul
  real ur
  real wquad(nquad)
  real xn(0:n)
  real xquad(nmax,nquad)
!
!  Given a set of N nodes (where N increases on each iteration),
!  compute the other geometric information.
!
  call geomet ( h, ibc, indx, n, nl, nmax, node, nquad, nu, wquad, xn, xquad )
!
!  Assemble the linear system.
!
  call assemb ( adiag, aleft, arite, f, h, n, indx, node, nu, nl, &
    nquad, nmax, ul, ur, wquad, xn, xquad )
!
!  Print out the linear system, just once.
!
  if ( kount == 1 ) then
    call prsys ( adiag, aleft, arite, f, nu )
  end if
!
!  Solve the linear system.
!
  call solve ( adiag, aleft, arite, f, nu )
!
!  Print out the solution.
!
  print *, ' '
  print *, 'Basic solution'

  call output ( f, ibc, indx, n, nu, ul, ur, xn )

  return
end
subroutine solvey ( eta, f, h, n, nu, ul, ur, xn )
!
!*******************************************************************************
!
!! SOLVEY computes error estimators for a finite element solution.
!
!
!  Discussion:
!
!    SOLVEY accepts information about the solution of a finite element
!    problem on a grid of nodes with coordinates XN.  It then starts
!    at node 0, and for each node, computes two "error estimators",
!    one for the left, and one for the right interval associated with the
!    node.  These estimators are found by solving a finite element problem
!    over the two intervals, using the known values of the original
!    solution as boundary data, and using a mesh that is "slightly"
!    refined over the original one.
!
!    Note that the computations at the 0-th and N-th nodes only involve
!    a single interval.
!
  implicit none
!
  integer, parameter :: nl = 2
  integer, parameter :: ny = 2
  integer, parameter :: nquad = 2
!
  integer, parameter :: nmay = 2*ny
!
  integer n
  integer nu
!
  real adiag(nmay)
  real aleft(nmay)
  real arite(nmay)
  real eta(n)
  real f(nu)
  real fy(nmay)
  real h(n)
  real hy(nmay)
  integer i
  integer ibcy
  integer indy(0:nmay)
  integer j
  integer jhi
  integer jlo
  integer jmid
  integer k
  integer m
  integer nodey(nmay,nl)
  integer nuy
  real pp
  real qq
  real sum
  real ul
  real uleft
  real ulval
  real uly
  real uprime
  real ur
  real urite
  real urval
  real ury
  real uval
  real vlval
  real vprime
  real vrval
  real vval
  real wquad(nquad)
  real xn(0:n)
  real xquady(nmay,nquad)
  real y
  real yl
  real ym
  real yn(0:nmay)
  real yr
!
!  Initialize the error estimators to zero.
!
  eta(1:n) = 0.0
!
!  Set the boundary conditions for each subproblem to be
!  known values of U at the left and right.
!
!
!  For each node, subdivide its left and right hand intervals
!  into NY subintervals.
!
!  Set up and solve the differential equation again on this
!  smaller region.
!
!  The 0-th and N-th nodes are special cases.
!
  ibcy = 3

  do j = 0, n

    if ( j == 0 ) then
      m = ny
      jlo = j
      jmid = j+1
      jhi = j+1
    else if ( j == n ) then
      m = ny
      jlo = j-1
      jmid = j
      jhi = j
    else
      m = 2*ny
      jlo = j-1
      jmid = j
      jhi = j+1
    end if
!
!  Set the location of the nodes in the subintervals.
!
    yl = xn(jlo)
    ym = xn(jmid)
    yr = xn(jhi)

    do i = 0, ny
      yn(i) = ((ny-i)*yl+i*ym)/real(ny)
    end do

    do i = ny+1, m
      yn(i) = ((m-i)*ym+(i-ny)*yr)/real(ny)
    end do
!
!  Set up the geometry of the sub-problem.
!
    call geomet(hy,ibcy,indy,m,nl,nmay,nodey,nquad,nuy,wquad,yn,xquady)
!
!  Set the boundary values for the sub-problem.
!
    if ( j <= 1 ) then
      uly = ul
    else
      uly = f(j-1)
    end if

    if ( j >= n-1 ) then
      ury = ur
    else
      ury = f(j+1)
    end if
!
!  Assemble the matrix for the sub-problem.
!
    call assemb(adiag,aleft,arite,fy,hy,m,indy,nodey,nuy,nl, &
      nquad,nmay,uly,ury,wquad,yn,xquady)
!
!  Solve the system.
!
    call solve(adiag,aleft,arite,fy,nuy)
!
!  Compute the weighted sum of the squares of the differences
!  of the original computed slope and the refined computed slopes.
!
!  Calculation for left interval.
!
    if ( j >= 1 ) then

      if ( j <= 1 ) then
        uleft = ul
        urite = f(1)
      else if ( j == n ) then
        uleft = f(j-1)
        urite = ur
      else
        uleft = f(j-1)
        urite = f(j)
      end if

      uprime = (urite-uleft)/h(j)

      sum = 0.0
      do i = 1,ny

        yl = yn(i-1)
        yr = yn(i)

        if ( i == 1 ) then
          vlval = uly
          vrval = fy(i)
        else if ( i == m ) then
          vlval = fy(i-1)
          vrval = ury
        else
          vlval = fy(i-1)
          vrval = fy(i)
        end if

        vprime = (vrval-vlval)/hy(i)

        ulval = ((ny+1-i)*uleft+(i-1)*urite)/real(ny)
        urval = ((ny-i)*uleft+i*urite)/real(ny)
!
!  Compute the integral of
!
!    p(x)*(u'(x)-v'(x))**2 + q(x)*(u(x)-v(x))**2
!
        do k = 1,nquad

          y  =  xquady(i,k)

          uval = ( (yl-y)*urval+(y-yr)*ulval)/(yl-yr)
          vval = ( (yl-y)*vrval+(y-yr)*vlval)/(yl-yr)

          sum = sum+0.5*wquad(k)*hy(i) * &
            (pp(y)*(uprime-vprime)**2 &
            +qq(y)*(uval-vval)**2)

        end do

      end do

      eta(j) = eta(j)+0.5*sqrt(sum)

    end if
!
!  Calculation for right interval.
!
    if ( j <= n-1 ) then

      if ( j == 0 ) then
        uleft = ul
        urite = f(j+1)
      else if ( j.ge.n-1 ) then
        uleft = f(j)
        urite = ur
      else
        uleft = f(j)
        urite = f(j+1)
      end if

      uprime = (urite-uleft)/h(j+1)

      sum = 0.0
      do i = m+1-ny, m

        yl = yn(i-1)
        yr = yn(i)

        if ( i == 1 ) then
          vlval = uly
          vrval = fy(i)
        else if ( i == m ) then
          vlval = fy(i-1)
          vrval = ury
        else
          vlval = fy(i-1)
          vrval = fy(i)
        end if

        vprime = (vrval-vlval)/hy(i)

        ulval = ((m+1-i)*uleft+(i+ny-m-1)*urite)/real(ny)
        urval = ((m-i)*uleft+(i+ny-m)*urite)/real(ny)
!
!  Compute the integral of
!
!    p(x)*(u'(x)-v'(x))**2 + q(x)*(u(x)-v(x))**2
!

        do k = 1,nquad

          y  =  xquady(i,k)

          uval = ( (yl-y)*urval+(y-yr)*ulval)/(yl-yr)
          vval = ( (yl-y)*vrval+(y-yr)*vlval)/(yl-yr)

          sum = sum+0.5*wquad(k)*hy(i) * &
            (pp(y)*(uprime-vprime)**2 &
            +qq(y)*(uval-vval)**2)

        end do

      end do

      eta(j+1) = eta(j+1) + 0.5 * sqrt ( sum )

    end if

  end do
!
!  Print out the error estimators.
!
  print *, ' '
  print *, 'ETA'
  print *, ' '
  do j = 1, n
    print *, eta(j)
  end do

  return
end
subroutine subdiv ( eta, jadd, kount, n, nmax, tol, xn, xt )
!
!*******************************************************************************
!
!! SUBDIV decides which intervals should be subdivided.
!
  implicit none
!
  integer n
!
  real ave
  real eta(n)
  integer j
  integer jadd(n)
  integer k
  integer kount
  integer nmax
  real sum
  real temp
  real tol
  real xn(0:nmax)
  real xt(0:nmax)
!
!  Add up the ETA's, and get their average.
!
  sum = 0.0
  do j = 1, n
    sum = sum+eta(j)
  end do
  ave = sum / real(n)
!
!  Look for intervals whose ETA value is relatively large,
!  and note in JADD that these intervals should be subdivided.
!
  k = 0
  temp = max ( 1.2E+00 * ave + 0.00001, tol**2 / real ( n ) )

  print *, ' '
  print *, 'Tolerance = ',temp
  print *, ' '

  do j = 1, n

    if ( eta(j) > temp ) then
      k = k+1
      jadd(j) = 1
      print *, 'Subdivide interval ',j
    else
      jadd(j) = 0
    end if

  end do
!
!  If no subdivisions needed, we're done.
!
  if ( k <= 0 ) then
    print *, 'Success on step ',kount
    stop
  end if
!
!  See if we're about to go over our limit.
!
  if ( n+k > nmax ) then
    print *, ' '
    print *, 'The iterations did not reach their goal.'
    print *, 'The next value of N is ',n+k
    print *, 'which exceeds NMAX = ',nmax
    stop
  end if
!
!  Insert new nodes where needed.
!
  k = 0
  xt(0) = xn(0)
  do j = 1, n

    if ( jadd(j) > 0 ) then
      xt(j+k) = 0.5*(xn(j)+xn(j-1))
      k = k+1
    end if

    xt(j+k) = xn(j)

  end do
!
!  Update the value of N, and copy the new nodes into XN.
!
  n = n+k

  xn(0:n) = xt(0:n)

  return
end
function uexact ( x )
!
!*******************************************************************************
!
!! UEXACT returns the value of the exact solution at any point X.
!
  implicit none
!
  real beta
  real eps
  integer iprob
  real, parameter :: pi = 3.14159265E+00
  real uexact
  real x
!
!  Find out which problem we're working on.
!
  call getprb ( iprob )

  if ( iprob == 1 ) then
    uexact = x
  else if ( iprob == 2 ) then
    uexact = x**2
  else if ( iprob == 3 ) then
    uexact = sin ( pi * x / 2.0 )
  else if ( iprob == 4 ) then
    uexact = cos ( pi * x / 2.0 )
  else if ( iprob == 5 ) then
    call getbet ( beta )
    uexact = ( x**(beta+2) ) / ( ( beta + 2.0 ) * ( beta + 1.0 ) )
  else if ( iprob == 6 ) then
    call geteps ( eps )
    uexact = atan ( ( x - 0.5 ) / eps )
  else
    uexact = 0.0E+00
  end if

  return
end
subroutine timestamp ( )
!
!*******************************************************************************
!
!! TIMESTAMP prints the current YMDHMS date as a time stamp.
!
!
!  Example:
!
!    May 31 2001   9:45:54.872 AM
!
!  Modified:
!
!    31 May 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    None
!
  implicit none
!
  character ( len = 8 ) ampm
  integer d
  character ( len = 8 ) date
  integer h
  integer m
  integer mm
  character ( len = 9 ), parameter, dimension(12) :: month = (/ &
    'January  ', 'February ', 'March    ', 'April    ', &
    'May      ', 'June     ', 'July     ', 'August   ', &
    'September', 'October  ', 'November ', 'December ' /)
  integer n
  integer s
  character ( len = 10 )  time
  integer values(8)
  integer y
  character ( len = 5 ) zone
!
  call date_and_time ( date, time, zone, values )

  y = values(1)
  m = values(2)
  d = values(3)
  h = values(5)
  n = values(6)
  s = values(7)
  mm = values(8)

  if ( h < 12 ) then
    ampm = 'AM'
  else if ( h == 12 ) then
    if ( n == 0 .and. s == 0 ) then
      ampm = 'Noon'
    else
      ampm = 'PM'
    end if
  else
    h = h - 12
    if ( h < 12 ) then
      ampm = 'PM'
    else if ( h == 12 ) then
      if ( n == 0 .and. s == 0 ) then
        ampm = 'Midnight'
      else
        ampm = 'AM'
      end if
    end if
  end if

  write ( *, '(a,1x,i2,1x,i4,2x,i2,a1,i2.2,a1,i2.2,a1,i3.3,1x,a)' ) &
    trim ( month(m) ), d, y, h, ':', n, ':', s, '.', mm, trim ( ampm )

  return
end