subroutine ctof ( nc, uc, nf, uf )

c*********************************************************************72
c                                                    
cc ctof() transfers data from a coarse to a finer grid.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    25 November 2011
c
c  Author:
c
c    John Burkardt
c
c  Reference:
c
c    William Hager,
c    Applied Numerical Linear Algebra,
c    Prentice-Hall, 1988,
c    ISBN13: 978-0130412942,
c    LC: QA184.H33.
c
c  Input:
c
c    integer NC, the number of coarse nodes.
c
c    double precision UC(NC), the coarse correction data.
c
c    integer NF, the number of fine nodes.
c
c    double precision UF(NF): the fine grid data.
c
c  Output:
c
c    double precision UF(NF): the updated fine grid data.
c
      implicit none

      integer nc
      integer nf

      integer ic
      integer if
      double precision uc(nc)
      double precision uf(nf)

      do ic = 1, nc
        if = 2 * ic - 1
        uf(if) = uf(if) + uc(ic)
      end do

      do ic = 1, nc - 1
        if = 2 * ic
        uf(if) = uf(if) + 0.5D+00 * ( uc(ic) + uc(ic+1) )
      end do

      return
      end
      subroutine ftoc ( nf, uf, rf, nc, uc, rc )

c*********************************************************************72
c                                                    
cc ftoc() transfers data from a fine grid to a coarser grid.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    07 December 2011
c
c  Author:
c
c    John Burkardt
c
c  Reference:
c
c    William Hager,
c    Applied Numerical Linear Algebra,
c    Prentice-Hall, 1988,
c    ISBN13: 978-0130412942,
c    LC: QA184.H33.
c
c  Input:
c
c    integer NF, the number of fine nodes.
c
c    double precision UF(NF), the fine data.
c
c    double precision RF(NF), the right hand side for the fine grid.
c
c    integer NC, the number of coarse nodes.
c
c  Output:
c
c    double precision UC(NC), the coarse grid data, set to zero.
c
c    double precision RC(NC), the right hand side for the coarse grid.
c
      implicit none

      integer nc
      integer nf

      integer ic
      integer if
      double precision rc(nc)
      double precision rf(nf)
      double precision uc(nc)
      double precision uf(nf)

      uc(1) = 0.0D+00
c     rc(1) = rf(1) - uf(1)
      rc(1) = 0.0D+00
      do ic = 2, nc - 1
        if = 2 * ic - 1
        uc(ic) = 0.0D+00
        rc(ic) = 4.0D+00 
     &    * ( rf(if) + uf(if-1) - 2.0D+00 * uf(if) + uf(if+1) )
      end do
      uc(nc) = 0.0D+00
c     rc(nc) = rf(nc) - uf(nc)
      rc(nc) = 0.0D+00

      return
      end
      subroutine gauss_seidel ( n, r, u, dif_l1 )

c*********************************************************************72
c                                                    
cc gauss_seidel() carries out one step of a Gauss-Seidel iteration.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    25 November 2011
c
c  Author:
c
c    John Burkardt
c
c  Reference:
c
c    William Hager,
c    Applied Numerical Linear Algebra,
c    Prentice-Hall, 1988,
c    ISBN13: 978-0130412942,
c    LC: QA184.H33.
c
c  Input:
c
c    integer N, the number of unknowns.
c
c    double precision R(N), the right hand side.
c
c    double precision U(N), the estimated solution.
c
c  Output:
c
c    double precision U(N), the updated estimated solution.
c
c    double precision DIF_L1, the L1 norm of the difference between the
c    input and output solution estimates.
c
      implicit none

      integer n

      double precision dif_l1
      integer i
      double precision r(n)
      double precision u(n)
      double precision u_old

      dif_l1 = 0.0D+00
c
c  Setting U(1) = R(1), U(N) = R(N) seems right, but it makes the
c  code fail.
c
c     u(1) = r(1)
      do i = 2, n - 1
        u_old = u(i)
        u(i) = 0.5D+00 * ( u(i-1) + u(i+1) + r(i) )
        dif_l1 = dif_l1 + abs ( u(i) - u_old )
      end do
c     u(n) = r(n)

      return
      end
      function i4_log_2 ( i )

c*********************************************************************72
c
cc i4_log_2() returns the integer part of the logarithm base 2 of an I4.
c
c  Discussion:
c
c    For positive I4_LOG_2(I), it should be true that
c      2^I4_LOG_2(X) .le. |I| < 2^(I4_LOG_2(I)+1).
c    The special case of I4_LOG_2(0) returns -HUGE().
c
c    An I4 is an integer value.
c
c  Example:
c
c     I  I4_LOG_2
c
c     0  -1
c     1,  0
c     2,  1
c     3,  1
c     4,  2
c     5,  2
c     6,  2
c     7,  2
c     8,  3
c     9,  3
c    10,  3
c   127,  6
c   128,  7
c   129,  7
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    28 October 2007
c
c  Author:
c
c    John Burkardt
c
c  Input:
c
c    integer I, the number whose logarithm base 2
c    is desired.
c
c  Output:
c
c    integer I4_LOG_2, the integer part of the
c    logarithm base 2 of the absolute value of I.
c
      implicit none

      integer i
      integer i_abs
      integer i4_log_2
      integer i4_huge
      parameter ( i4_huge = 2147483647 )

      if ( i .eq. 0 ) then

        i4_log_2 = - i4_huge

      else

        i4_log_2 = 0

        i_abs = abs ( i )

10      continue

        if ( 2 .le. i_abs ) then
          i_abs = i_abs / 2
          i4_log_2 = i4_log_2 + 1
          go to 10
        end if

      end if

      return
      end
      subroutine poisson_1d_multigrid ( n, a, b, ua, ub, force, 
     &  it_num, u )

c*********************************************************************72
c                                                    
cc poisson_1d_multigrid() solves a 1D PDE using the multigrid method.
c
c  Discussion:
c
c    This routine solves a 1D boundary value problem of the form
c
c      - U''(X) = F(X) for A < X < B,
c
c    with boundary conditions U(A) = UA, U(B) = UB.
c
c    The multigrid method is used. 
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    26 July 2014
c
c  Author:
c
c    Original Fortran77 version by William Hager.
c    This version by John Burkardt.
c
c  Reference:
c
c    William Hager,
c    Applied Numerical Linear Algebra,
c    Prentice-Hall, 1988,
c    ISBN13: 978-0130412942,
c    LC: QA184.H33.
c
c  Input:
c
c    integer N, the number of intervals.
c    N must be a power of 2.
c
c    double precision A, B, the endpoints.
c
c    double precision UA, UB, the boundary values at the endpoints.
c
c    double precision external FORCE, the name of the function 
c    which evaluates the right hand side.
c
c  Output:
c
c    integer IT_NUM, the number of iterations.
c
c    integer U(N+1), the computed solution.
c
      implicit none

      integer n
      integer nl_max
      parameter ( nl_max = 518 )

      double precision a
      double precision b
      double precision d0
      double precision d1
      external force
      double precision force
      double precision h
      integer i
      integer i4_log_2
      integer it
      integer it_num
      integer j
      integer k
      integer l
      integer ll
      integer m
      integer nl
      double precision r(nl_max)
      double precision tol
      double precision u(n+1)
      double precision ua
      double precision ub
      double precision utol
      double precision uu(nl_max)
      double precision x(n+1)
c
c  Determine if we have enough storage.
c
      k = i4_log_2 ( n )

      if ( n .ne. 2**k ) then
        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) 'poisson_1d_multigrid(): Fatal error!'
        write ( *, '(a)' ) '  Interval number N is not a power of 2.'
        stop 1
      end if

      nl = n + n + k - 2

      if ( nl_max < nl ) then
        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) 'poisson_1d_multigrid(): Fatal error!'
        write ( *, '(a)' ) '  Grid index K requires NL = ', nl
        write ( *, '(a)' ) '  Internal parameter NL_MAX = ', nl_max
        stop 1
      end if
c
c  Initialization.
c
      it = 4
      it_num = 0
      tol = 0.0001D+00
      utol = 0.7D+00
      m = n
c
c  Set the nodes.
c
      call r8vec_linspace ( n + 1, a, b, x )
c 
c  Set the right hand side.
c 
      r(1) = ua
      h = ( b - a ) / dble ( n )
      do i = 2, n
        r(i) = h**2 * force ( x(i) )
      end do
      r(n+1) = ub

      do i = 1, nl
        uu(i) = 0.0D+00
      end do
c
c  L points to first entry of solution
c  LL points to penultimate entry.
c
      l = 1
      ll = n
c 
c  Gauss-Seidel iteration
c
      d1 = 0.0D+00
      j = 0

10    continue

        d0 = d1
        j = j + 1
        call gauss_seidel ( n + 1, r(l), uu(l), d1 )
        it_num = it_num + 1
!
!  Must do at least 4 iterations at each level.
!
        if ( j .lt. it ) then

          go to 10
!
!  If enough iterations, and satisfactory decrease, and
!  on finest grid, exit.
!
        else if ( d1 .lt. tol .and. n .eq. m ) then

          go to 20
!
!  Enough iterations, satisfactory convergence, go finer.
!
        else if ( d1 .lt. tol ) then

          call ctof ( n + 1, uu(l), n + n + 1, uu(l-1-n-n) )

          n = n + n
          ll = l - 2
          l = l - 1 - n
          j = 0
!
!  Enough iterations, slow convergence, and 2 < n, go coarser.
!
        else if ( utol * d0 .le. d1 .and. 2 .lt. n ) then

          call ftoc ( n + 1, uu(l), r(l), (n/2)+1, uu(l+n+1), r(l+n+1) )

          n = n / 2
          l = ll + 2
          ll = ll + n + 1
          j = 0

        end if

      go to 10
c
c  Computation complete
c
20    continue

      do i = 1, n + 1
        u(i) = uu(i)
      end do

      return
      end
      subroutine r8vec_linspace ( n, a, b, x )

c*********************************************************************72
c
cc r8vec_linspace() creates a vector of linearly spaced values.
c
c  Discussion:
c
c    An R8VEC is a vector of R8's.
c
c    4 points evenly spaced between 0 and 12 will yield 0, 4, 8, 12.
c
c    In other words, the interval is divided into N-1 even subintervals,
c    and the endpoints of intervals are used as the points.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    14 March 2011
c
c  Author:
c
c    John Burkardt
c
c  Input:
c
c    integer N, the number of entries in the vector.
c
c    double precision A, B, the first and last entries.
c
c  Output:
c
c    double precision X(N), a vector of linearly spaced data.
c
      implicit none

      integer n

      double precision a
      double precision b
      integer i
      double precision x(n)

      if ( n .eq. 1 ) then

        x(1) = ( a + b ) / 2.0D+00

      else

        do i = 1, n
          x(i) = ( dble ( n - i     ) * a
     &           + dble (     i - 1 ) * b )
     &           / dble ( n     - 1 )
        end do

      end if

      return
      end