program main

c*********************************************************************72
c
cc jacobi_test() tests jacobi().
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    13 January 2013
c
c  Author:
c
c    John Burkardt
c
      implicit none

      call timestamp ( )
      write ( *, '(a)' ) ''
      write ( *, '(a)' ) 'jacobi_test():'
      write ( *, '(a)' ) '  FORTRAN77 version.'
      write ( *, '(a)' ) '  Test jacobi().'

      call test01 ( )
c
c  Terminate.
c
      write ( *, '(a)' ) ''
      write ( *, '(a)' ) 'JACOBI_test():'
      write ( *, '(a)' ) '  Normal end of execution.'
      write ( *, '(a)' ) ''
      call timestamp ( )

      return
      end
      subroutine test01 ( )

c*********************************************************************72
c
c  Purpose:
c
c    jacobi_test01() tests jacobi1().
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    13 January 2013
c
c  Author:
c
c    John Burkardt
c
      implicit none

      integer n
      parameter ( n = 33 )

      double precision a(n,n)
      double precision b(n)
      double precision d_max
      double precision e_max
      integer i
      integer it
      integer it_max
      double precision r(n)
      double precision r_max
      double precision t
      double precision tol
      double precision x(n)
      double precision x_exact(n)
      double precision x_new(n)

      write ( *, '(a)' ) ''
      write ( *, '(a)' ) 'jacobi_test01():'
      write ( *, '(a)' ) '  Try the Jacobi iteration on the second'
      write ( *, '(a)' ) '  difference matrix'

      it_max = 2000
      tol = 1.0D-05

      write ( *, '(a)' ) ' '
      write ( *, '(a,i4)' ) '  Matrix order N = ', n
      write ( *, '(a,i6)' ) '  Maximum number of iterations = ', it_max
c
c  Set the matrix A.
c
      call dif2 ( n, n, a )
c
c  Determine the right hand side vector B.
c
      do i = 1, n
        t = dble ( i - 1 ) / dble ( n - 1 )
        x_exact(i) = exp ( t ) * ( t - 1 ) * t
      end do

      b = matmul ( a(1:n,1:n), x_exact(1:n) )
c
c  Set the initial estimate for the solution.
c
      it = 0

      x(1:n) = 0.0D+00

      r = matmul ( a(1:n,1:n), x(1:n) )
      r(1:n) = b(1:n) - r(1:n)
      r_max = maxval ( r );

      e_max = maxval ( abs ( x - x_exact ) )

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 
     &  '       I    Resid           X-Error         X-Change'
      write ( *, '(a)' ) ' '

      write ( *, '(2x,i6,2x,g14.6,2x,g14.6)' ) it, r_max, e_max
c
c  Carry out the iteration.
c
      do it = 1, it_max

        call jacobi1 ( n, a, b, x, x_new )

        r = matmul ( a(1:n,1:n), x(1:n) )
        r(1:n) = b(1:n) - r(1:n)
        r_max = maxval ( abs ( r ) )
c
c  Compute the average point motion.
c
        d_max = maxval ( abs ( x - x_new ) )
c
c  Compute the average point motion.
c
        e_max = maxval ( abs ( x - x_exact ) )
c
c  Update the solution
c
        x(1:n) = x_new(1:n)

        write ( *, '(2x,i6,2x,g14.6,2x,g14.6,2x,g14.6)' ) 
     &    it, r_max, e_max, d_max

        if ( r_max <= tol ) then
          write ( *, '(a)' ) ' '
          write ( *, '(a,i6)' ) '  Convergence on step ', it
          exit
        end if

      end do

      call r8vec_print ( n, x, "  Estimated solution:" )

      return
      end