Interfaces, Generics and Operator Overloading

! NON-CONFORMING ON PURPOSE: an external procedure with an implicit interface
! cannot be called with argument keywords (15.4.2.2).
program implicit_keyword
  implicit none
  external :: smooth
  real :: u(3)
  u = [1.0, 2.0, 3.0]
  call smooth(n=3, u=u)
end program implicit_keyword

err_04_implicit_keyword.f90:8:16:

    8 |   call smooth(n=3, u=u)
      |                1
Error: Keyword argument requires explicit interface for procedure 'smooth' at (1)

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_blas
  use heatflow_kinds, only: dp
  implicit none
  private
  public :: ddot

  ! An interface body gives an external procedure an explicit interface.
  ! BLAS is written in Fortran, so no bind(c) is needed; the dummy argument
  ! names and the intents document the contract of the routine.
  interface
    function ddot(n, x, incx, y, incy) result(d)
      import :: dp
      integer, intent(in) :: n, incx, incy
      real(dp), intent(in) :: x(*), y(*)
      real(dp) :: d
    end function ddot
  end interface
end module heatflow_blas

program interface_body_demo
  use heatflow_kinds, only: dp
  use heatflow_blas, only: ddot
  implicit none
  real(dp) :: u(4), v(4)

  u = [1.0_dp, 2.0_dp, 3.0_dp, 4.0_dp]
  v = [4.0_dp, 3.0_dp, 2.0_dp, 1.0_dp]
  print '(a, f6.1)', 'ddot(u, v)        = ', ddot(size(u), u, 1, v, 1)
  print '(a, f6.1)', 'with keywords     = ', ddot(x=u, y=v, n=4, incx=1, incy=1)
end program interface_body_demo

ddot(u, v)        =   20.0
with keywords     =   20.0

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_stats
  use heatflow_kinds, only: sp, dp
  implicit none
  private
  public :: rms

  ! A generic interface block: the generic name is rms, the specifics keep
  ! their own names (rms_1d, rms_2d, rms_1d_sp) and stay callable.
  interface rms
    module procedure rms_1d
    module procedure rms_2d
    module procedure rms_1d_sp
  end interface rms

  ! The GENERIC statement (Fortran 2023) does the same without a block.
  generic, public :: mean_value => mean_1d, mean_2d

contains

  pure function rms_1d(x) result(r)
    real(dp), intent(in) :: x(:)
    real(dp) :: r
    r = sqrt(sum(x**2)/real(size(x), dp))
  end function rms_1d

  pure function rms_2d(x) result(r)
    real(dp), intent(in) :: x(:, :)
    real(dp) :: r
    r = sqrt(sum(x**2)/real(size(x), dp))
  end function rms_2d

  pure function rms_1d_sp(x) result(r)
    real(sp), intent(in) :: x(:)
    real(sp) :: r
    r = sqrt(sum(x**2)/real(size(x), sp))
  end function rms_1d_sp

  pure function mean_1d(x) result(r)
    real(dp), intent(in) :: x(:)
    real(dp) :: r
    r = sum(x)/real(size(x), dp)
  end function mean_1d

  pure function mean_2d(x) result(r)
    real(dp), intent(in) :: x(:, :)
    real(dp) :: r
    r = sum(x)/real(size(x), dp)
  end function mean_2d

end module heatflow_stats

program generic_rms_demo
  use heatflow_kinds, only: sp, dp
  use heatflow_stats, only: rms, mean_value
  implicit none
  real(dp) :: row(4), plate(3, 3)
  real(sp) :: row_sp(4)
  integer :: k

  row = [1.0_dp, 2.0_dp, 3.0_dp, 4.0_dp]
  plate = reshape([(real(k, dp), k=1, 9)], [3, 3])
  row_sp = real(row, sp)

  print '(a, f8.4)', 'rms of a rank-1 dp array : ', rms(row)
  print '(a, f8.4)', 'rms of a rank-2 dp array : ', rms(plate)
  print '(a, f8.4)', 'rms of a rank-1 sp array : ', rms(row_sp)
  print '(a, f8.4)', 'mean of a rank-1 dp array: ', mean_value(row)
  print '(a, f8.4)', 'mean of a rank-2 dp array: ', mean_value(plate)
end program generic_rms_demo

rms of a rank-1 dp array :   2.7386
rms of a rank-2 dp array :   5.6273
rms of a rank-1 sp array :   2.7386
mean of a rank-1 dp array:   2.5000
mean of a rank-2 dp array:   5.0000

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_adjust
  use heatflow_kinds, only: dp
  implicit none
  private
  public :: adjust

  interface adjust
    module procedure adjust_by_factor
    module procedure adjust_by_passes
    module procedure adjust_vector
    module procedure adjust_elemental
  end interface adjust

contains

  ! The second argument differs in TYPE (real versus integer), so the two
  ! specifics are distinguishable by position and by keyword.
  pure function adjust_by_factor(t, factor) result(r)
    real(dp), intent(in) :: t
    real(dp), intent(in) :: factor
    real(dp) :: r
    r = t*factor
  end function adjust_by_factor

  pure function adjust_by_passes(t, passes) result(r)
    real(dp), intent(in) :: t
    integer, intent(in) :: passes
    real(dp) :: r
    r = t*0.5_dp**passes
  end function adjust_by_passes

  ! A rank-1 specific: preferred over the elemental one for rank-1 actuals.
  pure function adjust_vector(t) result(r)
    real(dp), intent(in) :: t(:)
    real(dp) :: r(size(t))
    r = t - sum(t)/real(size(t), dp)
  end function adjust_vector

  ! An elemental specific: used for every other rank.
  elemental function adjust_elemental(t) result(r)
    real(dp), intent(in) :: t
    real(dp) :: r
    r = max(t, 0.0_dp)
  end function adjust_elemental

end module heatflow_adjust

program resolution_demo
  use heatflow_kinds, only: dp
  use heatflow_adjust, only: adjust
  implicit none
  real(dp) :: line(3), plate(2, 2)

  line = [10.0_dp, 20.0_dp, 60.0_dp]
  plate = reshape([-1.0_dp, 2.0_dp, -3.0_dp, 4.0_dp], [2, 2])

  print '(a, f6.1)', 'factor= : ', adjust(10.0_dp, factor=3.0_dp)
  print '(a, f6.1)', 'passes= : ', adjust(10.0_dp, passes=2)
  print '(a, f6.1)', 'by type : ', adjust(8.0_dp, 3)
  print '(a, 3f7.1)', 'rank 1  : ', adjust(line)
  print '(a, 4f7.1)', 'rank 2  : ', adjust(plate)
  print '(a, f6.1)', 'scalar  : ', adjust(-5.0_dp)
end program resolution_demo

factor= :   30.0
passes= :    2.5
by type :    1.0
rank 1  :   -20.0  -10.0   30.0
rank 2  :     0.0    2.0    0.0    4.0
scalar  :    0.0

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_flux
  use heatflow_kinds, only: dp
  implicit none
  private
  public :: flux_t, operator(+), operator(-), operator(*), operator(==)
  public :: operator(.dot.), operator(.norm.)

  real(dp), parameter :: tolerance = 1.0e-12_dp

  ! Heat flux through one cell face: a two-component vector.
  type :: flux_t
    real(dp) :: qx = 0.0_dp
    real(dp) :: qy = 0.0_dp
  end type flux_t

  ! A generic name equal to the type name overrides nothing: the default
  ! structure constructor stays available and this adds a second one.
  interface flux_t
    module procedure flux_from_array
  end interface flux_t

  interface operator(+)           ! binary +
    module procedure flux_add
  end interface operator(+)

  interface operator(-)           ! unary and binary -: both are generic
    module procedure flux_negate
    module procedure flux_subtract
  end interface operator(-)

  interface operator(*)           ! scalar * flux and flux * scalar
    module procedure real_times_flux
    module procedure flux_times_real
  end interface operator(*)

  interface operator(==)
    module procedure flux_equal
  end interface operator(==)

  interface operator(.dot.)       ! a new binary operator
    module procedure flux_dot
  end interface operator(.dot.)

  interface operator(.norm.)      ! a new unary operator
    module procedure flux_norm
  end interface operator(.norm.)

contains

  pure function flux_from_array(q) result(f)
    real(dp), intent(in) :: q(2)
    type(flux_t) :: f
    f%qx = q(1)
    f%qy = q(2)
  end function flux_from_array

  pure function flux_add(a, b) result(f)
    type(flux_t), intent(in) :: a, b
    type(flux_t) :: f
    f = flux_t(a%qx + b%qx, a%qy + b%qy)
  end function flux_add

  pure function flux_negate(a) result(f)
    type(flux_t), intent(in) :: a
    type(flux_t) :: f
    f = flux_t(-a%qx, -a%qy)
  end function flux_negate

  pure function flux_subtract(a, b) result(f)
    type(flux_t), intent(in) :: a, b
    type(flux_t) :: f
    f = flux_t(a%qx - b%qx, a%qy - b%qy)
  end function flux_subtract

  pure function real_times_flux(s, a) result(f)
    real(dp), intent(in) :: s
    type(flux_t), intent(in) :: a
    type(flux_t) :: f
    f = flux_t(s*a%qx, s*a%qy)
  end function real_times_flux

  pure function flux_times_real(a, s) result(f)
    type(flux_t), intent(in) :: a
    real(dp), intent(in) :: s
    type(flux_t) :: f
    f = flux_t(s*a%qx, s*a%qy)
  end function flux_times_real

  pure function flux_equal(a, b) result(same)
    type(flux_t), intent(in) :: a, b
    logical :: same
    same = abs(a%qx - b%qx) <= tolerance .and. abs(a%qy - b%qy) <= tolerance
  end function flux_equal

  pure function flux_dot(a, b) result(d)
    type(flux_t), intent(in) :: a, b
    real(dp) :: d
    d = a%qx*b%qx + a%qy*b%qy
  end function flux_dot

  pure function flux_norm(a) result(n)
    type(flux_t), intent(in) :: a
    real(dp) :: n
    n = sqrt(a%qx**2 + a%qy**2)
  end function flux_norm

end module heatflow_flux

program flux_operators_demo
  use heatflow_kinds, only: dp
  use heatflow_flux
  implicit none
  type(flux_t) :: a, b

  a = flux_t(3.0_dp, 4.0_dp)          ! default structure constructor
  b = flux_t([1.0_dp, 2.0_dp])        ! the generic constructor

  call show('a + b          ', a + b)
  call show('-a             ', -a)
  call show('a - b          ', a - b)
  call show('2 * a          ', 2.0_dp*a)
  call show('a * 0.5        ', a*0.5_dp)
  print '(a, f6.1)', 'a .dot. b      = ', a .dot. b
  print '(a, f6.1)', '.norm. a       = ', .norm. a
  print '(a, l2)', 'a == a         = ', a == a
  print '(a, l2)', 'a .eq. b       = ', a .eq. b      ! the same operator
  ! Precedence is that of the operator's category (10.1.3): a defined unary
  ! operator binds tighter than +, a defined binary operator looser than *.
  print '(a, f6.1)', '.norm. a + 1   = ', .norm. a + 1.0_dp
  print '(a, f6.1)', 'a .dot. 2 * b  = ', a .dot. 2.0_dp*b

contains

  subroutine show(label, f)
    character(len=*), intent(in) :: label
    type(flux_t), intent(in) :: f
    print '(a, a, 2f6.1)', label, '= ', f%qx, f%qy
  end subroutine show

end program flux_operators_demo

a + b          =    4.0   6.0
-a             =   -3.0  -4.0
a - b          =    2.0   2.0
2 * a          =    6.0   8.0
a * 0.5        =    1.5   2.0
a .dot. b      =   11.0
.norm. a       =    5.0
a == a         =  T
a .eq. b       =  F
.norm. a + 1   =    6.0
a .dot. 2 * b  =   22.0

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_flux_tbp
  use heatflow_kinds, only: dp
  implicit none
  private
  public :: flux_t

  type :: flux_t
    real(dp) :: qx = 0.0_dp
    real(dp) :: qy = 0.0_dp
  contains
    procedure, private :: add
    procedure, private :: times_real
    procedure, pass(b), private :: real_times
    procedure, private :: from_array
    procedure, private :: from_scalar
    procedure, pass(f), private :: to_array
    ! A type-bound generic needs a passed-object dummy argument for every
    ! operator or assignment binding (C779); pass(b) and pass(f) say which.
    generic :: operator(+) => add
    generic :: operator(*) => times_real, real_times
    generic :: assignment(=) => from_array, from_scalar, to_array
  end type flux_t

contains

  pure function add(a, b) result(f)
    class(flux_t), intent(in) :: a
    type(flux_t), intent(in) :: b
    type(flux_t) :: f
    f = flux_t(a%qx + b%qx, a%qy + b%qy)
  end function add

  pure function times_real(a, s) result(f)
    class(flux_t), intent(in) :: a
    real(dp), intent(in) :: s
    type(flux_t) :: f
    f = flux_t(s*a%qx, s*a%qy)
  end function times_real

  pure function real_times(s, b) result(f)
    real(dp), intent(in) :: s
    class(flux_t), intent(in) :: b
    type(flux_t) :: f
    f = flux_t(s*b%qx, s*b%qy)
  end function real_times

  ! Defined assignment: the first argument is intent(out) or intent(inout),
  ! the second intent(in) or value; the right-hand side is passed in
  ! parentheses. A pure procedure cannot have a polymorphic intent(out)
  ! dummy argument, so the passed-object lhs is intent(inout) here.
  elemental subroutine from_scalar(lhs, rhs)
    class(flux_t), intent(inout) :: lhs
    real(dp), intent(in) :: rhs
    lhs%qx = rhs
    lhs%qy = rhs
  end subroutine from_scalar

  pure subroutine from_array(lhs, rhs)
    class(flux_t), intent(inout) :: lhs
    real(dp), intent(in) :: rhs(2)
    lhs%qx = rhs(1)
    lhs%qy = rhs(2)
  end subroutine from_array

  pure subroutine to_array(lhs, f)
    real(dp), intent(out) :: lhs(2)
    class(flux_t), intent(in) :: f
    lhs = [f%qx, f%qy]
  end subroutine to_array

end module heatflow_flux_tbp

program type_bound_demo
  use heatflow_kinds, only: dp
  use heatflow_flux_tbp, only: flux_t
  implicit none
  type(flux_t) :: a, b, cells(3)
  real(dp) :: raw(2)

  a = [3.0_dp, 4.0_dp]          ! defined assignment, array on the right
  b = 1.5_dp                    ! defined assignment, scalar on the right
  raw = a + b                   ! operator(+), then defined assignment to an array
  print '(a, 2f6.1)', 'raw   = a + b   : ', raw

  raw = 2.0_dp*a                ! operator(*) bound with pass(b)
  print '(a, 2f6.1)', 'raw   = 2 * a   : ', raw
  raw = a*0.5_dp                ! operator(*) bound with the default pass
  print '(a, 2f6.1)', 'raw   = a * 0.5 : ', raw

  cells = 0.5_dp               ! the elemental specific fills the whole array
  raw = cells(3)
  print '(a, 2f6.1)', 'cells(3)        : ', raw
end program type_bound_demo

raw   = a + b   :    4.5   5.5
raw   = 2 * a   :    6.0   8.0
raw   = a * 0.5 :    1.5   2.0
cells(3)        :    0.5   0.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_rates
  use heatflow_kinds, only: dp
  implicit none
  private
  public :: rate_fn, advance, newton_cooling, sqrt_cooling

  ! An abstract interface has no procedure behind it: it only names a set
  ! of characteristics that dummy arguments and procedure pointers can share.
  abstract interface
    function rate_fn(u) result(du)
      import :: dp
      real(dp), intent(in) :: u
      real(dp) :: du
    end function rate_fn
  end interface

contains

  ! Explicit Euler with a rate supplied by the caller (a dummy procedure).
  function advance(f, u0, dt, nsteps) result(u)
    procedure(rate_fn) :: f
    real(dp), intent(in) :: u0, dt
    integer, intent(in) :: nsteps
    real(dp) :: u
    integer :: n
    u = u0
    do n = 1, nsteps
      u = u + dt*f(u)
    end do
  end function advance

  ! A module procedure with the characteristics of rate_fn.
  function newton_cooling(u) result(du)
    real(dp), intent(in) :: u
    real(dp) :: du
    real(dp), parameter :: ambient = 20.0_dp, k = 0.5_dp
    du = -k*(u - ambient)
  end function newton_cooling

  ! The specific intrinsic names (DSQRT and friends) are obsolescent in
  ! Fortran 2023 (Annex B.3.12): wrap the generic intrinsic in a procedure
  ! of your own instead of passing the intrinsic itself.
  function sqrt_cooling(u) result(du)
    real(dp), intent(in) :: u
    real(dp) :: du
    du = -sqrt(u)
  end function sqrt_cooling

end module heatflow_rates

program dummy_procedure_demo
  use heatflow_kinds, only: dp
  use heatflow_rates, only: rate_fn, advance, newton_cooling, sqrt_cooling
  implicit none
  procedure(rate_fn), pointer :: chosen

  print '(a, f8.4)', 'module procedure   : ', advance(newton_cooling, 100.0_dp, 0.1_dp, 10)
  chosen => newton_cooling
  print '(a, f8.4)', 'procedure pointer  : ', advance(chosen, 100.0_dp, 0.1_dp, 10)
  print '(a, f8.4)', 'wrapped intrinsic  : ', advance(sqrt_cooling, 100.0_dp, 0.1_dp, 10)
end program dummy_procedure_demo

module procedure   :  67.8990
procedure pointer  :  67.8990
wrapped intrinsic  :  90.2256

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_hooks
  use heatflow_kinds, only: dp
  implicit none
  private
  public :: rate_fn, report_sub, probe_t, select_rate, constant_rate, linear_rate

  type :: probe_t
    character(len=12) :: label = 'probe'
    ! A procedure pointer component with NOPASS: no passed-object dummy.
    procedure(rate_fn), pointer, nopass :: rate => null()
    ! With PASS the object itself is passed as the named dummy argument;
    ! the interface (report_sub) has self as its first argument.
    procedure(report_sub), pointer, pass(self) :: report => null()
  end type probe_t

  ! The abstract interfaces come after the type: report_sub refers to probe_t
  ! and probe_t refers to report_sub, and the type is imported by name.
  abstract interface
    function rate_fn(u) result(du)
      import :: dp
      real(dp), intent(in) :: u
      real(dp) :: du
    end function rate_fn

    subroutine report_sub(self, value)
      import :: dp, probe_t
      class(probe_t), intent(in) :: self
      real(dp), intent(in) :: value
    end subroutine report_sub
  end interface

contains

  function constant_rate(u) result(du)
    real(dp), intent(in) :: u
    real(dp) :: du
    du = 0.0_dp*u - 1.0_dp
  end function constant_rate

  function linear_rate(u) result(du)
    real(dp), intent(in) :: u
    real(dp) :: du
    du = -0.5_dp*u
  end function linear_rate

  ! A function can return a procedure pointer.
  function select_rate(name) result(p)
    character(len=*), intent(in) :: name
    procedure(rate_fn), pointer :: p
    select case (name)
    case ('constant')
      p => constant_rate
    case ('linear')
      p => linear_rate
    case default
      p => null()
    end select
  end function select_rate

end module heatflow_hooks

module heatflow_printers
  use heatflow_kinds, only: dp
  use heatflow_hooks, only: probe_t
  implicit none
  private
  public :: print_with_label
contains
  subroutine print_with_label(self, value)
    class(probe_t), intent(in) :: self
    real(dp), intent(in) :: value
    print '(a, a, a, f7.2)', 'probe "', trim(self%label), '" reads ', value
  end subroutine print_with_label
end module heatflow_printers

program procedure_pointer_demo
  use heatflow_kinds, only: dp
  use heatflow_hooks, only: probe_t, rate_fn, select_rate, constant_rate, linear_rate
  use heatflow_printers, only: print_with_label
  implicit none
  type(probe_t) :: probe
  procedure(rate_fn), pointer :: p
  real(dp) :: u

  probe%label = 'centre'
  probe%rate => linear_rate                 ! pointer assignment to a component
  probe%report => print_with_label

  u = 80.0_dp
  call probe%report(probe%rate(u))          ! probe is passed as self; rate gets no object
  print '(a, l2)', 'rate points to linear_rate  :', associated(probe%rate, linear_rate)
  print '(a, l2)', 'rate points to constant_rate:', associated(probe%rate, constant_rate)

  p => select_rate('constant')
  print '(a, f7.2)', 'selected constant rate      :', p(u)
  p => select_rate('unknown')
  print '(a, l2)', 'unknown name gives a pointer:', associated(p)

  probe%rate => null()
  print '(a, l2)', 'after => null(), associated :', associated(probe%rate)
end program procedure_pointer_demo

probe "centre" reads  -40.00
rate points to linear_rate  : T
rate points to constant_rate: F
selected constant rate      :  -1.00
unknown name gives a pointer: F
after => null(), associated : F