program main

c*********************************************************************72
c
cc poisson_1d_multigrid_test() tests poisson_1d_multigrid().
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
      implicit none

      call timestamp ( )
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'poisson_1d_multigrid_test():'
      write ( *, '(a)' ) '  FORTRAN77 version'
      write ( *, '(a)' ) '  Test poisson_1d_multigrid().'

      call test01_mono ( ) 
      call test01_multi ( ) 
      call test02_mono ( ) 
      call test02_multi ( ) 
c
c  Terminate.
c
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'poisson_1d_multigrid_test():'
      write ( *, '(a)' ) '  Normal end of execution.'
      write ( *, '(a)' ) ' '
      call timestamp ( )

      stop
      end
      subroutine test01_mono ( ) 

c*********************************************************************72
c
cc test01_mono() tests poisson_1d_monogrid() on test case 1.
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    John Burkardt
c
      implicit none

      integer n_max
      parameter ( n_max = 32 )

      double precision a
      double precision b
      double precision difmax
      external exact1
      double precision exact1
      external force1
      double precision force1
      integer i
      integer it_num
      integer k
      integer n
      double precision u(n_max+1)
      double precision ua
      double precision ub
      double precision x(n_max+1)

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'test01_mono():'
      write ( *, '(a)' ) 
     &  '  poisson_1d_monogrid() solves a 1D Poisson BVP'
      write ( *, '(a)' ) '  using the Gauss-Seidel method.'

      a = 0.0D+00
      b = 1.0D+00
      ua = 0.0D+00
      ub = 0.0D+00

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) '  -u"(x) = 1, for 0 < x < 1'
      write ( *, '(a)' ) '  u(0) = u(1) = 0.'
      write ( *, '(a)' ) '  Solution is u(x) = ( -x^2 + x ) / 2'

      do k = 5, 5

        n = 2**k
        call r8vec_linspace ( n + 1, a, b, x )

        write ( *, '(a)' ) ' '
        write ( *, '(a,i4)' ) '  Mesh index K = ', k
        write ( *, '(a,i6)' ) '  Number of intervals N=2^K = ', n
        write ( *, '(a,i6)' ) '  Number of nodes = 2^K+1 =   ', 2**k + 1

        call poisson_1d_monogrid ( n, a, b, ua, ub, force1, exact1, 
     &    it_num, u )

        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) 
     &    '     I        X(I)      U(I)         U Exact(X(I))'
        write ( *, '(a)' ) ' '
        do i = 1, n + 1
          write ( *, '(2x,i4,2x,f10.4,2x,g14.6,2x,g14.6)' ) 
     &      i, x(i), u(i), exact1 ( x(i) )
        end do

        write ( *, '(a)' ) ' '

        difmax = 0.0D+00
        do i = 1, n + 1
          difmax = max ( difmax, abs ( u(i) - exact1 ( x(i) ) ) )
        end do 
        write ( *, '(a,g14.6)' ) '  Maximum error = ', difmax
        write ( *, '(a,i6)' ) 
     &    '  Number of Gauss-Seidel iterations = ', it_num

      end do

      return
      end
      subroutine test01_multi ( ) 

c*********************************************************************72
c
cc test01_multi() tests poisson_1d_multigrid() on test case 1.
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    John Burkardt
c
      implicit none

      integer n_max
      parameter ( n_max = 128 )

      double precision a
      double precision b
      double precision difmax
      external exact1
      double precision exact1
      external force1
      double precision force1
      integer i
      integer it_num
      integer k
      integer n
      double precision u(n_max+1)
      double precision ua
      double precision ub
      double precision x(n_max+1)

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'test01_multi():'
      write ( *, '(a)' ) 
     &  '  poisson_1d_multigrid() solves a 1D Poisson BVP'
      write ( *, '(a)' ) '  using the multigrid method.'

      a = 0.0D+00
      b = 1.0D+00
      ua = 0.0D+00
      ub = 0.0D+00

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) '  -u"(x) = 1, for 0 < x < 1'
      write ( *, '(a)' ) '  u(0) = u(1) = 0.'
      write ( *, '(a)' ) '  Solution is u(x) = ( -x^2 + x ) / 2'

      do k = 5, 5

        n = 2**k
        call r8vec_linspace ( n + 1, a, b, x )

        write ( *, '(a)' ) ' '
        write ( *, '(a,i4)' ) '  Mesh index K = ', k
        write ( *, '(a,i6)' ) '  Number of intervals N=2^K = ', n
        write ( *, '(a,i6)' ) '  Number of nodes = 2^K+1 =   ', n + 1

        call poisson_1d_multigrid ( n, a, b, ua, ub, force1, 
     &    exact1, it_num, u )

        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) 
     &    '     I        X(I)      U(I)         U Exact(X(I))'
        write ( *, '(a)' ) ' '
        do i = 1, n + 1
          write ( * , '(2x,i4,2x,f10.4,2x,g14.6,2x,g14.6)' ) 
     &      i, x(i), u(i), exact1 ( x(i) )
        end do

        write ( *, '(a)' ) ' '

        difmax = 0.0D+00
        do i = 1, n + 1
          difmax = max ( difmax, abs ( u(i) - exact1 ( x(i) ) ) )
        end do 
        write ( *, '(a,g14.6)' ) '  Maximum error = ', difmax
        write ( *, '(a,i6)' ) '  Number of iterations = ', it_num

      end do

      return
      end
      function exact1 ( x )

c*********************************************************************72
c                                                    
cc exact1() evaluates the exact solution for test case #1.
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    double precision X, the evaluation point.
c
c  Output:
c
c    double precision EXACT1, the value of the exact solution at X.
c
      double precision exact1
      double precision x

      exact1 = 0.5D+00 * ( - x * x + x )

      return
      end
      function force1 ( x )

c*********************************************************************72
c                                                    
cc force1() evaluates the forcing function for test case #1.
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    double precision X, the evaluation point.
c
c  Output:
c
c    double precision FORCE1, the value of the forcing function at X.
c
      double precision force1
      double precision x

      force1 = 1.0D+00

      return
      end
      subroutine test02_mono ( ) 

c*********************************************************************72
c
cc test02_mono() tests poisson_1d_monogrid() on test case 2.
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    John Burkardt
c
      implicit none

      integer n_max
      parameter ( n_max = 32 )

      double precision a
      double precision b
      double precision difmax
      external exact2
      double precision exact2
      external force2
      double precision force2
      integer i
      integer it_num
      integer k
      integer n
      double precision u(n_max+1)
      double precision ua
      double precision ub
      double precision x(n_max+1)

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'test02_mono()'
      write ( *, '(a)' ) 
     &  '  poisson_1d_monogrid() solves a 1D Poisson BVP'
      write ( *, '(a)' ) '  using the Gauss-Seidel method.'

      a = 0.0D+00
      b = 1.0D+00
      ua = 0.0D+00
      ub = 0.0D+00

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 
     &  '  -u"(x) = - x * (x+3) * exp(x), for 0 < x < 1'
      write ( *, '(a)' ) '  u(0) = u(1) = 0.'
      write ( *, '(a)' ) '  Solution is u(x) = x * (x-1) * exp(x)'

      do k = 5, 5

        n = 2**k
        call r8vec_linspace ( n + 1, a, b, x )

        write ( *, '(a)' ) ' '
        write ( *, '(a,i4)' ) '  Mesh index K = ', k
        write ( *, '(a,i6)' ) '  Number of intervals N=2^K = ', n
        write ( *, '(a,i6)' ) '  Number of nodes = 2^K+1 =   ', n + 1

        call poisson_1d_monogrid ( n, a, b, ua, ub, force2, exact2, 
     &    it_num, u )

        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) 
     &    '     I        X(I)      U(I)         U Exact(X(I))'
        write ( *, '(a)' ) ' '
        do i = 1, n + 1
          write ( *, '(2x,i4,2x,f10.4,2x,g14.6,2x,g14.6)' ) 
     &      i, x(i), u(i), exact2 ( x(i) )
        end do

        write ( *, '(a)' ) ' '

        difmax = 0.0D+00
        do i = 1, n + 1
          difmax = max ( difmax, abs ( u(i) - exact2 ( x(i) ) ) )
        end do 
        write ( *, '(a,g14.6)' ) '  Maximum error = ', difmax
        write ( *, '(a,i6)' ) 
     &    '  Number of Gauss-Seidel iterations = ', it_num

      end do

      return
      end
      subroutine test02_multi ( ) 

c*********************************************************************72
c
cc test02_multi() tests poisson_1d_multigrid() on test case 2.
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    John Burkardt
c
      implicit none

      integer n_max
      parameter ( n_max = 128 )

      double precision a
      double precision b
      double precision difmax
      external exact2
      double precision exact2
      external force2
      double precision force2
      integer i
      integer it_num
      integer k
      integer n
      double precision u(n_max+1)
      double precision ua
      double precision ub
      double precision x(n_max+1)

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'test02_multi():'
      write ( *, '(a)' ) 
     &  '  poisson_1d_multigrid() solves a 1D Poisson BVP'
      write ( *, '(a)' ) '  using the multigrid method.'

      a = 0.0D+00
      b = 1.0D+00
      ua = 0.0D+00
      ub = 0.0D+00

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 
     &  '  -u"(x) = - x * (x+3) * exp(x), for 0 < x < 1'
      write ( *, '(a)' ) '  u(0) = u(1) = 0.'
      write ( *, '(a)' ) '  Solution is u(x) = x * (x-1) * exp(x)'

      do k = 5, 5

        n = 2**k
        call r8vec_linspace ( n + 1, a, b, x )

        write ( *, '(a)' ) ' '
        write ( *, '(a,i4)' ) '  Mesh index K = ', k
        write ( *, '(a,i6)' ) '  Number of intervals N=2^K = ', n
        write ( *, '(a,i6)' ) '  Number of nodes = 2^K+1 =   ', n + 1

        call poisson_1d_multigrid ( n, a, b, ua, ub, force2, exact2, 
     &    it_num, u )

        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) 
     &    '     I        X(I)      U(I)         U Exact(X(I))'
        write ( *, '(a)' ) ' '
        do i = 1, n + 1
          write ( *, '(2x,i4,2x,f10.4,2x,g14.6,2x,g14.6)' ) 
     &      i, x(i), u(i), exact2 ( x(i) )
        end do

        write ( *, '(a)' ) ' '

        difmax = 0.0D+00
        do i = 1, n + 1
          difmax = max ( difmax, abs ( u(i) - exact2 ( x(i) ) ) )
        end do 
        write ( *, '(a,g14.6)' ) '  Maximum error = ', difmax
        write ( *, '(a,i6)' ) '  Number of iterations = ', it_num

      end do

      return
      end
      function exact2 ( x )

c*********************************************************************72
c                                                    
cc exact2() evaluates the exact solution for test case #2.
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    double precision X, the evaluation point.
c
c  Output:
c
c    double precision EXACT2, the value of the exact solution at X.
c
      double precision exact2
      double precision x

      exact2 = x * ( x - 1.0D+00 ) * exp ( x )

      return
      end
      function force2 ( x )

c*********************************************************************72
c                                                    
cc force2() evaluates the forcing function for test case #2.
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    double precision X, the evaluation point.
c
c  Output:
c
c    double precision FORCE2, the value of the forcing function at X.
c
      double precision force2
      double precision x

      force2 = - x * ( x + 3.0D+00 ) * exp ( x )

      return
      end
      subroutine timestamp ( )

c*********************************************************************72
c
cc timestamp() prints out the current YMDHMS date as a timestamp.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    12 January 2007
c
c  Author:
c
c    John Burkardt
c
c  Parameters:
c
c    None
c
      implicit none

      character * ( 8 ) ampm
      integer d
      character * ( 8 ) date
      integer h
      integer m
      integer mm
      character * ( 9 ) month(12)
      integer n
      integer s
      character * ( 10 ) time
      integer y

      save month

      data month /
     &  'January  ', 'February ', 'March    ', 'April    ', 
     &  'May      ', 'June     ', 'July     ', 'August   ', 
     &  'September', 'October  ', 'November ', 'December ' /

      call date_and_time ( date, time )

      read ( date, '(i4,i2,i2)' ) y, m, d
      read ( time, '(i2,i2,i2,1x,i3)' ) h, n, s, mm

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

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

      return
      end