program main

c*********************************************************************72
c
cc MAIN is the main program for JACOBI.
c
c  Licensing:
c
c    This code is distributed under the GNU LGPL license. 
c
      include 'omp_lib.h'

      integer n
      parameter ( n = 500 )

      double precision a(n,n)
      double precision b(n)
      integer i
      integer j
      double precision row_sum
      integer seed
      double precision x(n)
      double precision x_old(n)

      write ( *, * ) ' '
      write ( *, * ) 'JACOBI'
      write ( *, * ) '  FORTRAN77 version'
      write ( *, '(a,i8)' ) '  Number of threads = ', 
     &  omp_get_max_threads ( )
      write ( *, '(a,i8)' ) '  Problem size N = ', n

      seed = 123456789

      call r8vec_uniform_01 ( n * n, seed, a )

      do i = 1, n
        row_sum = 0.0
        do j = 1, n
          row_sum = row_sum + abs ( a(i,j) )
        end do
        a(i,i) = 1.0 + 2.0 * row_sum
      end do

      do i = 1, n
        x(i) = dble ( i )
      end do

      do i = 1, n
        b(i) = 0.0
        do j = 1, n
          b(i) = b(i) + a(i,j) * x(j)
        end do
      end do

      call r8vec_uniform_01 ( n, seed, x )

      call jacobi ( n, a, b, x, x_old )

      stop
      end
      subroutine jacobi ( n, a, b, x, x_old )

c*********************************************************************72
c
cc JACOBI carries out the Jacobi iteration.
c
c  Licensing:
c
c    This code is distributed under the GNU LGPL license. 
c
      include 'omp_lib.h'

      integer n

      double precision a(n,n)
      double precision axi
      double precision b(n)
      double precision diff
      double precision diff_tol
      integer i
      integer it
      integer it_max
      parameter ( it_max = 100 )
      integer j
      double precision r8_epsilon
      parameter ( r8_epsilon = 1.0D-14 )
      double precision wtime
      double precision wtime1
      double precision wtime2
      double precision x(n)
      double precision x_old(n)

      diff_tol = 0.0
      do i = 1, n
        diff_tol = diff_tol + abs( b(i) )
      end do
      diff_tol = ( diff_tol + 1.0 ) * r8_epsilon

      wtime1 = omp_get_wtime ( )

      do it = 1, it_max

        do i = 1, n
          x_old(i) = x(i)
        end do

        do i = 1, n
          axi = 0.0
          do j = 1, n
            axi = axi + a(i,j) * x_old(j)
          end do
          x(i) = x_old(i) + ( b(i) - axi ) / a(i,i)
        end do

        diff = 0.0
        do i = 1, n
          diff = diff + abs ( x(i) - x_old(i) )
        end do

        write ( *, * ) it, diff
        if ( diff .le. diff_tol ) then
          go to 10
        end if

      end do

10    continue

      wtime2 = omp_get_wtime ( )
      wtime = wtime2 - wtime1

      write ( *, * ) ' '
      write ( *, * ) 'DIFF     = ', diff
      write ( *, * ) 'DIFF_TOL = ', diff_tol
      write ( *, * ) 'Time     = ', wtime

      write ( *, * ) ' '
      write ( *, * ) '  First 10 elements of estimated solution'
      write ( *, * ) ' '

      do i = 1, 10
        write ( *, * ) i, x(i)
      end do

      return
      end
      subroutine r8vec_uniform_01 ( n, seed, r )

c*********************************************************************72
c
cc R8VEC_UNIFORM_01 returns a unit pseudorandom vector.
c
c  Licensing:
c
c    This code is distributed under the GNU LGPL license. 
c
c  Parameters:
c
c    Input, integer N, the number of entries in the vector.
c
c    Input/output, integer SEED, the "seed" value, which 
c    should NOT be 0.  On output, SEED has been updated.
c
c    Output, double precision R(N), the vector of pseudorandom values.
c
      implicit none

      integer n

      integer i
      integer i4_huge
      parameter ( i4_huge = 2147483647 )
      integer k
      integer seed
      double precision r(n)

      if ( seed .eq. 0 ) then
        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) 'R8VEC_UNIFORM_01 - Fatal errorc'
        write ( *, '(a)' ) '  Input value of SEED = 0.'
        stop
      end if

      do i = 1, n

        k = seed / 127773

        seed = 16807 * ( seed - k * 127773 ) - k * 2836

        if ( seed .lt. 0 ) then
          seed = seed + i4_huge
        end if

        r(i) = dble ( seed ) * 4.656612875D-10

      end do

      return
      end