1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
|
!!
!! Liebmann's method to compute heat transfer across a 2-D surface
!! J. Overbey 8/27/08, updated for OpenACC 5/3/12
!!
!! Use liebmann-viz.sh in the main project directory to visualize the resulting
!! table using gnuplot
!!
!! This version is based on the code in src-liebmann-3-loops, but
!! (1) the "iterate" subroutine is inlined,
!! (2) the number of iterations is always a multiple of 2 (delta is
!! only computed on even-numbered iterations), and
!! (3) with an OpenACC-compliant compiler (e.g., Cray Fortran 8.0), the
!! stencil computations will be moved to an accelerator device/GPU
!!
program liebmann_example
implicit none
integer, parameter :: SIZE = 4096
integer, parameter :: INTERIOR_SIZE = SIZE - 2
integer, parameter :: OUTPUT_SIZE = 128
real, parameter :: BOUNDARY_VALUE = 5.0
real, parameter :: EPSILON = 0.001
call main()
contains
subroutine main()
real :: surface(SIZE, SIZE)
integer :: i, j
integer :: start_time, end_time, count_rate
call system_clock(start_time, count_rate)
call liebmann(surface)
call system_clock(end_time)
do i = 1, SIZE, SIZE/OUTPUT_SIZE
do j = 1, SIZE-1, SIZE/OUTPUT_SIZE
write (*,'(F4.2" ")',advance="no") surface(i, j)
end do
write (*,'(F4.2" ")') surface(i, SIZE)
end do
print *, "Elapsed Time (seconds):", real(end_time - start_time) / real(count_rate)
end subroutine
subroutine liebmann(surface)
real, dimension(SIZE, SIZE), intent(out) :: surface
real, dimension(SIZE, SIZE) :: prev, next
real :: delta
integer :: i, j
logical :: done
call init_with_boundaries(prev)
call init_with_boundaries(next)
!$acc data copyin(prev,next) copyout(prev)
done = .false.
do while (.not. done)
!$acc parallel loop
do j = 2, SIZE-1
do i = 2, SIZE-1
next(i,j) = &
(prev(i-1, j) + &
prev(i+1, j) + &
prev(i, j-1) + &
prev(i, j+1)) / 4.0
end do
end do
!$acc end parallel loop
delta = 0.0
!$acc parallel loop reduction(max:delta)
do j = 2, SIZE-1
do i = 2, SIZE-1
prev(i,j) = &
(next(i-1, j) + &
next(i+1, j) + &
next(i, j-1) + &
next(i, j+1)) / 4.0
delta = max(delta, abs(prev(i,j)-next(i,j)))
end do
end do
!$acc end parallel loop
if (delta < EPSILON) then
done = .true.
end if
end do
!$acc end data
surface = prev
end subroutine
subroutine init_with_boundaries(surface)
real, dimension(SIZE, SIZE), intent(out) :: surface
surface = 0.0
surface(1, :) = BOUNDARY_VALUE
surface(SIZE, :) = BOUNDARY_VALUE
surface(:, 1) = BOUNDARY_VALUE
surface(:, SIZE) = BOUNDARY_VALUE
end subroutine
end program
|
