Generic and Rank Independent Programming

module heatflow_kinds
  use, intrinsic :: iso_fortran_env, only: real32, real64
  implicit none
  private
  public :: sp, dp
  integer, parameter :: sp = real32
  integer, parameter :: dp = real64
end module heatflow_kinds

module heatflow_norms
  use heatflow_kinds, only: sp, dp
  implicit none
  private
  public :: l2_norm

  ! One public name, one specific per (kind, rank) pair: 2 kinds x 2 ranks.
  interface l2_norm
    module procedure l2_norm_sp_1, l2_norm_dp_1, l2_norm_sp_2, l2_norm_dp_2
  end interface l2_norm

contains

  pure function l2_norm_sp_1(u) result(r)
    real(sp), intent(in) :: u(:)
    real(sp) :: r
    r = sqrt(sum(u**2))
  end function l2_norm_sp_1

  pure function l2_norm_dp_1(u) result(r)
    real(dp), intent(in) :: u(:)
    real(dp) :: r
    r = sqrt(sum(u**2))
  end function l2_norm_dp_1

  pure function l2_norm_sp_2(u) result(r)
    real(sp), intent(in) :: u(:,:)
    real(sp) :: r
    r = sqrt(sum(u**2))
  end function l2_norm_sp_2

  pure function l2_norm_dp_2(u) result(r)
    real(dp), intent(in) :: u(:,:)
    real(dp) :: r
    r = sqrt(sum(u**2))
  end function l2_norm_dp_2
end module heatflow_norms

program ex01_generic_by_kind_and_rank
  use heatflow_kinds, only: sp, dp
  use heatflow_norms, only: l2_norm
  implicit none
  real(dp) :: field(2,2)
  real(sp) :: row(3)
  field = reshape([3.0_dp, 0.0_dp, 0.0_dp, 4.0_dp], [2, 2])
  row = [1.0_sp, 2.0_sp, 2.0_sp]
  print '(a, f6.3)', 'rank 2, dp: ', l2_norm(field)
  print '(a, f6.3)', 'rank 1, sp: ', l2_norm(row)
  print '(a, f6.3)', 'rank 1, dp: ', l2_norm(real(row, dp))
  ! l2_norm(real(1, dp)) would be rejected: no specific accepts a scalar.
end program ex01_generic_by_kind_and_rank

rank 2, dp:  5.000
rank 1, sp:  3.000
rank 1, dp:  3.000

module heatflow_kinds
  use, intrinsic :: iso_fortran_env, only: real32, real64, int32
  implicit none
  private
  public :: sp, dp, ip
  integer, parameter :: sp = real32, dp = real64, ip = int32
end module heatflow_kinds

module heatflow_unlimited
  use heatflow_kinds, only: sp, dp, ip
  implicit none
  private
  public :: celsius_to_kelvin, describe

contains

  ! elemental: written for scalars, works for every rank (and is pure).
  elemental function celsius_to_kelvin(t) result(k)
    real(dp), intent(in) :: t
    real(dp) :: k
    k = t + 273.15_dp
  end function celsius_to_kelvin

  ! class(*): one specific for every type, at the price of select type
  ! and of a run-time dispatch.
  function describe(x) result(s)
    class(*), intent(in) :: x
    character(:), allocatable :: s
    character(len=32) :: buf
    select type (x)
    type is (integer(ip))
      write (buf, '(i0)') x
      s = 'integer ' // trim(buf)
    type is (real(sp))
      write (buf, '(f0.2)') x
      s = 'real(sp) ' // trim(buf)
    type is (real(dp))
      write (buf, '(f0.2)') x
      s = 'real(dp) ' // trim(buf)
    type is (character(*))
      s = 'character "' // x // '"'
    class default
      s = 'unsupported type'
    end select
  end function describe
end module heatflow_unlimited

program ex02_elemental_and_class_star
  use heatflow_kinds, only: sp, dp, ip
  use heatflow_unlimited, only: celsius_to_kelvin, describe
  implicit none
  real(dp) :: t0, t1(3), t2(2,2)
  t0 = 20.0_dp
  t1 = [0.0_dp, 20.0_dp, 100.0_dp]
  t2 = reshape([0.0_dp, 10.0_dp, 20.0_dp, 30.0_dp], [2, 2])
  print '(f7.2)', celsius_to_kelvin(t0)
  print '(3f8.2)', celsius_to_kelvin(t1)
  print '(2f8.2)', celsius_to_kelvin(t2)
  print '(a)', describe(42_ip)
  print '(a)', describe(1.5_sp)
  print '(a)', describe(2.5_dp)
  print '(a)', describe('heat')
  print '(a)', describe(.true.)
end program ex02_elemental_and_class_star

 293.15
  273.15  293.15  373.15
  273.15  283.15
  293.15  303.15
integer 42
real(sp) 1.50
real(dp) 2.50
character "heat"
unsupported type

! stats_body.inc -- included inside a module that has defined the kind "wp".
pure function mean_1d(u) result(m)
  real(wp), intent(in) :: u(:)
  real(wp) :: m
  m = sum(u) / real(size(u), wp)
end function mean_1d

pure function range_1d(u) result(r)
  real(wp), intent(in) :: u(:)
  real(wp) :: r
  r = maxval(u) - minval(u)
end function range_1d

module heatflow_kinds
  use, intrinsic :: iso_fortran_env, only: real32, real64
  implicit none
  private
  public :: sp, dp
  integer, parameter :: sp = real32, dp = real64
end module heatflow_kinds

! The same body is compiled twice, once per kind: the include file is the template.
module heatflow_stats_sp
  use heatflow_kinds, only: wp => sp
  implicit none
  private
  public :: mean_1d, range_1d
contains
  include 'stats_body.inc'
end module heatflow_stats_sp

module heatflow_stats_dp
  use heatflow_kinds, only: wp => dp
  implicit none
  private
  public :: mean_1d, range_1d
contains
  include 'stats_body.inc'
end module heatflow_stats_dp

! A third module renames the clashing specifics and merges them into generics.
module heatflow_stats
  use heatflow_stats_sp, only: mean_sp => mean_1d, range_sp => range_1d
  use heatflow_stats_dp, only: mean_dp => mean_1d, range_dp => range_1d
  implicit none
  private
  public :: mean, range_of
  interface mean
    module procedure mean_sp, mean_dp
  end interface mean
  interface range_of
    module procedure range_sp, range_dp
  end interface range_of
end module heatflow_stats

program ex03_include_template
  use heatflow_kinds, only: sp, dp
  use heatflow_stats, only: mean, range_of
  implicit none
  real(sp) :: a(4) = [1.0_sp, 2.0_sp, 3.0_sp, 6.0_sp]
  real(dp) :: b(3) = [10.0_dp, 20.0_dp, 60.0_dp]
  print '(2(a, f6.2))', 'sp mean=', mean(a), '  range=', range_of(a)
  print '(2(a, f6.2))', 'dp mean=', mean(b), '  range=', range_of(b)
end program ex03_include_template

sp mean=  3.00  range=  5.00
dp mean= 30.00  range= 50.00

! Preprocessor genericity: the file suffix .F90 makes the compiler run cpp first.
module heatflow_kinds
  use, intrinsic :: iso_fortran_env, only: real32, real64
  implicit none
  private
  public :: sp, dp
  integer, parameter :: sp = real32, dp = real64
end module heatflow_kinds

module heatflow_scale
  use heatflow_kinds, only: sp, dp
  implicit none
  private
  public :: scaled
  interface scaled
    module procedure scaled_sp, scaled_dp
  end interface scaled

#define DEFINE_SCALED(NAME, KIND) \
  pure function NAME(u, f) result(r); \
    real(KIND), intent(in) :: u(:), f; \
    real(KIND) :: r(size(u)); \
    r = f * u; \
  end function NAME

contains
  DEFINE_SCALED(scaled_sp, sp)
  DEFINE_SCALED(scaled_dp, dp)
end module heatflow_scale

program ex04_cpp_template
  use heatflow_kinds, only: sp, dp
  use heatflow_scale, only: scaled
  implicit none
  print '(3f6.1)', scaled([1.0_sp, 2.0_sp, 3.0_sp], 2.0_sp)
  print '(3f6.1)', scaled([1.0_dp, 2.0_dp, 3.0_dp], 0.5_dp)
end program ex04_cpp_template

   2.0   4.0   6.0
   0.5   1.0   1.5

module heatflow_kinds
  use, intrinsic :: iso_fortran_env, only: real64
  implicit none
  private
  public :: dp
  integer, parameter :: dp = real64
end module heatflow_kinds

module heatflow_inspect
  use heatflow_kinds, only: dp
  implicit none
  private
  public :: summarize, total

contains

  ! An assumed-rank dummy has the rank of its actual argument; the rank can be
  ! queried, but the elements can only be reached through select rank.
  subroutine summarize(a)
    real(dp), intent(in) :: a(..)
    integer :: k
    print '(a, i0, a, i0)', 'rank=', rank(a), ' size=', size(a)
    do k = 1, rank(a)
      print '(a, i0, a, i0, a, i0)', '  dim ', k, ': lbound=', lbound(a, k), ' extent=', size(a, k)
    end do
  end subroutine summarize

  function total(a) result(s)
    real(dp), intent(in) :: a(..)
    real(dp) :: s
    select rank (a)
    rank (0)
      s = a
    rank (1)
      s = sum(a)
    rank (2)
      s = sum(a)
    rank (3)
      s = sum(a)
    rank default
      error stop 'total: rank above 3 is not supported'
    end select
  end function total
end module heatflow_inspect

program ex05_assumed_rank
  use heatflow_kinds, only: dp
  use heatflow_inspect, only: summarize, total
  implicit none
  real(dp) :: scalar, line(5), plane(3, 4), block(2, 2, 2)
  real(dp), allocatable :: w(:, :)
  scalar = 7.0_dp
  line = 1.0_dp
  plane = 2.0_dp
  block = 0.5_dp
  allocate (w(0:2, -1:1))
  w = 1.0_dp
  call summarize(scalar)
  call summarize(plane)
  call summarize(w)
  call summarize(line(2:4))
  print '(4(a, f6.2))', 'scalar ', total(scalar), '  line ', total(line), '  plane ', total(plane), '  block ', total(block)
end program ex05_assumed_rank

rank=0 size=1
rank=2 size=12
  dim 1: lbound=1 extent=3
  dim 2: lbound=1 extent=4
rank=2 size=9
  dim 1: lbound=1 extent=3
  dim 2: lbound=1 extent=3
rank=1 size=3
  dim 1: lbound=1 extent=3
scalar   7.00  line   5.00  plane  24.00  block   4.00

module heatflow_kinds
  use, intrinsic :: iso_fortran_env, only: real32, real64
  implicit none
  private
  public :: sp, dp
  integer, parameter :: sp = real32, dp = real64
end module heatflow_kinds

module heatflow_typeof
  use heatflow_kinds, only: sp, dp
  implicit none
  private
  public :: grid_t, mean_of, clone

  type :: grid_t
    integer :: nx = 0
    real(dp), allocatable :: u(:, :)
  end type grid_t

  interface mean_of
    module procedure mean_sp, mean_dp
  end interface mean_of

contains

  ! typeof(u): the accumulator and the result follow the kind of u. The
  ! declaration is written once; only the dummy argument names the kind.
  function mean_sp(u) result(m)
    real(sp), intent(in) :: u(:)
    typeof(u) :: m
    typeof(u) :: acc
    acc = sum(u)
    m = acc / size(u)
  end function mean_sp

  function mean_dp(u) result(m)
    real(dp), intent(in) :: u(:)
    typeof(u) :: m
    typeof(u) :: acc
    acc = sum(u)
    m = acc / size(u)
  end function mean_dp

  ! classof(g): polymorphic with the declared type of g; with class(grid_t)
  ! this would be written twice, here it follows the argument.
  function clone(g) result(copy)
    class(grid_t), intent(in) :: g
    classof(g), allocatable :: copy
    allocate (copy, source=g)
  end function clone
end module heatflow_typeof

program ex07_typeof_classof
  use heatflow_kinds, only: sp, dp
  use heatflow_typeof, only: grid_t, mean_of, clone
  implicit none
  type(grid_t) :: g
  class(grid_t), allocatable :: h
  g%nx = 2
  allocate (g%u(2, 2))
  g%u = 5.0_dp
  h = clone(g)
  h%u = 1.0_dp
  print '(a, f5.1, a, f5.1)', 'original ', g%u(1, 1), '  clone ', h%u(1, 1)
  print '(a, f5.2)', 'mean sp ', mean_of([1.0_sp, 2.0_sp, 6.0_sp])
  print '(a, f5.2)', 'mean dp ', mean_of([1.0_dp, 2.0_dp, 6.0_dp])
end program ex07_typeof_classof

original   5.0  clone   1.0
mean sp  3.00
mean dp  3.00

! F2023 typeof: accepted by Flang 23, rejected by gfortran 16.2.
program err05_typeof
  implicit none
  real(kind(1.0d0)) :: x(3) = [1.0d0, 2.0d0, 3.0d0]
  typeof(x) :: y
  y = x(2)
  print *, y
end program err05_typeof

err05_typeof.f90:5:7:

    5 |   typeof(x) :: y
      |       1
Error: Mangled derived type definition at (1)
err05_typeof.f90:6:3:

    6 |   y = x(2)
      |   1
Error: Symbol 'y' at (1) has no IMPLICIT type

! F2023 classof: accepted by Flang 23, rejected by gfortran 16.2.
module shapes
  implicit none
  type :: shape_t
    real :: s = 1.0
  end type shape_t
contains
  function clone(a) result(b)
    class(shape_t), intent(in) :: a
    classof(a), allocatable :: b
    allocate (b, source=a)
  end function clone
end module shapes

err06_classof.f90:10:5:

   10 |     classof(a), allocatable :: b
      |     1
Error: Unclassifiable statement at (1)
err06_classof.f90:11:15:

   11 |     allocate (b, source=a)
      |               1
Error: Allocate-object at (1) is neither a data pointer nor an allocatable variable
err06_classof.f90:8:29:

    8 |   function clone(a) result(b)
      |                             1
Error: Result 'b' of contained function 'clone' at (1) has no IMPLICIT type

module heatflow_kinds
  use, intrinsic :: iso_fortran_env, only: real64
  implicit none
  private
  public :: dp
  integer, parameter :: dp = real64
end module heatflow_kinds

module heatflow_rank
  use heatflow_kinds, only: dp
  implicit none
  private
  public :: relabel

contains

  ! rank(rank(ref)): the dummy and the allocatable result have the rank of ref,
  ! decided by a constant expression instead of a column of colons.
  subroutine relabel(ref, copy)
    real(dp), rank(2), intent(in) :: ref
    real(dp), rank(rank(ref)), allocatable, intent(out) :: copy
    allocate (copy, mold=ref)
    copy = -ref
  end subroutine relabel
end module heatflow_rank

program ex08_rank_clause
  use heatflow_kinds, only: dp
  use heatflow_rank, only: relabel
  implicit none
  real(dp) :: u(0:1, 0:2)
  real(dp), allocatable :: v(:, :)
  u = 3.0_dp
  call relabel(u, v)
  print '(a, 2i2, a, f5.1)', 'shape(v) =', shape(v), '  v(1,1) =', v(1, 1)
end program ex08_rank_clause

shape(v) = 2 3  v(1,1) = -3.0

! F2023 rank clause: accepted by Flang 23, rejected by gfortran 16.2.
program err07_rank_clause
  implicit none
  real, rank(3), allocatable :: x
  allocate (x(2, 2, 2))
  print *, rank(x)
end program err07_rank_clause

err07_rank_clause.f90:4:8:

    4 |   real, rank(3), allocatable :: x
      |        1
Error: Invalid character in name at (1)
err07_rank_clause.f90:5:13:

    5 |   allocate (x(2, 2, 2))
      |             1
Error: Allocate-object at (1) is neither a data pointer nor an allocatable variable
err07_rank_clause.f90:6:17:

    6 |   print *, rank(x)
      |                 1
Error: Symbol 'x' at (1) has no IMPLICIT type
err07_rank_clause.f90:6:17:

    6 |   print *, rank(x)
      |                 1
Error: Symbol 'x' at (1) has no IMPLICIT type

! UNCHECKED F2023 SYNTAX: multiple subscripts (9.5.3.2). Rejected by gfortran 16.2 and Flang 23.
program err01_at_subscript
  implicit none
  real :: a(3, 3, 3)
  integer, parameter :: idx(3) = [2, 3, 1]
  a = 0.0
  a(@idx) = 1.0
  print *, sum(a(@[1, 1]:[2, 2], 1))
end program err01_at_subscript

err01_at_subscript.f90:7:5:

    7 |   a(@idx) = 1.0
      |     1
Error: Expected array subscript at (1)
err01_at_subscript.f90:8:18:

    8 |   print *, sum(a(@[1, 1]:[2, 2], 1))
      |                  1
Error: Expected array subscript at (1)

error: Semantic errors in err01_at_subscript.f90
err01_at_subscript.f90:7:5: error: No explicit type declared for '@idx'
    a(@idx) = 1.0
      ^^^^
err01_at_subscript.f90:8:18: error: No explicit type declared for '@'
    print *, sum(a(@[1, 1]:[2, 2], 1))
                   ^
err01_at_subscript.f90:8:26: error: Must be a scalar value, but is a rank-1 array
    print *, sum(a(@[1, 1]:[2, 2], 1))
                           ^^^^^^

! UNCHECKED F2023 SYNTAX: ALLOCATE with a bounds array (9.7.1.1). Rejected by gfortran 16.2 and Flang 23.
program err02_allocate_bounds
  implicit none
  real, allocatable :: a(:, :, :)
  integer :: sh(3)
  sh = [2, 3, 4]
  allocate (a(sh))
  print *, shape(a)
end program err02_allocate_bounds

err02_allocate_bounds.f90:7:13:

    7 |   allocate (a(sh))
      |             1
Error: Rank mismatch in array reference at (1) (1/3)

error: Semantic errors in err02_allocate_bounds.f90
err02_allocate_bounds.f90:7:13: error: The number of shape specifications, when they appear, must match the rank of allocatable object
    allocate (a(sh))
              ^
err02_allocate_bounds.f90:4:24: Declared here with rank 3
    real, allocatable :: a(:, :, :)
                         ^
err02_allocate_bounds.f90:7:15: error: Must be a scalar value, but is a rank-1 array
    allocate (a(sh))
                ^^

! UNCHECKED F2023 SYNTAX: pointer bounds from arrays (10.2.2.3). Rejected by gfortran 16.2 and Flang 23.
program err03_pointer_bounds
  implicit none
  real, target :: t(24)
  real, pointer :: p(:, :, :)
  integer, parameter :: lb(3) = [0, 0, 0], ub(3) = [1, 2, 3]
  p(lb:ub) => t
  print *, shape(p)
end program err03_pointer_bounds

err03_pointer_bounds.f90:7:4:

    7 |   p(lb:ub) => t
      |    1
Error: Array index at (1) must be scalar

error: Semantic errors in err03_pointer_bounds.f90
err03_pointer_bounds.f90:7:5: error: Must be a scalar value, but is a rank-1 array
    p(lb:ub) => t
      ^^
err03_pointer_bounds.f90:7:8: error: Must be a scalar value, but is a rank-1 array
    p(lb:ub) => t
         ^^

! UNCHECKED F2023 SYNTAX: bounds from an array in a declaration (8.5.8.2). Rejected by gfortran 16.2 and Flang 23.
program err04_shape_decl
  implicit none
  integer, parameter :: sh(3) = [2, 3, 4]
  real :: a(sh)
  print *, shape(a)
end program err04_shape_decl

err04_shape_decl.f90:5:15:

    5 |   real :: a(sh)
      |               1
Error: Explicit shaped array with nonconstant bounds at (1)
err04_shape_decl.f90:6:18:

error: Semantic errors in err04_shape_decl.f90
err04_shape_decl.f90:5:3: error: not yet implemented: TODO: Analyze overload for ExplicitShapeBoundsSpec
    real :: a(sh)
    ^^^^^^^^^^^^^