From cfde90ff81cd5828fc7f01ee41b38dd225744e16 Mon Sep 17 00:00:00 2001
From: Damian Rouson <rouson@lbl.gov>
Date: Thu, 27 Apr 2023 02:24:11 -0700
Subject: [PATCH 1/2] feat(example): add procedural alternative

---
 example/procedural.f90 | 137 +++++++++++++++++++++++++++++++++++++++++
 1 file changed, 137 insertions(+)
 create mode 100644 example/procedural.f90

diff --git a/example/procedural.f90 b/example/procedural.f90
new file mode 100644
index 0000000..0e0c677
--- /dev/null
+++ b/example/procedural.f90
@@ -0,0 +1,137 @@
+program heat
+  implicit none
+
+  real dx, dy, dt
+  real, parameter :: alpha = 1.
+  real, allocatable :: T(:,:), halo_x(:,:)[:]
+  integer, parameter :: west=1, east=2, nx=11, ny=nx, steps = 1000
+  integer me, my_nx, num_subdomains, my_internal_west, my_internal_east, i
+
+  me = this_image()
+  num_subdomains = num_images()
+
+  call define(side=1., boundary_val=1., internal_val=2., n=nx) ! 2D step function
+  sync all
+
+  dt = dx*dy/(4*alpha)
+
+  do i = 1, steps
+    call step(T, dt)
+  end do
+
+  critical
+    do i = 1, size(T,2)
+      print *,"image ",me,": ", T(:,i)
+    end do
+  end critical
+
+contains
+
+  pure subroutine assert(assertion, description)
+    logical, intent(in) :: assertion
+    character(len=*), intent(in) :: description
+    if (.not. assertion) error stop description
+  end subroutine
+
+  subroutine exchange_halo(s)
+    real, intent(in) :: s(:,:)
+    if (me>1) halo_x(east,:)[me-1] = s(1,:)
+    if (me<num_subdomains) halo_x(west,:)[me+1] = s(my_nx,:)
+  end subroutine
+
+  subroutine define(side, boundary_val, internal_val, n)
+    real, intent(in) :: side, boundary_val, internal_val
+    integer, intent(in) :: n !! number of grid points in each coordinate direction
+    integer, parameter :: nx_boundaries = 2 
+
+    call assert(num_subdomains <= n-nx_boundaries, "define: num_subdomains <= n-nx_boundaries")
+
+    dx = side/(nx-1)
+    dy = dx 
+    my_nx = nx/num_subdomains + merge(1, 0, me <= mod(nx, num_subdomains))
+    my_internal_west = merge(2, 1, me==1)
+    my_internal_east = merge(my_nx-1, my_nx, me==num_subdomains)
+
+    allocate(T(my_nx, ny))
+
+    T(my_internal_west:my_internal_east, 1) = boundary_val ! bottom subdomain boundary
+    T(my_internal_west:my_internal_east, ny) = boundary_val ! top subdomain boundary
+    T(my_internal_west:my_internal_east, 2:ny-1) = internal_val ! internal points
+    T(1, 2:ny-1) = merge(boundary_val, internal_val, me==1) ! west subdomain boundary
+    T(my_nx, 2:ny-1) = merge(boundary_val, internal_val, me==num_subdomains) ! east subdomain boundary
+
+    if (allocated(halo_x)) deallocate(halo_x)
+    allocate(halo_x(west:east, ny)[*])
+    call exchange_halo(T)
+  end subroutine
+
+  subroutine internal_points(alpha_dt_Laplacian_T)
+    real, intent(inout) :: alpha_dt_Laplacian_T(:,:)
+    integer i, j
+
+    do concurrent(i=my_internal_west+1:my_internal_east-1, j=2:ny-1)
+      alpha_dt_Laplacian_T(i,j) = alpha*dt*( &
+        (T(i-1,j) - 2*T(i,j) + T(i+1,j))/dx**2 + &
+        (T(i,j-1) - 2*T(i,j) + T(i,j+1))/dy**2 &
+      )
+    end do
+
+  end subroutine
+
+  subroutine edge_points(ds)
+    real, intent(inout) :: ds(:,:)
+    real, allocatable :: halo_west(:), halo_east(:)
+    integer i, j
+
+    halo_west = merge(halo_x(west,:), T(1,:), me/=1)
+    halo_east = merge(halo_x(east,:), T(my_nx,:), me/=num_subdomains)
+
+    i = my_internal_west
+    do concurrent(j=2:ny-1)
+      ds(i,j) = alpha*dt*( &
+        (halo_west(j)   - 2*T(i,j) + T(i+1,j))/dx**2 + &
+        (T(i,j-1)  - 2*T(i,j) + T(i,j+1))/dy**2 &
+      )
+    end do
+
+    i = my_internal_east
+    do concurrent(j=2:ny-1)
+      ds(i,j) = alpha*dt*( &
+        (T(i-1,j)  - 2*T(i,j) + halo_east(j))/dx**2 + &
+        (T(i,j-1) - 2*T(i,j) + T(i,j+1))/dy**2 &
+      )
+    end do
+  end subroutine
+
+  subroutine apply_boundary_condition(ds)
+    real, intent(inout) :: ds(:,:)
+    ds(:, 1) = 0.
+    ds(:,ny) = 0.
+    if (me==1) ds(1,:) = 0.
+    if (me==num_subdomains) ds(my_nx,:) = 0.
+  end subroutine
+
+  subroutine step(T, dt)
+    real, intent(inout), allocatable :: T(:,:)
+    real, intent(in) :: dt
+    real, allocatable :: increment(:,:)
+
+    call assert(allocated(T), "laplacian: allocated(T)")
+    call assert(allocated(halo_x), "laplacian: allocated(halo_x)")
+    call assert(my_internal_west+1<=my_nx, "laplacian: westernmost size")
+    call assert(my_internal_east-1>0, "laplacian: subdomain size")
+
+    allocate(increment(my_nx,ny))
+ 
+    call internal_points(increment)
+    call edge_points(increment)
+    call apply_boundary_condition(increment)
+
+    sync all
+    T = T + increment
+    sync all
+    call exchange_halo(T)
+
+  end subroutine
+
+end program

From 480a806dd8e869e4dea18ece96efef720f659ae4 Mon Sep 17 00:00:00 2001
From: Damian Rouson <rouson@lbl.gov>
Date: Thu, 27 Apr 2023 03:14:00 -0700
Subject: [PATCH 2/2] feat(example/procedural): simplify

---
 example/procedural.f90 | 120 ++++++++++++++++-------------------------
 1 file changed, 46 insertions(+), 74 deletions(-)

diff --git a/example/procedural.f90 b/example/procedural.f90
index 0e0c677..d31f937 100644
--- a/example/procedural.f90
+++ b/example/procedural.f90
@@ -1,29 +1,34 @@
-program heat
+program main
   implicit none
 
-  real dx, dy, dt
-  real, parameter :: alpha = 1.
-  real, allocatable :: T(:,:), halo_x(:,:)[:]
-  integer, parameter :: west=1, east=2, nx=11, ny=nx, steps = 1000
+  real dx, dy, dt, t_start, t_end, T_min, T_max
+  real, parameter :: alpha = 1., T_boundary = 1., T_initial = 2.
+  real, allocatable :: T(:,:), halo_x(:,:)[:], laplacian(:,:)
+  integer, parameter :: west=1, east=2, nx=51, ny=nx, steps = 3000
   integer me, my_nx, num_subdomains, my_internal_west, my_internal_east, i
 
   me = this_image()
   num_subdomains = num_images()
-
-  call define(side=1., boundary_val=1., internal_val=2., n=nx) ! 2D step function
+  call initialize(T, side=1., boundary_val=T_boundary, internal_val=T_initial, n=nx) ! 2D step function
+  allocate(halo_x(west:east, ny)[*]) ! implicit synchronization
+  call exchange_halo(T)
   sync all
-
   dt = dx*dy/(4*alpha)
 
+  call cpu_time(t_start)
   do i = 1, steps
     call step(T, dt)
+    call exchange_halo(T)
+    sync all
   end do
+  call cpu_time(t_end)
 
-  critical
-    do i = 1, size(T,2)
-      print *,"image ",me,": ", T(:,i)
-    end do
-  end critical
+  T_min = minval(T)
+  T_max = maxval(T)
+  call co_min(T_min, result_image=1)
+  call co_min(T_max, result_image=1)
+  if (me==1) print *,"T_boundary, T_ininitial, T_min, T_max: ", T_boundary, T_initial, T_min, T_max
+  print *,"cpu_time: ", t_end - t_start
 
 contains
 
@@ -39,7 +44,8 @@ subroutine exchange_halo(s)
     if (me<num_subdomains) halo_x(west,:)[me+1] = s(my_nx,:)
   end subroutine
 
-  subroutine define(side, boundary_val, internal_val, n)
+  subroutine initialize(T, side, boundary_val, internal_val, n)
+    real, intent(out), allocatable :: T(:,:)
     real, intent(in) :: side, boundary_val, internal_val
     integer, intent(in) :: n !! number of grid points in each coordinate direction
     integer, parameter :: nx_boundaries = 2 
@@ -59,79 +65,45 @@ subroutine define(side, boundary_val, internal_val, n)
     T(my_internal_west:my_internal_east, 2:ny-1) = internal_val ! internal points
     T(1, 2:ny-1) = merge(boundary_val, internal_val, me==1) ! west subdomain boundary
     T(my_nx, 2:ny-1) = merge(boundary_val, internal_val, me==num_subdomains) ! east subdomain boundary
-
-    if (allocated(halo_x)) deallocate(halo_x)
-    allocate(halo_x(west:east, ny)[*])
-    call exchange_halo(T)
-  end subroutine
-
-  subroutine internal_points(alpha_dt_Laplacian_T)
-    real, intent(inout) :: alpha_dt_Laplacian_T(:,:)
-    integer i, j
-
-    do concurrent(i=my_internal_west+1:my_internal_east-1, j=2:ny-1)
-      alpha_dt_Laplacian_T(i,j) = alpha*dt*( &
-        (T(i-1,j) - 2*T(i,j) + T(i+1,j))/dx**2 + &
-        (T(i,j-1) - 2*T(i,j) + T(i,j+1))/dy**2 &
-      )
-    end do
-
-  end subroutine
-
-  subroutine edge_points(ds)
-    real, intent(inout) :: ds(:,:)
-    real, allocatable :: halo_west(:), halo_east(:)
-    integer i, j
-
-    halo_west = merge(halo_x(west,:), T(1,:), me/=1)
-    halo_east = merge(halo_x(east,:), T(my_nx,:), me/=num_subdomains)
-
-    i = my_internal_west
-    do concurrent(j=2:ny-1)
-      ds(i,j) = alpha*dt*( &
-        (halo_west(j)   - 2*T(i,j) + T(i+1,j))/dx**2 + &
-        (T(i,j-1)  - 2*T(i,j) + T(i,j+1))/dy**2 &
-      )
-    end do
-
-    i = my_internal_east
-    do concurrent(j=2:ny-1)
-      ds(i,j) = alpha*dt*( &
-        (T(i-1,j)  - 2*T(i,j) + halo_east(j))/dx**2 + &
-        (T(i,j-1) - 2*T(i,j) + T(i,j+1))/dy**2 &
-      )
-    end do
-  end subroutine
-
-  subroutine apply_boundary_condition(ds)
-    real, intent(inout) :: ds(:,:)
-    ds(:, 1) = 0.
-    ds(:,ny) = 0.
-    if (me==1) ds(1,:) = 0.
-    if (me==num_subdomains) ds(my_nx,:) = 0.
   end subroutine
 
   subroutine step(T, dt)
     real, intent(inout), allocatable :: T(:,:)
     real, intent(in) :: dt
-    real, allocatable :: increment(:,:)
+    real, allocatable :: halo_west(:), halo_east(:)
+    integer i, j
 
     call assert(allocated(T), "laplacian: allocated(T)")
     call assert(allocated(halo_x), "laplacian: allocated(halo_x)")
     call assert(my_internal_west+1<=my_nx, "laplacian: westernmost size")
     call assert(my_internal_east-1>0, "laplacian: subdomain size")
+    if (.not. allocated(laplacian)) allocate(laplacian(my_nx,ny))
 
-    allocate(increment(my_nx,ny))
- 
-    call internal_points(increment)
-    call edge_points(increment)
-    call apply_boundary_condition(increment)
+    internal_points: &
+    do concurrent(i=my_internal_west+1:my_internal_east-1, j=2:ny-1)
+      laplacian(i,j) = ((T(i-1,j) - 2*T(i,j) + T(i+1,j))/dx**2 + &
+                        (T(i,j-1) - 2*T(i,j) + T(i,j+1))/dy**2 )
+    end do internal_points
 
-    sync all
-    T = T + increment
-    sync all
-    call exchange_halo(T)
+    i = my_internal_west
+    halo_west = merge(halo_x(west,:), T(1,:), me/=1)
+
+    west_boundary_adjacent_points: &
+    do concurrent(j=2:ny-1)
+      laplacian(i,j) = ( (halo_west(j) - 2*T(i,j) + T(i+1,j))/dx**2 + &
+                         (T(i,j-1)     - 2*T(i,j) + T(i,j+1))/dy**2 )
+    end do west_boundary_adjacent_points
+    
+    i = my_internal_east
+    halo_east = merge(halo_x(east,:), T(my_nx,:), me/=num_subdomains)
+
+    east_boundary_adjacent_points: &
+    do concurrent(j=2:ny-1)
+      laplacian(i,j) = ( (T(i-1,j) - 2*T(i,j) + halo_east(j))/dx**2 + &
+                         (T(i,j-1) - 2*T(i,j) + T(i,j+1)    )/dy**2 )
+    end do east_boundary_adjacent_points
 
+    T = T + alpha*dt*laplacian
   end subroutine
 
 end program