Array Intrinsics and Reductions

This content was generated with the assistance of AI and should be verified against the official documentation before being relied on in production.

The reference material consulted while preparing these pages is listed in the section bibliography.

This page is written against Fortran 2023 (ISO/IEC 1539-1:2023, J3 final draft J3/23-007r1) and was checked on 2026-10-18 with GNU Fortran 16.2 and LLVM Flang 23. The array intrinsic functions are what lets Fortran code say "the sum of the temperatures", "where is the hottest cell" or "shift the field one cell to the left" without writing a loop, and a compiler is free to turn each such statement into the best loop, vectorised or not, for the machine. The page covers the whole family in the order in which a numerical program meets it: the rules shared by all of them (arguments, dim, mask, kind, array element order), the reductions, the location functions, the inquiry functions, the construction and manipulation functions, and matmul and dot_product. It ends with the diagnostics each compiler gives for misuse, the cost of the temporaries that some of these functions create, and a table of where each function came from.

The arrays themselves (declaration, sections, constructors, whole-array operations and where) are on Arrays Fundamentals; allocatable arrays on Allocatable Arrays and Dynamic Memory; how array arguments reach a procedure on Array Arguments: Assumed Shape, Size and Rank; and the catalogue of all other intrinsic procedures on Intrinsic Procedures Overview.

Every example belongs to the heatflow scenario (a two-dimensional heat-diffusion solver) and uses its names (dp for real64, a temperature field u), but each one is a complete, self-contained file that can be extracted and run on its own. Unless a paragraph says otherwise, every example compiles without warnings under gfortran -std=f2023 -Wall -Wextra -pedantic -fimplicit-none -fcheck=all and was also built and run with LLVM Flang 23, and the output shown is what it printed. Examples that are deliberately wrong are labelled as non-conforming in the text and in a comment on their first line, and the compiler messages or run-time failures they produce are shown.

The rules shared by every array intrinsic

Most of this page concerns transformational functions: functions that take whole arrays and return a value (or an array) computed from all of them, rather than working element by element. Three other classes appear alongside them. Inquiry functions (size, shape, lbound, ubound, rank, is_contiguous) return a property of an argument and do not even need its value to be defined; elemental functions (merge, and the bit and mathematical functions that are applied to arrays) work on each element independently; and a few array intrinsics are subroutines. The classes are listed for every procedure in J3/23-007r1, Table 16.1, and the specification of each function is in 16.9.

Source: J3/23-007r1 (Fortran 2023), 16.2 Arguments to intrinsic procedures and 16.7 Standard generic intrinsic procedures

The rules that apply to all of them:

  • Argument keywords. Every intrinsic can be called with positional or keyword arguments (16.2.1), using the names in the standard (array, dim, mask, kind, back, vector, shift, boundary, ncopies, …​). Because many functions have several optional arguments, keywords are the readable way to call them: maxloc(u, back=.true.), eoshift(v, 2, boundary=-1).

  • The dim argument selects one dimension to collapse. For an array of rank n, dim is an integer scalar from 1 to n. A reduction called with dim returns an array of rank n - 1 whose shape is that of the argument with dimension dim removed; without dim the result is a scalar (or, for the location functions, a rank-one array of subscripts). A rank-one argument with dim gives a scalar (16.2.4). A dim that is not a valid dimension is non-conforming; a compiler rejects a constant one and, with run-time checking, a variable one (see the diagnostics section).

  • The mask argument selects the elements that take part. It is a logical array conformable with the argument (same shape, or scalar). Elements that the mask does not select need not even be defined (16.2.3), but the mask does not change how the arguments, which may be array expressions, are evaluated before the call. sum(u, mask=m) is defined to equal sum(pack(u, m)); the function form just avoids the temporary.

  • The kind argument sets the kind of an integer result. count, maxloc, minloc, findloc, size, shape, lbound and ubound return default integers unless kind= names another integer kind; it must be a constant expression. Use kind=int64 for arrays whose extents or subscripts can exceed huge(0).

  • Array element order. When a function reads or writes an array "in order" (reshape, pack, unpack, and reduce unless ordered= says otherwise) it uses array element order: the first subscript varies fastest, so the elements of u(4, 3) are u(1,1), u(2,1), u(3,1), u(4,1), u(1,2), and so on (9.5.3.3). This is the memory order of Fortran arrays, and the opposite of C.

  • Shapes must be defined. The arguments of transformational and elemental functions must have a defined shape (16.2.2), so an unallocated allocatable array or a disassociated pointer cannot be passed.

The table lists every function of the family with its class, the argument that gets reduced or shifted, and the shape of the result. Source names the clause of J3/23-007r1 where the function is specified.

Function Class Result Source

sum, product

transformational

scalar, or rank n - 1 with dim

16.9.201, 16.9.165

maxval, minval

transformational

scalar, or rank n - 1 with dim

16.9.138, 16.9.144

count, any, all, parity

transformational

scalar (integer or logical), or rank n - 1 with dim

16.9.66, 16.9.15, 16.9.12, 16.9.159

iall, iany, iparity

transformational

scalar integer, or rank n - 1 with dim

16.9.99, 16.9.101, 16.9.112

norm2

transformational

scalar real, or rank n - 1 with dim

16.9.153

reduce

transformational

scalar of the array’s type, or rank n - 1

16.9.173

maxloc, minloc, findloc

transformational

rank-one array of subscripts, or rank n - 1 with dim

16.9.137, 16.9.143, 16.9.88

size, shape, lbound, ubound, rank, is_contiguous

inquiry

scalar, or rank-one array

16.9.194, 16.9.185, 16.9.119, 16.9.215, 16.9.171, 16.9.115

reshape, spread, pack, unpack

transformational

a new array whose shape is computed

16.9.175, 16.9.197, 16.9.158, 16.9.217

cshift, eoshift, transpose

transformational

an array of the same shape (reversed for transpose)

16.9.68, 16.9.77, 16.9.213

merge

elemental

same shape as the conformable arguments

16.9.139

matmul, dot_product

transformational

matrix, vector or scalar

16.9.134, 16.9.73

Reductions

A reduction collapses an array along one dimension or all of them. The simplest are sum, product, maxval and minval (for numeric and, for the last two, character data), count (how many true elements of a logical array), any and all. They are all called as f(array [, dim] [, mask]), with dim and mask as keywords whenever one is skipped, and count also takes kind. The first example applies each of them to a 4 by 3 temperature field and prints the result, so the effect of dim=1 (collapse the first dimension, one result per column), dim=2 (one result per row) and mask can be read from the numbers. It ends with the results for zero-sized arrays and with count called with a kind.

program reductions
  use, intrinsic :: iso_fortran_env, only: real64, int8, int64
  implicit none
  integer, parameter :: dp = real64
  real(dp) :: u(4, 3)
  logical :: hot(4, 3)
  integer :: i, j

  ! temperatures on a 4 x 3 grid; u(i, j) = 10*i + j makes element order visible
  do j = 1, 3
    do i = 1, 4
      u(i, j) = real(10*i + j, dp)
    end do
  end do
  hot = u > 25.0_dp

  print '(a, f7.1)',  'sum(u)               = ', sum(u)
  print '(a, 3f7.1)', 'sum(u, dim=1)        = ', sum(u, dim=1)
  print '(a, 4f7.1)', 'sum(u, dim=2)        = ', sum(u, dim=2)
  print '(a, f7.1)',  'sum(u, mask=hot)     = ', sum(u, mask=hot)
  print '(a, 3f7.1)', 'product(u(1:2,:),1)  = ', product(u(1:2, :), dim=1)
  print '(a, f7.1)',  'maxval(u)            = ', maxval(u)
  print '(a, 3f7.1)', 'minval(u, dim=1)     = ', minval(u, dim=1)
  print '(a, f7.1)',  'maxval(u, u < 25)    = ', maxval(u, mask=u < 25.0_dp)
  print '(a, i0)',    'count(hot)           = ', count(hot)
  print '(a, 3i2)',   'count(hot, dim=1)    = ', count(hot, dim=1)
  print '(a, l1, a, l1)', 'any(hot) = ', any(hot), '  all(hot) = ', all(hot)
  print '(a, 3l2)',   'all(hot, dim=1)     = ', all(hot, dim=1)
  print '(a, 4l2)',   'any(hot, dim=2)     = ', any(hot, dim=2)

  ! empty arrays: the identity element of each reduction
  print '(a, f4.1, a, f4.1, a, i0, a, l1, a, l1)', 'empty: sum=', sum(u(1:0, 1)), &
    ' product=', product(u(1:0, 1)), ' count=', count(hot(1:0, 1)), ' any=', any(hot(1:0, 1)), &
    ' all=', all(hot(1:0, 1))
  print '(a, es14.4e3, a, es14.4e3)', 'empty: maxval=', maxval(u(1:0, 1)), ' minval=', minval(u(1:0, 1))

  ! kind= on count and the integer overflow trap of small kinds
  print '(a, i0, a, i0)', 'kind(count(hot))=', kind(count(hot)), '  kind(count(hot, kind=int64))=', &
    kind(count(hot, kind=int64))
  block
    integer(int8) :: small(3)
    small = [100_int8, 100_int8, 100_int8]
    print '(a, i0)', 'sum of three int8 values of 100 (int64 widening): ', sum(int(small, int64))
  end block
end program reductions
sum(u)               =   324.0
sum(u, dim=1)        =   104.0  108.0  112.0
sum(u, dim=2)        =    36.0   66.0   96.0  126.0
sum(u, mask=hot)     =   222.0
product(u(1:2,:),1)  =   231.0  264.0  299.0
maxval(u)            =    43.0
minval(u, dim=1)     =    11.0   12.0   13.0
maxval(u, u < 25)    =    23.0
count(hot)           = 6
count(hot, dim=1)    =  2 2 2
any(hot) = T  all(hot) = F
all(hot, dim=1)     =  F F F
any(hot, dim=2)     =  F F T T
empty: sum= 0.0 product= 1.0 count=0 any=F all=T
empty: maxval=  -1.7977E+308 minval=   1.7977E+308
kind(count(hot))=4  kind(count(hot, kind=int64))=8
sum of three int8 values of 100 (int64 widening): 300

Source: J3/23-007r1 (Fortran 2023), 16.9.201 SUM, 16.9.165 PRODUCT, 16.9.138 MAXVAL, 16.9.144 MINVAL, 16.9.12 ALL, 16.9.15 ANY and 16.9.66 COUNT

What the output shows and the rules behind it:

  • sum(u, dim=1) has three elements, one per column of u(4, 3); sum(u, dim=2) has four, one per row. The dimension named in dim is the one that disappears.

  • The mask is applied before reducing: sum(u, mask=hot) adds only the six elements greater than 25, and maxval(u, mask=u < 25) is the maximum of the elements that satisfy the mask (23), not of the whole array.

  • count(hot, dim=1) and all(hot, dim=1) follow the same shape rule as sum: one value per column.

  • Empty arrays return the identity of the operation. sum of nothing is 0, product is 1, count is 0, any is false and all is true. maxval of a zero-sized integer or real array is the negative number of largest magnitude the type supports (-huge), and minval is +huge; for a zero-sized character array they return a string of char(0) (16.9.138). A loop that starts from the first element would instead read outside the array, which is why these functions are safer than handwritten reductions on slices that can be empty (an inactive boundary strip, for example).

  • The result has the type and kind of the argument. sum of a default-integer array is a default integer and can overflow; the standard does not say what happens then, and compilers do not trap it even with -fcheck=all. Convert the array to a wider kind first (sum(int(small, int64))) when sums may be large. For real data the same applies to accuracy: summing a million real32 values in a real32 accumulator loses digits that a real64 one would keep.

  • The order of evaluation is not specified. The standard lets a processor evaluate any mathematically equivalent expression (10.1.5.2.4), so the last bits of a floating-point sum can depend on the optimisation level, on vectorisation and on OpenMP. A test that compares two sums for exact equality is wrong; compare with a tolerance.

  • count returns a default integer unless kind= is given, so count(mask, kind=int64) is the way to count more than huge(0) elements.

maxval and minval have no back argument and return values, not positions; the positions come from maxloc and minloc (below). The result of maxval for an array that contains a NaN is not defined by the standard and can differ between processors, so clean the data first or reduce with reduce and an operation that tests ieee_is_nan (see Floating-Point and IEEE Arithmetic).

Norms, bit reductions and parity

norm2(x [, dim]) is the Euclidean (L2) norm of a real array. It is not simply sqrt(sum(x**2)): the standard recommends that the processor compute it without undue overflow or underflow (16.9.153), and both compilers do. dot_product(x, y) is the inner product of two vectors of the same size. iall, iany and iparity reduce an integer array with the bitwise iand, ior and ieor; parity reduces a logical array with .neqv., so it is true when an odd number of elements are true. All of them accept dim and, except dot_product, mask.

The second example uses them on the heatflow data, adds reduce, and is described in the next subsection. The bit reductions are applied to four small flag words (one bit per kind of boundary condition seen by a cell).

module reduce_ops
  use, intrinsic :: iso_fortran_env, only: real64
  implicit none
  private
  public :: dp, absmax, hotter
  integer, parameter :: dp = real64

  type, public :: sample_t
    real(dp) :: t = 0.0_dp
    integer :: cell = 0
  end type sample_t

contains

  pure function absmax(a, b) result(r)       ! associative and commutative
    real(dp), intent(in) :: a, b
    real(dp) :: r
    r = max(abs(a), abs(b))
  end function absmax

  pure function hotter(a, b) result(r)       ! keeps the sample with the larger t
    type(sample_t), intent(in) :: a, b
    type(sample_t) :: r
    if (b%t > a%t) then
      r = b
    else
      r = a
    end if
  end function hotter

end module reduce_ops

program norms_bits_reduce
  use, intrinsic :: iso_fortran_env, only: int32
  use reduce_ops
  implicit none
  real(dp) :: r(5), g(2, 3), big(2)
  integer(int32) :: flags(4)
  logical :: odd(4)
  type(sample_t) :: s(4), best
  integer :: k

  r = [3.0_dp, -4.0_dp, 0.0_dp, 12.0_dp, -84.0_dp]
  g = reshape([1.0_dp, 3.0_dp, 2.0_dp, 4.0_dp, 5.0_dp, 6.0_dp], [2, 3])
  print '(a, f10.4)', 'norm2(r)          = ', norm2(r)
  print '(a, f10.4)', 'sqrt(sum(r**2))   = ', sqrt(sum(r**2))
  print '(a, 3f9.4)', 'norm2(g, dim=1)   = ', norm2(g, dim=1)
  print '(a, 2f9.4)', 'norm2(g, dim=2)   = ', norm2(g, dim=2)
  big = 1.0e200_dp
  print '(a, es10.3e3)', 'norm2 of [1e200, 1e200]: ', norm2(big)
  print '(a, es10.3e3)', 'sqrt(sum(x**2))        : ', sqrt(sum(big**2))
  print '(a, f10.4)', 'dot_product(r, r) = ', dot_product(r, r)

  ! one bit per boundary condition seen in a cell
  flags = [3_int32, 6_int32, 5_int32, 7_int32]   ! 0011 0110 0101 0111
  print '(a, b4.4, a, b4.4, a, b4.4)', 'iall=', iall(flags), '  iany=', iany(flags), '  iparity=', iparity(flags)
  odd = [(btest(flags(k), 0), k = 1, 4)]
  print '(a, l1, a, l1)', 'parity(odd)=', parity(odd), '  count(odd) odd? ', mod(count(odd), 2) == 1

  ! reduce: user-defined operation
  print '(a, f6.1)', 'reduce(r, absmax)          = ', reduce(r, absmax)
  print '(a, f6.1)', 'reduce(r, absmax, r > 0)   = ', reduce(r, absmax, mask=r > 0.0_dp, identity=0.0_dp)
  print '(a, f6.1)', 'nothing selected, identity 0.0 = ', reduce(r, absmax, mask=r > 100.0_dp, identity=0.0_dp)
  print '(a, 3f6.1)', 'reduce(g, absmax, dim=1)   = ', reduce(g, absmax, dim=1)
  s = [sample_t(1.5_dp, 1), sample_t(9.0_dp, 2), sample_t(4.0_dp, 3), sample_t(9.0_dp, 4)]
  best = reduce(s, hotter, ordered=.true.)
  print '(a, f4.1, a, i0)', 'hottest sample: t=', best%t, ' cell=', best%cell
end program norms_bits_reduce
norm2(r)          =    85.0000
sqrt(sum(r**2))   =    85.0000
norm2(g, dim=1)   =    3.1623   4.4721   7.8102
norm2(g, dim=2)   =    5.4772   7.8102
norm2 of [1e200, 1e200]: 1.414E+200
sqrt(sum(x**2))        :   Infinity
dot_product(r, r) =  7225.0000
iall=0000  iany=0111  iparity=0111
parity(odd)=T  count(odd) odd? T
reduce(r, absmax)          =   84.0
reduce(r, absmax, r > 0)   =   12.0
nothing selected, identity 0.0 =    0.0
reduce(g, absmax, dim=1)   =    3.0   4.0   6.0
hottest sample: t= 9.0 cell=2

Source: J3/23-007r1 (Fortran 2023), 16.9.153 NORM2, 16.9.99 IALL, 16.9.101 IANY, 16.9.112 IPARITY, 16.9.159 PARITY and 16.9.173 REDUCE

From the output: norm2(g, dim=1) is the norm of each column of the 2 by 3 array g, and norm2(g, dim=2) the norm of each row; norm2 of two values of 1e200 is 1.414e200 where the naive sqrt(sum(big**2)) overflows to infinity; iall of 3, 6, 5 and 7 is 0 (the three-way iand leaves no common bit), iany is 7, iparity is the xor 7, and the parity of the four low bits is true because three of the four flags are odd. Bit manipulation in general, including the bit model of 16.3, is on Numeric Model, Bit Manipulation and Conversion.

reduce: a reduction with your own operation

reduce (Fortran 2018) is the general reduction. reduce(array, operation [, mask, identity, ordered]) or reduce(array, operation, dim [, mask, identity, ordered]) applies a pure binary function repeatedly to adjacent elements until one value is left. It works for arrays of any type, including derived types, which none of the other reductions do. The rules, from J3/23-007r1, 16.9.173:

  • operation is a pure function of exactly two arguments, both scalar, non-allocatable, non-pointer, non-polymorphic, and both of the same declared type and type parameters as array; the result has that type as well. A function that does not fit this description does not conform.

  • operation should be associative, because the processor may group the elements in any way, and is not required to be commutative. If it is not associative, or if it depends on the order, give ordered=.true., which makes the reduction combine the first two elements of the remaining sequence at each step, in array element order. The price is that the compiler can no longer evaluate in parallel.

  • mask selects the elements, as for the other reductions, and dim reduces one dimension.

  • If the selected sequence is empty, the result is identity when it is present; otherwise the program is error terminated. That is why the identity argument matters whenever a mask can select nothing.

  • Fortran 2023 adds the reduce locality specifier to do concurrent for loops written as reductions; that is a separate feature (see do concurrent and GPU Offload), and the coarray collective for the same job is co_reduce (see Coarray Teams, Events, Collectives and Atomics).

In the example, absmax finds the element of largest magnitude, reduce(r, absmax, mask=r > 0.0_dp, identity=0.0_dp) restricts it to positive elements (12), the masked call that selects nothing returns the identity, reduce(g, absmax, dim=1) reduces columns, and hotter reduces an array of a derived type sample_t (temperature and cell number) to the sample with the highest temperature, which sum or maxval could not do.

The compilers differ here, and the difference is a bug rather than an extension, so it is worth showing in a program of its own. With GNU Fortran 16.2, reduce of a zero-sized array without a mask returns an undefined value instead of the identity, whereas the masked form (a mask that selects nothing) works and LLVM Flang 23 returns the identity in all three cases. The program is conforming; the first two lines should print 7.

program reduce_empty
  use, intrinsic :: iso_fortran_env, only: real64
  implicit none
  integer, parameter :: dp = real64
  real(dp) :: r(5), none(0)
  r = [3.0_dp, -4.0_dp, 0.0_dp, 12.0_dp, -84.0_dp]
  print '(a, g0)', 'zero-sized array, identity=7: ', reduce(none, absmax, identity=7.0_dp)
  print '(a, g0)', 'zero-sized section, identity=7: ', reduce(r(1:0), absmax, identity=7.0_dp)
  print '(a, g0)', 'nothing selected by mask, identity=7: ', reduce(r, absmax, mask=r > 100.0_dp, identity=7.0_dp)
contains
  pure function absmax(a, b) result(c)
    real(dp), intent(in) :: a, b
    real(dp) :: c
    c = max(abs(a), abs(b))
  end function absmax
end program reduce_empty
$ gfortran -std=f2023 -Wall -Wextra -pedantic reduce_empty.f90 && ./a.out    (GNU Fortran 16.2)
zero-sized array, identity=7: 0.13906697895580856E-308
zero-sized section, identity=7: 3.0000000000000000
nothing selected by mask, identity=7: 7.0000000000000000

$ flang reduce_empty.f90 && ./a.out    (LLVM Flang 23)
zero-sized array, identity=7: 7.
zero-sized section, identity=7: 7.
nothing selected by mask, identity=7: 7.

Source: J3/23-007r1 (Fortran 2023), 16.9.173 REDUCE

The first two values on the GNU Fortran side change from build to build (it printed 0.13262473693532952E-314 at -O2) and, in the second case, are the first element of the array that should have been ignored. Until the compiler is fixed, guard the call (if (size(x) > 0)) or pass a mask. GNU Fortran’s linker step also printed "requires executable stack" for programs that call reduce; that is a linker warning about the generated code, not a compiler warning, and it does not appear with Flang. GNU Fortran’s manual lists the other intrinsics of this page individually but has no page for reduce at the time of writing.

Location functions

maxloc and minloc return the position of the largest or smallest element; findloc returns the position of an element equal to a given value (Fortran 2008). All three return subscripts, not values and not offsets, and their result form is easy to get wrong:

  • Without dim, the result is a rank-one integer array with one element per dimension of the argument, even for a rank-one argument (an array of size 1). maxloc(u) is [i, j] for a rank-2 u; assigning it to a scalar is an error. Declare an integer array of size rank(u) or use maxloc(u, dim=1) for rank-one arrays.

  • With dim, the result is the position along that dimension for every vector of the other dimensions, with rank n - 1. maxloc(u, dim=1) has one element per column. For a rank-one argument the result with dim is a scalar.

  • Ties return the first match in array element order, and with back=.true. the last (Fortran 2008). The back argument makes maxloc(v) and maxloc(v, back=.true.) different when the maximum occurs more than once.

  • No match returns zero: if the array has size zero, if the mask selects nothing, or if findloc finds no element equal to value, every element of the result is 0 (16.9.137, 16.9.88). Zero is never a valid subscript for a lower bound of 1, but it can be a valid subscript for an array declared with 0:, so test with any(loc == 0) only for arrays whose lower bounds are 1.

  • The subscripts are relative to the array as if its lower bounds were 1, not to its declared bounds. For an array declared w(0:3, -1:1), maxloc(w) still counts from one; to get the subscript to use with w, add the lower bound and subtract one, as in loc + lbound(w) - 1. This is the usual source of off-by-one errors with these functions.

  • maxloc and minloc take integer, real and character arrays; findloc takes arrays of any intrinsic type (numeric, logical or character) and a scalar value of a type that can be compared. Because it tests equality, findloc on real data is only as exact as == is; use it with integers, logicals and characters, or with reals that are known to be bit-equal. To locate the first element satisfying a condition, use findloc(condition, .true., dim=1).

  • kind= sets the integer kind of the result, as for count.

The example uses a 4 by 3 field in which the maximum 43 occurs twice (at u(2,3) and u(4,3)), and a rank-one array with a repeated maximum.

program locations
  use, intrinsic :: iso_fortran_env, only: real64, int64
  implicit none
  integer, parameter :: dp = real64
  real(dp) :: u(4, 3)
  real(dp), allocatable :: w(:, :)
  real(dp) :: v(6)
  integer :: i, j, loc(2)
  integer(int64) :: loc64(2)

  do j = 1, 3
    do i = 1, 4
      u(i, j) = real(10*i + j, dp)
    end do
  end do
  u(2, 3) = u(4, 3)                       ! the maximum 43.0 now appears twice: (2,3) and (4,3)
  v = [2.0_dp, 6.0_dp, 4.0_dp, 6.0_dp, 1.0_dp, 6.0_dp]

  loc = maxloc(u)                         ! rank-one array of size 2, not a scalar
  print '(a, 2i3)', 'maxloc(u)            = ', loc
  print '(a, 2i3)', 'maxloc(u, back=.true.) = ', maxloc(u, back=.true.)
  print '(a, 2i3)', 'minloc(u)            = ', minloc(u)
  print '(a, 3i3)', 'maxloc(u, dim=1)     = ', maxloc(u, dim=1)
  print '(a, 4i3)', 'minloc(u, dim=2)     = ', minloc(u, dim=2)
  print '(a, 2i3)', 'maxloc(u, mask=u<40) = ', maxloc(u, mask=u < 40.0_dp)
  print '(a, i0)',  'maxloc(v, dim=1)     = ', maxloc(v, dim=1)
  print '(a, 3i2)', 'maxloc(v), maxloc(v, back=.true.), size(maxloc(v)) = ', maxloc(v), maxloc(v, back=.true.), size(maxloc(v))
  loc64 = maxloc(u, kind=int64)
  print '(a, i0)',  'kind of maxloc(u, kind=int64): ', kind(loc64)
  print '(a, 2i3)', 'maxloc of nothing selected  = ', maxloc(u, mask=u < 0.0_dp)

  ! findloc: where is a given value (equality, so use it for integers, logicals, characters)
  print '(a, 2i3)', 'findloc(u, 43.0)       = ', findloc(u, 43.0_dp)
  print '(a, 2i3)', 'findloc(u, 43.0, back) = ', findloc(u, 43.0_dp, back=.true.)
  print '(a, 3i3)', 'findloc(u, 43.0, dim=1)= ', findloc(u, 43.0_dp, dim=1)
  print '(a, 2i3)', 'findloc(u, 99.0)       = ', findloc(u, 99.0_dp)
  print '(a, i0, a, i0)', 'first .true. in v>5: ', findloc(v > 5.0_dp, .true., dim=1), &
    '  last: ', findloc(v > 5.0_dp, .true., dim=1, back=.true.)

  ! subscripts come back relative to the array as an expression, not the declared bounds
  allocate (w(0:3, -1:1))
  w = u
  loc = maxloc(w)
  print '(a, 2i3, a, 2i3)', 'maxloc(w) with w(0:3,-1:1) = ', loc, '   declared subscripts of that element: ', &
    loc + lbound(w) - 1
  print '(a, f5.1, a, f5.1)', 'w(1, 1) = ', w(1, 1), '   maxval(w) = ', maxval(w)
end program locations
maxloc(u)            =   2  3
maxloc(u, back=.true.) =   4  3
minloc(u)            =   1  1
maxloc(u, dim=1)     =   4  4  2
minloc(u, dim=2)     =   1  1  1  1
maxloc(u, mask=u<40) =   3  3
maxloc(v, dim=1)     = 2
maxloc(v), maxloc(v, back=.true.), size(maxloc(v)) =  2 6 1
kind of maxloc(u, kind=int64): 8
maxloc of nothing selected  =   0  0
findloc(u, 43.0)       =   2  3
findloc(u, 43.0, back) =   4  3
findloc(u, 43.0, dim=1)=   0  0  2
findloc(u, 99.0)       =   0  0
first .true. in v>5: 2  last: 6
maxloc(w) with w(0:3,-1:1) =   2  3   declared subscripts of that element:   1  1
w(1, 1) =  43.0   maxval(w) =  43.0

Source: J3/23-007r1 (Fortran 2023), 16.9.137 MAXLOC, 16.9.143 MINLOC and 16.9.88 FINDLOC

The output confirms the rules: maxloc(u) is 2 3 (the first 43 in element order) and maxloc(u, back=.true.) is 4 3; maxloc(u, dim=1) gives 4 4 2 (the first row that holds each column’s maximum); maxloc(v) has size 1; findloc(u, 43.0, dim=1) gives 0 0 2 because columns 1 and 2 contain no 43; and for w(0:3, -1:1) the element maxloc reports as 2 3 is w(1, 1).

Inquiry functions

The inquiry functions describe an array without computing anything from its values, which is why the standard lets them be applied to arrays whose elements are undefined (16.2.2).

Function Result

size(array [, dim, kind])

The number of elements, or the extent along dim. The result is a default integer unless kind= is given.

shape(source [, kind])

A rank-one array of the extents; for a scalar it is an array of size zero.

lbound(array [, dim, kind]), ubound(array [, dim, kind])

The lower and upper bounds. Without dim the result is a rank-one array of size rank(array).

rank(a)

The rank of any data object as a scalar integer (Fortran 2018); 0 for a scalar. Works on assumed-rank dummies, where it is the way to find the rank at run time.

is_contiguous(array)

True if the array is contiguous in memory (Fortran 2008). For a pointer it asks whether its target is contiguous.

Details that matter in practice (16.9.119, 16.9.215, 16.9.185, 16.9.194):

  • lbound and ubound of a whole array return the declared bounds, so for u(0:4, -1:2) they are [0, -1] and [4, 2]. Anything that is not a whole array has lower bounds of 1: an array section, an array expression and a function result. An assumed-shape dummy argument also has lower bounds of 1 unless the dummy declares another lower bound, whatever the bounds of the actual argument were.

  • For a dimension of extent zero, lbound is 1 and ubound is 0, so ubound(a, dim) - lbound(a, dim) + 1 is always the extent (zero), and ubound is smaller than lbound.

  • size(a) is the product of shape(a) and is zero when any extent is zero.

  • shape and size take an optional kind, but rank does not.

  • size(shape(x)) is the rank of x, and is 0 for a scalar.

module inquiry_demo
  use, intrinsic :: iso_fortran_env, only: real64
  implicit none
  integer, parameter :: dp = real64
contains
  subroutine describe(a)
    real(dp), intent(in) :: a(:, :)           ! assumed shape: lower bounds are 1
    print '(a, 2i3, a, 2i3, a, 2i3)', '  in describe: shape', shape(a), '  lbound', lbound(a), '  ubound', ubound(a)
  end subroutine describe
end module inquiry_demo

program inquiry
  use, intrinsic :: iso_fortran_env, only: int8, int64
  use inquiry_demo
  implicit none
  real(dp), allocatable :: u(:, :)
  real(dp) :: tile(2, 2)
  real(dp), pointer :: p(:, :)
  real(dp), allocatable, target :: store(:, :)
  integer :: n = 0

  allocate (u(0:4, -1:2), source=0.0_dp)           ! 5 x 4, declared bounds 0:4 and -1:2
  print '(a, 2i3)', 'shape(u)          =', shape(u)
  print '(a, 2i3)', 'lbound(u)         =', lbound(u)
  print '(a, 2i3)', 'ubound(u)         =', ubound(u)
  print '(a, 3i3)', 'size(u), size(u,1), size(u,2) =', size(u), size(u, 1), size(u, 2)
  print '(a, i0)',  'rank(u)           = ', rank(u)
  print '(a, i0)',  'rank(3.0)         = ', rank(3.0)
  print '(a, 2i3)', 'lbound(u, dim=2), ubound(u, dim=2) =', lbound(u, dim=2), ubound(u, dim=2)
  print '(a, i0, a, i0)', 'kind(shape(u)) = ', kind(shape(u)), '   kind(shape(u, kind=int8)) = ', kind(shape(u, kind=int8))
  print '(a, i0, a, i0)', 'kind(size(u, kind=int64)) = ', kind(size(u, kind=int64)), '   kind(size(u)) = ', kind(size(u))

  print '(a)', 'the same array passed to an assumed-shape dummy:'
  call describe(u)
  print '(a)', 'an array section and an array expression:'
  call describe(u(1:3, 0:1))
  call describe(u + 1.0_dp)
  print '(a, 2i3)', '  lbound of the section u(1:3,0:1) in place =', lbound(u(1:3, 0:1))
  print '(a, 2i3)', '  lbound of the expression (u + 1)          =', lbound(u + 1.0_dp)

  ! zero-sized dimensions
  deallocate (u)
  allocate (u(3, 0))
  print '(a, 2i3, a, 2i3, a, i0)', 'zero size: shape', shape(u), '  ubound', ubound(u), '  size ', size(u)

  ! scalars, and inquiry without a defined value
  print '(a, i0, a, i0)', 'size(shape(tile)) = ', size(shape(tile)), '   size(shape(n)) = ', size(shape(n))

  ! pointer association and contiguity
  allocate (store(6, 6), source=0.0_dp)
  p => store(::2, ::2)
  print '(a, l1, a, l1, a, l1)', 'is_contiguous(store) ', is_contiguous(store), '   (store(::2,::2)) ', &
    is_contiguous(p), '   (store(:,2:3)) ', is_contiguous(store(:, 2:3))
  print '(a, 2i3)', 'shape(p) = ', shape(p)
end program inquiry
shape(u)          =  5  4
lbound(u)         =  0 -1
ubound(u)         =  4  2
size(u), size(u,1), size(u,2) = 20  5  4
rank(u)           = 2
rank(3.0)         = 0
lbound(u, dim=2), ubound(u, dim=2) = -1  2
kind(shape(u)) = 4   kind(shape(u, kind=int8)) = 1
kind(size(u, kind=int64)) = 8   kind(size(u)) = 4
the same array passed to an assumed-shape dummy:
  in describe: shape  5  4  lbound  1  1  ubound  5  4
an array section and an array expression:
  in describe: shape  3  2  lbound  1  1  ubound  3  2
  in describe: shape  5  4  lbound  1  1  ubound  5  4
  lbound of the section u(1:3,0:1) in place =  1  1
  lbound of the expression (u + 1)          =  1  1
zero size: shape  3  0  ubound  3  0  size 0
size(shape(tile)) = 2   size(shape(n)) = 0
is_contiguous(store) T   (store(::2,::2)) F   (store(:,2:3)) T
shape(p) =   3  3

Source: J3/23-007r1 (Fortran 2023), 16.9.194 SIZE, 16.9.185 SHAPE, 16.9.119 LBOUND, 16.9.215 UBOUND, 16.9.171 RANK and 16.9.115 IS_CONTIGUOUS

The output shows the declared bounds 0 -1 and 4 2 of u(0:4, -1:2), the dummy argument in describe seeing lower bounds 1 1 for the whole array, a section and an expression alike, and the zero-sized u(3, 0) with ubound equal to 3 0. The pointer p => store(::2, ::2) is not contiguous, while store(:, 2:3), a block of whole columns, is. The consequences for argument passing, copy-in and copy-out, and the contiguous attribute are on Array Arguments: Assumed Shape, Size and Rank.

Construction and manipulation

These functions build a new array from existing ones. None of them changes its arguments.

Call What it returns

reshape(source, shape [, pad, order])

The elements of source, in array element order, as an array of shape shape. If source has too few elements, pad (an array of the same type) is used cyclically to fill the rest; without pad it is an error. order is a permutation of 1..rank(shape) that lists the dimensions from the one that varies fastest to the slowest, so order=[2, 1] fills by rows. shape is a rank-one integer array of constant size whose elements are not negative.

spread(source, dim, ncopies)

source replicated ncopies times along a new dimension dim; the result has rank rank(source) + 1. A scalar source gives a rank-one array, which makes spread the explicit form of broadcasting a scalar or a row to a matrix.

pack(array, mask [, vector])

The elements of array selected by mask, in array element order, as a rank-one array. If vector is present, the result has the size of vector: the selected elements come first and the rest of vector follows. vector must have at least as many elements as the mask selects.

unpack(vector, mask, field)

An array of the shape of mask in which the true positions receive consecutive elements of vector and the others the corresponding element of field (an array or a scalar). The inverse of pack. vector must have at least as many elements as mask has true positions.

merge(tsource, fsource, mask)

Elemental: element-wise, tsource where mask is true and fsource where it is false. Unlike an if, both sources are evaluated, so the unselected one must still be a valid value (a division by zero in it is still a division by zero).

cshift(array, shift [, dim])

A circular shift: element i of the result is array(i + shift), wrapping around. A positive shift moves the elements toward lower subscripts. For a rank-two array, shift may be an array that gives a different amount to each row or column, and dim defaults to 1.

eoshift(array, shift [, boundary, dim])

An end-off shift: the elements that fall off one end are lost and boundary (default 0, blank for characters, false for logicals) fills the vacated positions. boundary can be an array of one rank lower for rank-two array.

transpose(matrix)

The rank-two array with rows and columns exchanged. It has no dim; for higher ranks use reshape with order.

Source: J3/23-007r1 (Fortran 2023), 16.9.175 RESHAPE, 16.9.197 SPREAD, 16.9.158 PACK, 16.9.217 UNPACK, 16.9.139 MERGE, 16.9.68 CSHIFT, 16.9.77 EOSHIFT and 16.9.213 TRANSPOSE

The example runs each of them on small integer arrays so the result can be checked by hand. The two lines at the end build a one-dimensional second difference (the 1D heat Laplacian) twice: with cshift, which treats the array as periodic, and with eoshift, which keeps the edges fixed by supplying the boundary value.

program construction
  use, intrinsic :: iso_fortran_env, only: real64
  implicit none
  integer, parameter :: dp = real64
  integer :: a(2, 3), b(3, 2), m(3, 3)
  integer :: v(6) = [1, 2, 3, 4, 5, 6]
  logical :: keep(2, 3)
  integer, allocatable :: packed(:)
  real(dp) :: u(5), lap(5)

  ! reshape: array element order is column-major; order= permutes it; pad= fills the rest
  a = reshape(v, [2, 3])
  print '(a, 6i3)', 'reshape(v,[2,3])            ', a
  b = reshape(v, [3, 2], order=[2, 1])
  print '(a, 6i3)', 'reshape(v,[3,2],order=[2,1])', b
  m = reshape([1, 2, 3, 4], [3, 3], pad=[0], order=[2, 1])
  print '(a, 9i3)', 'with pad=[0]                ', m
  print '(a, 6i3)', 'transpose(a) = ', transpose(a)

  ! spread: replicate along a new dimension
  print '(a, 6i3)', 'spread([1,2,3],dim=1,ncopies=2) ', spread([1, 2, 3], dim=1, ncopies=2)
  print '(a, 6i3)', 'spread([1,2,3],dim=2,ncopies=2) ', spread([1, 2, 3], dim=2, ncopies=2)

  ! pack and unpack
  keep = a > 2
  packed = pack(a, keep)
  print '(a, 4i3)', 'pack(a, a>2)                ', packed
  print '(a, 6i3)', 'pack(a, a>2, vector=v)      ', pack(a, keep, vector=v)
  print '(a, 6i3)', 'unpack(packed, keep, 0)     ', unpack(packed, keep, 0)
  print '(a, 6i3)', 'unpack(packed, keep, field=-a)', unpack(packed, keep, field=-a)

  ! merge: elementwise choice
  print '(a, 6i3)', 'merge(a, -a, keep)          ', merge(a, -a, keep)
  print '(a, 6i3)', 'merge(a, 0, mod(a,2)==0)    ', merge(a, 0, mod(a, 2) == 0)

  ! cshift and eoshift: periodic and fixed boundaries of a 1D stencil
  print '(a, 6i3)', 'cshift(v, 1)                ', cshift(v, 1)
  print '(a, 6i3)', 'cshift(v, -2)               ', cshift(v, -2)
  print '(a, 6i3)', 'eoshift(v, 2)               ', eoshift(v, 2)
  print '(a, 6i3)', 'eoshift(v, -2, boundary=-1) ', eoshift(v, -2, boundary=-1)
  print '(a, 6i3)', 'cshift(a, shift=[1,-1], dim=2)', cshift(a, shift=[1, -1], dim=2)
  print '(a, 6i3)', 'eoshift(a, 1, dim=1)        ', eoshift(a, 1, dim=1)

  u = [0.0_dp, 1.0_dp, 4.0_dp, 9.0_dp, 16.0_dp]
  lap = cshift(u, 1) - 2.0_dp*u + cshift(u, -1)               ! periodic second difference
  print '(a, 5f6.1)', 'periodic  laplacian  ', lap
  lap = eoshift(u, 1, boundary=u(5)) - 2.0_dp*u + eoshift(u, -1, boundary=u(1))   ! fixed edges
  print '(a, 5f6.1)', 'fixed     laplacian  ', lap
end program construction
reshape(v,[2,3])              1  2  3  4  5  6
reshape(v,[3,2],order=[2,1])  1  3  5  2  4  6
with pad=[0]                  1  4  0  2  0  0  3  0  0
transpose(a) =   1  3  5  2  4  6
spread([1,2,3],dim=1,ncopies=2)   1  1  2  2  3  3
spread([1,2,3],dim=2,ncopies=2)   1  2  3  1  2  3
pack(a, a>2)                  3  4  5  6
pack(a, a>2, vector=v)        3  4  5  6  5  6
unpack(packed, keep, 0)       0  0  3  4  5  6
unpack(packed, keep, field=-a) -1 -2  3  4  5  6
merge(a, -a, keep)           -1 -2  3  4  5  6
merge(a, 0, mod(a,2)==0)      0  2  0  4  0  6
cshift(v, 1)                  2  3  4  5  6  1
cshift(v, -2)                 5  6  1  2  3  4
eoshift(v, 2)                 3  4  5  6  0  0
eoshift(v, -2, boundary=-1)  -1 -1  1  2  3  4
cshift(a, shift=[1,-1], dim=2)  3  6  5  2  1  4
eoshift(a, 1, dim=1)          2  0  4  0  6  0
periodic  laplacian    17.0   2.0   2.0   2.0 -23.0
fixed     laplacian     1.0   2.0   2.0   2.0  -7.0

Source: J3/23-007r1 (Fortran 2023), 16.9.175 RESHAPE, 16.9.197 SPREAD, 16.9.158 PACK, 16.9.217 UNPACK, 16.9.139 MERGE, 16.9.68 CSHIFT, 16.9.77 EOSHIFT and 16.9.213 TRANSPOSE

The cases to read carefully: reshape(v, [3, 2], order=[2, 1]) returns 1 3 5 2 4 6 in memory order, which is the matrix [[1, 2], [3, 4], [5, 6]] filled row by row; spread([1,2,3], dim=1, ncopies=2) gives a 2 by 3 array whose rows are the vector (printed column by column); pack(a, a>2, vector=v) keeps the selected 3 4 5 6 and then the unused tail 5 6 of v; unpack with a scalar 0 as field zeroes the elements that were not selected; cshift(v, 1) moves 1 from the front to the back and cshift(v, -2) the last two to the front; eoshift(v, 2) drops the first two and appends two zeros; and cshift(a, shift=[1, -1], dim=2) shifts the first row one way and the second the other. In the Laplacian lines, the periodic version wraps the end values, giving 17 and -23 at the edges where the eoshift version, with the end values repeated as the boundary, gives 1 and -7.

Shifts, sections and array temporaries

cshift and eoshift are convenient and they are also the functions where the convenience has a cost. Each call creates a complete shifted copy of the array, so a stencil written with four shifts allocates and fills four temporaries per sweep. Writing the same update with array sections needs no temporaries: the compiler reads the original array through the shifted views. The next program computes the five-point Laplacian of a 6 by 5 field both ways; the interior values are identical and the difference is what the compiler has to do.

program stencil
  use, intrinsic :: iso_fortran_env, only: real64
  implicit none
  integer, parameter :: dp = real64
  integer, parameter :: nx = 6, ny = 5
  real(dp) :: u(nx, ny), lap_shift(nx, ny), lap_slice(nx, ny)
  integer :: i, j

  do j = 1, ny
    do i = 1, nx
      u(i, j) = real((i - 3)**2 + (j - 2)**2, dp)
    end do
  end do

  ! (a) shifts: simple to write, but each cshift/eoshift is a whole-array temporary
  lap_shift = cshift(u, 1, dim=1) + cshift(u, -1, dim=1) + cshift(u, 1, dim=2) + cshift(u, -1, dim=2) &
              - 4.0_dp*u

  ! (b) sections: no temporaries, only the interior is updated
  lap_slice = 0.0_dp
  lap_slice(2:nx - 1, 2:ny - 1) = u(3:nx, 2:ny - 1) + u(1:nx - 2, 2:ny - 1) &
                                  + u(2:nx - 1, 3:ny) + u(2:nx - 1, 1:ny - 2) &
                                  - 4.0_dp*u(2:nx - 1, 2:ny - 1)

  print '(a, f6.1)', 'largest interior difference: ', &
    maxval(abs(lap_shift(2:nx - 1, 2:ny - 1) - lap_slice(2:nx - 1, 2:ny - 1)))
  print '(a, f6.1)', 'interior laplacian (all equal for this field): ', lap_slice(3, 3)
  print '(a, 2f6.1)', 'wrapped edge value from cshift vs untouched slice: ', lap_shift(1, 1), lap_slice(1, 1)
end program stencil
largest interior difference:    0.0
interior laplacian (all equal for this field):    4.0
wrapped edge value from cshift vs untouched slice:    9.0   0.0

gfortran -Warray-temporaries makes the cost visible. For this program it reports, for the line that calls cshift four times:

Warning: Creating array temporary at (1) [-Warray-temporaries]
Warning: Creating array temporary at (1) [-Warray-temporaries]
Warning: Creating array temporary at (1) [-Warray-temporaries]
Warning: Creating array temporary at (1) [-Warray-temporaries]

Source: J3/23-007r1 (Fortran 2023), 16.9.68 CSHIFT

and nothing for the statement written with sections. Use shifts for clarity or when the boundary condition is the point (periodic wrap-around is a one-liner with cshift); use sections in a kernel that runs every time step. The same trade-off is discussed for the full solver on Performance and Memory Layout and the sections themselves on Arrays Fundamentals. Other operations that produce temporaries are pack, spread, reshape, transpose and any array expression passed to an argument that needs contiguous storage.

matmul and dot_product

matmul(a, b) multiplies matrices or vectors. The shapes allowed are the ones that make mathematical sense (16.9.134): [n, m] times [m, k] gives [n, k]; [m] times [m, k] gives [k] (a row vector times a matrix); [n, m] times [m] gives [n] (a matrix times a column vector). The two arguments cannot both be rank one; use dot_product or sum(x * y) for that. The type of the result follows the rules of the * operator (so real times complex is complex), and logical arguments give a logical result in which element (i, j) is true when any(a(i, :) .and. b(:, j)).

dot_product(x, y) is the sum of the products of corresponding elements of two vectors of the same size. For complex vectors the first argument is conjugated, dot_product(z, w) is sum(conjg(z) * w), unlike sum(z * w). For logical vectors it is any(x .and. y). The third example below shows all of it.

program matrix
  use, intrinsic :: iso_fortran_env, only: real64
  implicit none
  integer, parameter :: dp = real64
  real(dp) :: a(2, 3), x(3), y(2), b(3, 2), c(2, 2)
  complex(dp) :: z(2), w(2)
  logical :: p(2, 2), q(2)

  a = reshape([1.0_dp, 4.0_dp, 2.0_dp, 5.0_dp, 3.0_dp, 6.0_dp], [2, 3])   ! rows [1 2 3] and [4 5 6]
  x = [1.0_dp, 0.0_dp, -1.0_dp]
  y = [1.0_dp, 2.0_dp]
  b = transpose(a)

  print '(a, 2f7.1)', 'matmul(a, x)      = ', matmul(a, x)
  print '(a, 3f7.1)', 'matmul(y, a)      = ', matmul(y, a)
  c = matmul(a, b)
  print '(a, 4f7.1)', 'matmul(a, b)      = ', c
  print '(a, f7.1)',  'dot_product(x, x) = ', dot_product(x, x)
  print '(a, f7.1)',  'sum(a * a) = ', sum(a*a)

  z = [(1.0_dp, 2.0_dp), (0.0_dp, 1.0_dp)]
  w = [(3.0_dp, 0.0_dp), (0.0_dp, 1.0_dp)]
  print '(a, 2f6.1)', 'dot_product(z, w) = ', dot_product(z, w)
  print '(a, 2f6.1)', 'sum(z * w)        = ', sum(z*w)
  print '(a, 2f6.1)', 'sum(conjg(z) * w) = ', sum(conjg(z)*w)

  p = reshape([.true., .false., .true., .true.], [2, 2])
  q = [.true., .true.]
  print '(a, 2l2)', 'matmul(p, q) (logical: or of ands) = ', matmul(p, q)
  print '(a, l2)',  'dot_product(q, q) (logical)        = ', dot_product(q, q)
end program matrix
matmul(a, x)      =    -2.0   -2.0
matmul(y, a)      =     9.0   12.0   15.0
matmul(a, b)      =    14.0   32.0   32.0   77.0
dot_product(x, x) =     2.0
sum(a * a) =    91.0
dot_product(z, w) =    4.0  -6.0
sum(z * w)        =    2.0   6.0
sum(conjg(z) * w) =    4.0  -6.0
matmul(p, q) (logical: or of ands) =  T T
dot_product(q, q) (logical)        =  T

Source: J3/23-007r1 (Fortran 2023), 16.9.134 MATMUL and 16.9.73 DOT_PRODUCT

In the output, matmul(a, x) with x = [1, 0, -1] gives -2 -2 (the first element of each row minus its last), matmul(y, a) gives 9 12 15, matmul(a, b) with b = transpose(a) is the symmetric 2 by 2 matrix printed in memory order, 14 32 32 77, and for complex input dot_product(z, w) equals sum(conjg(z) * w) (4, -6) but not the unconjugated sum(z * w) (2, 6).

How the compilers implement matmul is a performance decision, not a language one. GNU Fortran inlines small products (-finline-matmul-limit=n, default 30, with front-end optimisation on) and can call an external BLAS gemm for larger ones when asked with -fexternal-blas and a link against a BLAS library (-fblas-matmul-limit=n sets the size above which BLAS is used, and -fexternal-blas64 selects a library with 64-bit integers). Without those options the library routine in libgfortran is used, which is correct but not tuned for large matrices. For big problems, call BLAS and LAPACK directly; see BLAS, LAPACK and Numerical Libraries. Compiler options are covered on Compilers and Compiler Options.

Source: GNU Fortran manual: Options for code generation conventions (-fexternal-blas, -fblas-matmul-limit, -finline-matmul-limit)

Choosing the right function

flowchart TD Q[What do you need?] --> A{A single value or one value per row or column?} A -->|a value| B{What kind of value?} A -->|a position| L[maxloc, minloc or findloc, remember the lower bounds] A -->|a new array| C{How is it built?} B -->|arithmetic| R1[sum, product, norm2, dot_product, matmul] B -->|order| R2[maxval, minval] B -->|logical| R3[count, any, all, parity] B -->|bits| R4[iall, iany, iparity] B -->|anything else| R5[reduce with your own pure function] C -->|same elements, new shape| S1[reshape, transpose, spread] C -->|chosen or moved elements| S2[pack, unpack, merge, cshift, eoshift] C -->|just a section| S3[array sections, no function needed]

Diagnostics

Misuse of these functions is caught at compile time when the shapes or arguments are constants and at run time, with -fcheck=bounds (part of -fcheck=all) or by the library, otherwise. The following non-conforming programs, each of which is in its own file, show the messages. All of them are shown with GNU Fortran 16.2 and then LLVM Flang 23.

Non-conforming statement GNU Fortran 16.2 LLVM Flang 23

c = matmul(a, b) with a(2,3) and b(2,3) (constants)

Error: Different shape on dimension 2 for argument 'matrix_a' and dimension 1 for argument 'matrix_b' at (1) for intrinsic matmul

compiles; at run time fatal Fortran runtime error(...): MATMUL: unacceptable operand shapes (2x3, 2x3)

dot_product(x, y) with sizes 3 and 4 (constants)

Error: Different shape for arguments 'vector_a' and 'vector_b' at (1) for intrinsic 'dot_product'

at run time DOT_PRODUCT: SIZE(VECTOR_A) is 3 but SIZE(VECTOR_B) is 4

sum(a, dim=3) for a(2,3) (constant dim)

Error: 'dim' argument of 'sum' intrinsic at (1) is not a valid dimension index

error: The value of DIM= (3) may not be greater than 2

sum(a, dim=d) with d = 3 (variable dim)

Fortran runtime error: Dim argument incorrect in SUM intrinsic: is 3, should be between 1 and 2

fatal Fortran runtime error(...): SUM: bad DIM=3 for ARRAY with rank 2

reshape(v, [2, 3]) with v(5), no pad

Error: Without padding, there are not enough elements in the intrinsic RESHAPE source at (1) to match the shape

compiles; at run time RESHAPE: not enough elements, need 6 but only have 5

a + b with a(2,3), b(3,2)

Error: Shapes for operands at (1) and (2) are not conformable

error: Dimension 1 of left operand has extent 2, but right operand has extent 3

matmul of two allocatable 2 by 3 arrays (run time)

Fortran runtime error: Incorrect extent in argument B in MATMUL intrinsic in dimension 1: is 2, should be 3

fatal Fortran runtime error(...): MATMUL: unacceptable operand shapes (2x3, 2x3)

dot_product of allocatable vectors of sizes 3 and 4 (run time)

Fortran runtime error: Array bound mismatch for dimension 1 of array 'x' (3/4) (the bounds check on the inlined loop)

fatal Fortran runtime error(...): DOT_PRODUCT: SIZE(VECTOR_A) is 3 but SIZE(VECTOR_B) is 4

Not everything is diagnosed. Conformance of shapes is checked, but a maxloc result assigned to a scalar is rejected only because the types do not match, a reduce over an empty array without identity is an error termination whose message is compiler specific, a wrong lower bound applied to a maxloc result is not diagnosed by anything, and an integer sum that overflows is silent. Run development builds with -fcheck=all (GNU Fortran) and keep -Wall -Wextra enabled.

Where each function comes from

The array intrinsics of Fortran 90 (sum, product, maxval, minval, count, any, all, maxloc, minloc, size, shape, lbound, ubound, reshape, spread, pack, unpack, merge, cshift, eoshift, transpose, matmul, dot_product) are the same in every later standard, with the following changes:

Standard Addition Compiler support

Fortran 2003

kind= argument on count, maxloc, minloc, shape, size, lbound and ubound

Both compilers, all current versions.

Fortran 2008

findloc, norm2, parity, iall, iany, iparity, is_contiguous; the back argument of maxloc and minloc

Both compilers. Used here with GNU Fortran 16.2 and LLVM Flang 23.

Fortran 2018

reduce, rank; assumed-rank dummies accepted by size, shape, lbound, ubound, rank and is_contiguous

Both compilers; GNU Fortran 16.2 has the empty-array defect described above.

Fortran 2023

No new array intrinsic function. reduce locality for do concurrent; @ multiple subscripts and the rank clause for declarations act on arrays but are not intrinsics (see Generic and Rank Independent Programming)

See that page.

Nothing described above is obsolescent or deleted in Fortran 2023. This is a point where older books are out of date: texts written before Fortran 2008 describe only maxloc and minloc for locating elements (and compensate with maxval followed by a masked minval or a loop), implement norms as sqrt(sum(x**2)), write reductions with user-defined functions and loops for lack of reduce, and know nothing of findloc, parity, iall/iany/iparity or is_contiguous. Books written for Fortran 2018 are current on this subject; the additions of Fortran 2023 concern do concurrent, rank-independent declarations and the intrinsics of other topics, not these functions.

Pitfalls

Mistake Effect Fix

Assigning maxloc(u) to a scalar integer

Compile error (Incompatible ranks 0 and 1 in assignment in GNU Fortran).

Use an integer array of size rank(u), or maxloc(u, dim=1) for a rank-one u.

Using the maxloc result as a subscript of an array with non-unit lower bounds

The wrong element, silently.

Add lbound(w) - 1.

Treating a maxloc or findloc result of 0 as a valid subscript

Out-of-bounds access or an element near the origin.

Test for 0 before using it; see "no match returns zero".

findloc on floating-point data

It finds nothing, because == is exact.

Compare with a tolerance in a loop or findloc(abs(x - v) < tol, .true.).

sum of small integer kinds

Silent overflow.

Convert to a wider kind first.

Comparing two sum results with ==

A result that changes with the optimisation level.

Compare with a tolerance.

reduce over an array that can be empty, without identity

Error termination (and, in GNU Fortran 16.2, wrong values for a zero-sized array without a mask).

Pass identity, guard on size, or use a mask.

A non-associative operation in reduce

The result depends on grouping and differs between compilers and optimisation levels.

Use an associative operation or ordered=.true..

cshift/eoshift stencils in an inner loop

An array temporary per call, with the memory traffic that goes with it.

Use array sections.

merge with an unsafe fsource

merge(a/b, 0.0, b /= 0.0) still divides by zero.

Use where or an if.

Expecting lbound of a section or expression to keep the original bounds

It is 1.

Compute with the section’s own lbound or shape.

matmul for large matrices in a performance-sensitive loop

Correct but slow without an optimised BLAS.

Link a BLAS with -fexternal-blas or call gemm.

dot_product on complex vectors

The first argument is conjugated, which sum(z * w) does not do.

Use the one that matches the mathematics.

See also

References