program main !*****************************************************************************80 ! !! MAIN is the main program for JACOBI. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! use omp_lib integer n double precision, allocatable, dimension ( :, : ) :: a double precision, allocatable, dimension ( : ) :: b integer seed double precision, allocatable, dimension ( : ) :: x n = 500 seed = 123456789 write ( *, * ) ' ' write ( *, * ) 'JACOBI' write ( *, * ) ' FORTRAN90 version' write ( *, '(a,i8)' ) ' Number of threads = ', omp_get_max_threads ( ) write ( *, '(a,i8)' ) ' Problem size n = ', n allocate ( a(1:n,1:n) ) allocate ( b(1:n) ) allocate ( x(1:n) ) call r8vec_uniform_01 ( n * n, seed, a ) do i = 1, n a(i,i) = 1.0 + 2.0 * sum ( abs ( a(i,1:n) ) ) end do do i = 1, n x(i) = i end do b(1:n) = matmul ( a(1:n,1:n), x(1:n) ) call r8vec_uniform_01 ( n, seed, x ) call jacobi ( n, a, b, x ) deallocate ( a ) deallocate ( b ) deallocate ( x ) stop end subroutine jacobi ( n, a, b, x ) !*****************************************************************************80 ! !! JACOBI carries out the Jacobi iteration. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! use omp_lib integer n double precision a(n,n) double precision b(n) double precision diag(n) double precision diff double precision diff_tol integer i integer it integer, parameter :: it_max = 100 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 ) * epsilon ( diff_tol ) 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 <= diff_tol ) then exit end if end do 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 ) !*****************************************************************************80 ! !! R8VEC_UNIFORM_01 returns a unit pseudorandom vector. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Parameters: ! ! Input, integer ( kind = 4 ) N, the number of entries in the vector. ! ! Input/output, integer ( kind = 4 ) SEED, the "seed" value, which ! should NOT be 0. On output, SEED has been updated. ! ! Output, real ( kind = 8 ) R(N), the vector of pseudorandom values. ! implicit none integer ( kind = 4 ) n integer ( kind = 4 ) i integer ( kind = 4 ), parameter :: i4_huge = 2147483647 integer ( kind = 4 ) k integer ( kind = 4 ) seed real ( kind = 8 ) r(n) if ( seed == 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'R8VEC_UNIFORM_01 - Fatal error!' 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 < 0 ) then seed = seed + i4_huge end if r(i) = real ( seed, kind = 8 ) * 4.656612875D-10 end do return end