WIP: nagfor/gfortran internal compiler error (ICE)

gfortran 13.2.0
---------------
The following command reproduces the ICE:

gfortran test

NAG 7.1 (Build 7143)
--------------------
The following command reproduces the ICE:

fpm test --compiler nagfor --flag "-fpp -f2018"

src/matcha/do_concurrent_s.f90

@@ -30,7 +30,7 @@ pure module subroutine do_concurrent_my_velocities(nsteps, dir, sampled_speeds,
 
     call assert(all([size(my_velocities,1),size(sampled_speeds,2)] == shape(sampled_speeds)), &
       "do_concurrent_my_velocities: argument size match")
-    call assert(all(shape(my_velocities,1)==shape(dir)), "do_concurrent_my_velocities: argument shape match")
+    call assert(all(size(my_velocities,1)==shape(dir)), "do_concurrent_my_velocities: argument shape match")
 
     do concurrent(step=1:nsteps)
       my_velocities(:,step,1) = sampled_speeds(:,step)*dir(:,step,1)
@@ -82,16 +82,14 @@ module subroutine do_concurrent_speeds(history, speeds) bind(C)
          x(i,:,:) = positions
       end do
 
-      associate(t => history%time)
-        do concurrent(i = 1:npositions-1, j = 1:ncells)
-          associate( &
-            u => (x(i+1,j,:) - x(i,j,:))/(t(i+1) - t(i)), &
-            ij => i + (j-1)*(npositions-1) &
-           )   
-            speeds(ij) = sqrt(sum([(u(k)**2, k=1,nspacedims)]))
-          end associate
-        end do
-      end associate
+      do concurrent(i = 1:npositions-1, j = 1:ncells)
+        associate( &
+          u => (x(i+1,j,:) - x(i,j,:))/(history(i+1)%time - history(i)%time), &
+          ij => i + (j-1)*(npositions-1) &
+         )   
+          speeds(ij) = sqrt(sum([(u(k)**2, k=1,nspacedims)]))
+        end associate
+      end do
     end associate
 
   end subroutine

src/matcha/output_s.f90

@@ -24,23 +24,26 @@
 
     integer, parameter :: speed=1, freq=2 ! subscripts for speeds and frequencies
 
-    associate(npositions => size(history), ncells => history(1)%positions_shape(1))
-      allocate(speeds(ncells*(npositions-1)))
+    associate(npositions => size(self%history_))
+      allocate(speeds(self%my_num_cells()*(npositions-1)))
     end associate
     call do_concurrent_speeds(t_cell_collection_bind_C_t(self%history_), speeds)
 
-    associate(emp_distribution => self%input_%sample_distribution())
+    block
+      real(c_double), allocatable :: emp_distribution(:,:)
+
+      emp_distribution = self%input_%sample_distribution()
       associate(nintervals => size(emp_distribution(:,1)), dvel_half => (emp_distribution(2,speed)-emp_distribution(1,speed))/2.d0)
         vel = [emp_distribution(1,speed) - dvel_half, [(emp_distribution(i,speed) + dvel_half, i=1,nintervals)]]
         if (allocated(k)) deallocate(k)
-        allocate(k(nspeeds))
+        allocate(k(size(speeds)))
         call do_concurrent_k(speeds, vel, k)
         if(allocated(output_distribution)) deallocate(output_distribution)
         allocate(output_distribution(nintervals,2))
-        call do_concurrent_output_distribution(nintervals, speed, freq, emp_distribution, k, output_distribution)
+        call do_concurrent_output_distribution(speed, freq, emp_distribution, k, output_distribution)
         output_distribution(:,freq) = output_distribution(:,freq)/sum(output_distribution(:,freq))
       end associate
-    end associate
+    end block
 
   end procedure
 

src/matcha/subdomain_m.f90

@@ -1,42 +1,36 @@
 module subdomain_m
-  use assert_m, only : assert
   implicit none
 
   private
   public :: subdomain_t
-  public :: operator(.laplacian.)
-  public :: step
 
   type subdomain_t 
     private
     real, allocatable :: s_(:,:,:)
   contains
     procedure, pass(self) :: define
-    procedure, pass(rhs) :: multiply
-    generic :: operator(*) => multiply
-    generic :: operator(+) => add
-    generic :: assignment(=) => assign_
     procedure dx
     procedure dy
     procedure dz
     procedure values
+    generic :: operator(*) => multiply
+    generic :: operator(+) => add
+    generic :: operator(.laplacian.) => laplacian
+    generic :: assignment(=) => assign_
+    procedure, private, pass(rhs) :: multiply
+    procedure, private :: laplacian
     procedure, private :: add
     procedure, private :: assign_
   end type
 
-  interface operator(.laplacian.)
-
-    module procedure laplacian
-    !pure module function laplacian(rhs) result(laplacian_rhs)
-    !  implicit none
-    !  type(subdomain_t), intent(in) :: rhs[*]
-    !  type(subdomain_t) laplacian_rhs
-    !end function
-
-  end interface
-
   interface
 
+    pure module function laplacian(rhs) result(laplacian_rhs)
+      implicit none
+      class(subdomain_t), intent(in) :: rhs[*]
+      type(subdomain_t) laplacian_rhs
+    end function
+
     module subroutine define(side, boundary_val, internal_val, n, self)
       implicit none
       real, intent(in) :: side, boundary_val, internal_val
@@ -96,62 +90,4 @@ module subroutine assign_(lhs, rhs)
 
   end interface
 
-  real dx_, dy_, dz_
-  integer my_nx, nx, ny, nz, me, num_subdomains, my_internal_west, my_internal_east
-
-contains
-
-  pure module function laplacian(rhs) result(laplacian_rhs)
-    type(subdomain_t), intent(in) :: rhs[*]
-    type(subdomain_t) laplacian_rhs
-
-    integer i, j, k
-    real, allocatable :: halo_west(:,:), halo_east(:,:)
-
-    call assert(allocated(rhs%s_), "subdomain_t%laplacian: allocated(rhs%s_)")
-
-    allocate(laplacian_rhs%s_, mold=rhs%s_)
-
-    if (me==1) then
-      halo_west = rhs%s_(1,:,:)
-    else
-      halo_west = rhs[me-1]%s_(ubound(rhs[me-1]%s_,1),:,:)
-    end if
-    i = my_internal_west
-    call assert(i+1<=my_nx,"laplacian: westernmost subdomain too small")
-    do concurrent(j=2:ny-1, k=2:nz-1)
-      laplacian_rhs%s_(i,j,k) = ( halo_west(j,k  ) - 2*rhs%s_(i,j,k) + rhs%s_(i+1,j  ,k  ))/dx_**2 + &
-                                (rhs%s_(i,j-1,k  ) - 2*rhs%s_(i,j,k) + rhs%s_(i  ,j+1,k  ))/dy_**2 + &
-                                (rhs%s_(i,j  ,k-1) - 2*rhs%s_(i,j,k) + rhs%s_(i  ,j  ,k+1))/dz_**2
-    end do
-
-    do concurrent(i=my_internal_west+1:my_internal_east-1, j=2:ny-1, k=2:nz-1)
-      laplacian_rhs%s_(i,j,k) = (rhs%s_(i-1,j  ,k  ) - 2*rhs%s_(i,j,k) + rhs%s_(i+1,j  ,k  ))/dx_**2 + &
-                                (rhs%s_(i  ,j-1,k  ) - 2*rhs%s_(i,j,k) + rhs%s_(i  ,j+1,k  ))/dy_**2 + &
-                                (rhs%s_(i  ,j  ,k-1) - 2*rhs%s_(i,j,k) + rhs%s_(i  ,j  ,k+1))/dz_**2
-    end do
-
-    if (me==1) then
-      halo_east = rhs%s_(1,:,:)
-    else
-      halo_east = rhs[me+1]%s_(lbound(rhs[me+1]%s_,1),:,:)
-    end if
-    i = my_internal_east
-    call assert(i-1>0,"laplacian: easternmost subdomain too small")
-    do concurrent(j=2:ny-1, k=2:nz-1)
-      laplacian_rhs%s_(i,j,k) = (rhs%s_(i-1,j  ,k  ) - 2*rhs%s_(i,j,k) +  halo_east(j  ,k  ))/dx_**2 + &
-                                (rhs%s_(i  ,j-1,k  ) - 2*rhs%s_(i,j,k) + rhs%s_(i  ,j+1,k  ))/dy_**2 + &
-                                (rhs%s_(i  ,j  ,k-1) - 2*rhs%s_(i,j,k) + rhs%s_(i  ,j  ,k+1))/dz_**2
-    end do
-
-    laplacian_rhs%s_(:, 1,:) = 0.
-    laplacian_rhs%s_(:,ny,:) = 0.
-    laplacian_rhs%s_(:,:, 1) = 0.
-    laplacian_rhs%s_(:,:,nz) = 0.
-    if (me==1) laplacian_rhs%s_(1,:,:) = 0.
-    if (me==num_subdomains) laplacian_rhs%s_(my_nx,:,:) = 0.
-
-  end function
-
-
 end module

src/matcha/subdomain_s.f90

@@ -1,10 +1,11 @@
 submodule(subdomain_m) subdomain_s
+  use assert_m, only : assert, intrinsic_array_t
   use sourcery_m, only : data_partition_t
-  use intrinsic_array_m, only : intrinsic_array_t
   implicit none
 
   type(data_partition_t) data_partition
-
+  real dx_, dy_, dz_
+  integer my_nx, nx, ny, nz, me, num_subdomains, my_internal_west, my_internal_east
   real, allocatable :: increment(:,:,:)
 
 contains
@@ -144,4 +145,56 @@ subroutine apply_boundary_condition(ds)
 
   end procedure
 
-end submodule subdomain_s
+  pure module function laplacian(rhs) result(laplacian_rhs)
+    class(subdomain_t), intent(in) :: rhs[*]
+    type(subdomain_t) laplacian_rhs
+
+    integer i, j, k
+    real, allocatable :: halo_west(:,:), halo_east(:,:)
+
+    call assert(allocated(rhs%s_), "subdomain_t%laplacian: allocated(rhs%s_)")
+
+    allocate(laplacian_rhs%s_, mold=rhs%s_)
+
+    if (me==1) then
+      halo_west = rhs%s_(1,:,:)
+    else
+      halo_west = rhs[me-1]%s_(ubound(rhs[me-1]%s_,1),:,:)
+    end if
+    i = my_internal_west
+    call assert(i+1<=my_nx,"laplacian: westernmost subdomain too small")
+    do concurrent(j=2:ny-1, k=2:nz-1)
+      laplacian_rhs%s_(i,j,k) = ( halo_west(j,k  ) - 2*rhs%s_(i,j,k) + rhs%s_(i+1,j  ,k  ))/dx_**2 + &
+                                (rhs%s_(i,j-1,k  ) - 2*rhs%s_(i,j,k) + rhs%s_(i  ,j+1,k  ))/dy_**2 + &
+                                (rhs%s_(i,j  ,k-1) - 2*rhs%s_(i,j,k) + rhs%s_(i  ,j  ,k+1))/dz_**2
+    end do
+
+    do concurrent(i=my_internal_west+1:my_internal_east-1, j=2:ny-1, k=2:nz-1)
+      laplacian_rhs%s_(i,j,k) = (rhs%s_(i-1,j  ,k  ) - 2*rhs%s_(i,j,k) + rhs%s_(i+1,j  ,k  ))/dx_**2 + &
+                                (rhs%s_(i  ,j-1,k  ) - 2*rhs%s_(i,j,k) + rhs%s_(i  ,j+1,k  ))/dy_**2 + &
+                                (rhs%s_(i  ,j  ,k-1) - 2*rhs%s_(i,j,k) + rhs%s_(i  ,j  ,k+1))/dz_**2
+    end do
+
+    if (me==1) then
+      halo_east = rhs%s_(1,:,:)
+    else
+      halo_east = rhs[me+1]%s_(lbound(rhs[me+1]%s_,1),:,:)
+    end if
+    i = my_internal_east
+    call assert(i-1>0,"laplacian: easternmost subdomain too small")
+    do concurrent(j=2:ny-1, k=2:nz-1)
+      laplacian_rhs%s_(i,j,k) = (rhs%s_(i-1,j  ,k  ) - 2*rhs%s_(i,j,k) +  halo_east(j  ,k  ))/dx_**2 + &
+                                (rhs%s_(i  ,j-1,k  ) - 2*rhs%s_(i,j,k) + rhs%s_(i  ,j+1,k  ))/dy_**2 + &
+                                (rhs%s_(i  ,j  ,k-1) - 2*rhs%s_(i,j,k) + rhs%s_(i  ,j  ,k+1))/dz_**2
+    end do
+
+    laplacian_rhs%s_(:, 1,:) = 0.
+    laplacian_rhs%s_(:,ny,:) = 0.
+    laplacian_rhs%s_(:,:, 1) = 0.
+    laplacian_rhs%s_(:,:,nz) = 0.
+    if (me==1) laplacian_rhs%s_(1,:,:) = 0.
+    if (me==num_subdomains) laplacian_rhs%s_(my_nx,:,:) = 0.
+
+  end function
+
+end submodule subdomain_s

src/matcha/t_cell_collection_m.f90

@@ -42,8 +42,9 @@ pure module function construct(positions, time) result(t_cell_collection)
 
   interface t_cell_collection_bind_C_t
 
-    elemental module function construct_bind_C(t_cell_collection) result(t_cell_collection_bind_C)
+    impure elemental module function construct_bind_C(t_cell_collection) result(t_cell_collection_bind_C)
       !! Result is bind(C) representation of the data inside a t_cell_collection_t object
+      !! This function is impure because it invokes c_loc. Fortran 2023 compliance will allow this function to be pure.
       implicit none
       type(t_cell_collection_t), intent(in), target :: t_cell_collection
       type(t_cell_collection_bind_C_t) t_cell_collection_bind_C
@@ -60,7 +61,6 @@ pure module function positions(self) result(my_positions)
       double precision, allocatable :: my_positions(:,:)
     end function
 
-
     elemental module function time(self) result(my_time)
       !! Return the t_cell_collection_t object's time stamp
       implicit none

src/matcha_s.f90 → src/matcha_s.F90

@@ -32,7 +32,9 @@
         associate(me => this_image())
           associate(my_num_cells => data_partition%last(me) - data_partition%first(me) + 1)
 
+#ifndef NAGFOR
             call random_init(repeatable=.true., image_distinct=.true.)
+#endif
 
             allocate(random_positions(my_num_cells,ndim))
             call random_number(random_positions)