c===========================================================
c     Solves 1-d linear boundary value problem  
c
c         u''(x) = f(x)  on  x = [0,1]; u(0) = u0, u(1) = u1
c 
c     using second-order finite difference technique and 
c     LAPACK tridiagonal solver DGTSV.
c===========================================================
      program         bvp1d

      implicit        none

      integer         i4arg

c-----------------------------------------------------------
c     Extrema of problem domain; note that this approach 
c     of defining extrema as parameters makes it easier 
c     to generalize program to arbitrary domains.
c-----------------------------------------------------------
      real*8          xmin,          xmax
      parameter     ( xmin = 0.0d0,  xmax = 1.0d0 )
c-----------------------------------------------------------
c     Maximum problem size (2**16 + 1)
c-----------------------------------------------------------
      integer         maxn
      parameter     ( maxn = 65 537 )
c-----------------------------------------------------------
c     Storage for discrete x-values, unknowns, exact 
c     solution and right hand side values
c-----------------------------------------------------------
      real*8          x(maxn),       u(maxn),
     &                uexact(maxn),  f(maxn)
c-----------------------------------------------------------
c     Storage for main, upper and lower diagonals of 
c     tridiagonal system, and right-hand-side vector 
c     for use with LAPACK routine DGTSV
c-----------------------------------------------------------
      real*8          d(maxn),       du(maxn),
     &                dl(maxn),      rhs(maxn)
      integer         nrhs,          info
c-----------------------------------------------------------
c     Discretization level and size of system (# of discrete
c     unknowns)
c-----------------------------------------------------------
      integer         level,         n,             j,
     &                option
c-----------------------------------------------------------
c     Mesh spacing and related constants (1/h**2, -2/h**2)
c-----------------------------------------------------------
      real*8          h,             hm2,           m2hm2
      real*8          rmserr
c-----------------------------------------------------------
c     Argument parsing.
c-----------------------------------------------------------
      level = i4arg(1,-1)
      if( level .lt. 0 ) go to 900
      n = 2 ** level + 1
      if( n .gt. maxn ) then
         write(0,*) 'Insufficient internal storage'
         stop
      end if
      option = i4arg(2,0)
c-----------------------------------------------------------
c     Set up finite-difference 'mesh' (discrete x-values)
c     and define some useful constants.
c-----------------------------------------------------------
      h     =  1.0d0 / (n - 1)
      do j = 1 , n
         x(j) = xmin + (j - 1) * h
      end do
      hm2   =  1.0d0 / (h * h)
      m2hm2 = -2.0d0 / (h * h)
c-----------------------------------------------------------
c     This only ensures that x(n) = xmax EXACTLY and is not 
c     essential.
c-----------------------------------------------------------
      x(n) = xmax
c-----------------------------------------------------------
c     Set up exact solution and right hand side vector.
c-----------------------------------------------------------
      call exact(uexact,f,x,n)

c===========================================================
c     Set up tridiagonal system.  Note that indexing on 
c     lower diagonal is always (j-1) when implementing the 
c     j'th equation.
c===========================================================

c-----------------------------------------------------------
c     Left boundary: u(1) = u_0
c-----------------------------------------------------------
      d(1)       = 1.0d0 
      du(1)      = 0.0d0
      rhs(1)     = uexact(1)
c-----------------------------------------------------------
c     Interior: Second order FDA of ODE.
c-----------------------------------------------------------
      do j = 2 , n - 1
         dl(j-1) = hm2
         d(j)    = m2hm2
         du(j)   = hm2
         rhs(j)  = f(j)
      end do
c-----------------------------------------------------------
c     Right boundary: u(n) = u_1
c-----------------------------------------------------------
      dl(n-1)    = 0.0d0
      d(n)       = 1.0d0
      rhs(n)     = uexact(n)

c===========================================================
c     Solve tridiagonal system.
c===========================================================
      
      nrhs = 1
      call dgtsv( n, nrhs, dl, d, du, rhs, n, info )

      if( info .eq. 0 ) then
c-----------------------------------------------------------
c        Solver successful, output either (x_j, u_j) or 
c        (x_j, error_j) to stdout.  Also compute rms error 
c        and output to standard error.
c-----------------------------------------------------------
         rmserr = 0.0d0
         do j = 1 , n
            if( option .eq. 0 ) then
               write(*,*) x(j), rhs(j)
            else 
               write(*,*) x(j), (uexact(j) - rhs(j))
            end if
            rmserr = rmserr + (uexact(j) - rhs(j)) ** 2
         end do
         rmserr = sqrt(rmserr / n)
         write(0,*) 'rmserr = ', rmserr
      else
c-----------------------------------------------------------
c        Solver failed.
c-----------------------------------------------------------
         write(0,*) 'bvp1d: dgtsv() failed, info = ', info
      end if

      stop

 900  continue
         write(0,*) 'usage: bvp1d <level> [<option>]'
         write(0,*) 
         write(0,*) '       Specify option .ne. 0 for output'
         write(0,*) '       of error instead of solution'
      stop
   
      end

c===========================================================
c     Computes exact values for u(x) (unknown function)
c     and f(x) (right hand side function).  x array must 
c     have been previously defined.
c===========================================================
      subroutine exact(u,f,x,n)

         implicit        none
         integer         n
         real*8          u(n),     f(n),    x(n)

         real*8          pi2
         integer         j

         pi2 = 8.0d0 * atan(1.0d0)
         do j = 1 , n
            u(j) = sin(pi2 * x(j))
            f(j) = -pi2 * pi2 * u(j)
         end do

         return

      end