c===========================================================
c     Solves 1-d non-linear boundary value problem  
c
c         u''(x) + (u u')^2 + sin(u) = f(x)
c
c         on  x = [0,1]; u(0) = 0, u(1) = 0
c 
c     using second-order finite difference technique, 
c     Newton's method and LAPACK tridiagonal solver DGTSV.
c
c     Currently set-up for solution
c
c         u(x) = sin(4 Pi x)
c===========================================================
      program         nlbvp1d

      implicit        none

      integer         i4arg
      real*8          r8arg,        drelabs,      dvl2norm

      real*8          r8_never
      parameter     ( r8_never = -1.0d-60 )
c-----------------------------------------------------------
c     Extrema of problem domain.
c-----------------------------------------------------------
      real*8          xmin,          xmax
      parameter     ( xmin = 0.0d0,  xmax = 1.0d0 )
c-----------------------------------------------------------
      integer         maxn
      parameter     ( maxn = 32 769 )
c-----------------------------------------------------------
c     Storage for discrete x-values, unknowns, coefficient
c     exact 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 (Jacobian matrix) and 
c     right-hand-side vector (residual vector) for use with 
c     LAPACK routine DGTSV.  Other parameters needed for 
c     call to DGTSV.
c-----------------------------------------------------------
      real*8          d(maxn),       du(maxn),
     &                dl(maxn),      rhs(maxn)
      integer         nrhs 
      parameter     ( nrhs = 1 )
      integer         info

c-----------------------------------------------------------
c     Discretization level and size of system (# of discrete
c     unknowns)
c-----------------------------------------------------------
      integer         level,         n,             i,
     &                option
c-----------------------------------------------------------
c     Variables used in Newton iteration.
c-----------------------------------------------------------
      integer         mxiter
      parameter     ( mxiter = 50 )

      integer         iter
      real*8          guess_factor,                 tol,        
     &                nrm2res,       nrm2du,        nrm2u

c-----------------------------------------------------------
c     Enable following parameter for full tracing of 
c     Newton iteration.
c-----------------------------------------------------------
      logical         newton_trace
      parameter     ( newton_trace = .false. )

c-----------------------------------------------------------
c     Mesh spacing and related constants 
c-----------------------------------------------------------
      real*8          h,             hm2,           m2hm2,
     &                hm1by2,        hhm2,          qhm2

      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
      guess_factor = r8arg(2,r8_never)
      if( guess_factor .eq. r8_never ) go to 900
      option  = i4arg(3,0)
      tol     = r8arg(4,1.0d-8)
c-----------------------------------------------------------
c     Set up finite-difference 'mesh' (discrete x-values)
c     and define some useful constants.
c-----------------------------------------------------------
      h     =  1.0d0 / (n - 1)
      do i = 1 , n
         x(i) = xmin + (i - 1) * h
      end do
      hm2     =  1.0d0 / (h * h)
      m2hm2   = -2.0d0 / (h * h)
      hm1by2  =  0.50d0 / h       
      hhm2    =  0.50d0 * hm2
      qhm2    =  0.25d0 * hm2
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, coefficient functions and right 
c     hand side vector.
c-----------------------------------------------------------
      call exact(uexact,f,x,n)

c-----------------------------------------------------------
c     Initialize unknown (u) to constant (guess_factor) 
c     times exact solution. 
c-----------------------------------------------------------
      do i = 1 , n
         u(i) = guess_factor * uexact(i)
      end do

c===========================================================
c     N E W T O N   L O O P
c===========================================================
      do iter = 1 , mxiter

c-----------------------------------------------------------
c        Set up tridiagonal Jacobian matrix and evaluate 
c        right-hand-side (residuals)
c-----------------------------------------------------------

c-----------------------------------------------------------
c        Left boundary: Dirichlet boundary condition has
c        0 residual.
c-----------------------------------------------------------
         d(1)       = 1.0d0 
         du(1)      = 0.0d0
         rhs(1)     = 0.0d0
c-----------------------------------------------------------
c        Interior: J[i,j] = d(F_i)/d_(u_j) and has non-zero
c        elements only for j = i-1, i and i+1.
c-----------------------------------------------------------
         do i = 2 , n - 1
            dl(i-1) = hm2 - hhm2 * u(i)**2 * (u(i+1) - u(i-1))
            d(i)    = m2hm2 
     &                + hhm2 * u(i) * (u(i+1) - u(i-1))**2
     &                + cos(u(i))
            du(i)   = hm2 + hhm2 * u(i)**2 * (u(i+1) - u(i-1))
            rhs(i)  = hm2 * (u(i+1) - 2.0d0 * u(i) + u(i-1)) 
     &              + qhm2 * u(i)**2 * (u(i+1) - u(i-1))**2
     &              + sin(u(i)) - f(i)
         end do
c-----------------------------------------------------------
c        Right boundary: Dirichlet boundary condition has
c        0 residual.
c-----------------------------------------------------------
         dl(n-1)    = 0.0d0
         d(n)       = 1.0d0
         rhs(n)     = 0.0d0
c-----------------------------------------------------------
c        Compute l-2 norm of residuals.
c-----------------------------------------------------------
         nrm2res = dvl2norm(rhs,n)
         if( newton_trace ) then
            write(0,*) 'iter = ', iter 
            write(0,*) 'res  = ', nrm2res
         end if
c-----------------------------------------------------------
c        Solve tridiagonal system for Newton update, delu,
c        which satisfies 
c
c        J delu = residuals
c-----------------------------------------------------------
      
         call dgtsv( n, nrhs, dl, d, du, rhs, n, info )

         if( info .eq. 0 ) then
c-----------------------------------------------------------
c           Solver successful: compute norms of u and delu,
c           update solution and check for convergence.
c-----------------------------------------------------------
            nrm2u  = dvl2norm(u,n)   
            nrm2du = dvl2norm(rhs,n)
            if( newton_trace ) then
               write(0,*) 'du   = ', nrm2du
               write(0,*) 'u    = ', nrm2u
            end if
            do i = 1 , n
               u(i) = u(i) - rhs(i)
            end do
            if( drelabs(nrm2du,nrm2u,1.0d-6) .le. tol ) 
     &          go to 500
         else
c-----------------------------------------------------------
c           Solver failed, write error message and exit.
c-----------------------------------------------------------
            write(0,*) 'nlbvp1d: dgtsv() failed, info = ', 
     &                 info
         end if
      end do
c-----------------------------------------------------------
c     Newton iteration failed to converge: write error 
c     message and exit.
c-----------------------------------------------------------
      write(0,*) 'nlbvp1d: No convergence after ', mxiter,  
     &           ' iterations'
      stop

c-----------------------------------------------------------
c     Newton iteration converged: output solution or error 
c     to stdout, depending on output option.  Also compute 
c     rms error and output to stderr.
c-----------------------------------------------------------
500   continue

      rmserr = 0.0d0
      do i = 1 , n
         if( option .eq. 0 ) then
            write(*,*) x(i), u(i)
         else 
            write(*,*) x(i), (uexact(i) - u(i))
         end if
         rmserr = rmserr + (uexact(i) - u(i)) ** 2
      end do
      rmserr = sqrt(rmserr / n)
      write(0,*) 'rmserr = ', rmserr

      stop

 900  continue
         write(0,*) 'usage: nlbvp1d <level> <guess_factor> '//
     &              '[<option> <tol>]'
         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          pi4
         integer         i

         pi4 = 16.0d0 * atan(1.0d0)
         do i = 1 , n
            u(i) = sin(pi4 * x(i))
            f(i) = -pi4**2 * sin(pi4 * x(i)) + 
     &              pi4**2 * (sin(pi4 * x(i)) * 
     &                        cos(pi4 * x(i)))**2 +
     &              sin(sin(pi4 * x(i)))
         end do

         return

      end

c===========================================================
c     dvl2norm:  Returns l2-norm of double precision vector.
c===========================================================
      real*8 function dvl2norm(v,n)

         implicit      none

         integer       n
         real*8        v(n)
         integer       i

         dvl2norm = 0.0d0
         do i = 1 , n
            dvl2norm = dvl2norm + v(i) * v(i)
         end do
         if( n .gt. 0 ) then
            dvl2norm = sqrt(dvl2norm / n)
         end if

         return

      end

c===========================================================
c     drelabs:  Function useful for 'relativizing' quantity
c     being monitored for detection of convergence.  
c===========================================================
      real*8 function drelabs(dx,x,xfloor)

         implicit      none

         real*8        dx,     x,    xfloor

         if( abs(x) .lt. abs(xfloor) ) then
            drelabs = abs(dx)
         else
            drelabs = abs(dx/x)
         end if

         return

      end
