Much can be achieved using the array functions of Fortran alone. Imperative (do/for) loops, for example, can often be avoided. The following code performs a one-dimensional convolution using array *sections*:

a(2:n-1) = a(1:n-2) + a(3:n)

There are also a few array “reduction” functions; such as *all*, *any*, *count*, *maxval*, *minval*, *product*, and *sum*. A higher-order reduction function could implement all of these functions; alas (in Fortran) it doesn’t *a priori *exist. For example, the functionality of *product* could be achieved using syntax such as:

result = reduce( (*), 1, [1,2,3] )

where *result* would be assigned a value of 6.

A common and useful partner to the *reduce* function is the *scan* function. The *scan* function returns not *one* final value, but all intermediate values too. If it existed in Fortran, a use case might look like this:

result = scan( (*), 1, [1,2,3] )

Assuming the multiplicative *unit* (1) is not included in the output, this would produce a result of [ 1, 2, 6 ]. I’d now like to implement this particular (multiplicative) scan using only existing Fortran array intrinsic functions; i.e. no loops. The earlier input of [ 1, 2, 3 ] will be used as a working example. (We will return to this scan function in the following post on diagonal selection.)

## Creating the Scan Expression

The intrinsic function *spread* will, given an array of rank n, produce a new array with rank n+1. The new array is constructed by replicating the input a number of times specified by *spread*‘s third (and final) argument. The second argument to *spread* specifies which dimension (>=1) to replicate over. We can visualise this choice with a pencil representing our 1D input. If we wish to replicate the pencil to make a square, we can begin *either* with the pencil positioned horizontally; *or,* vertically. Similarly, if we instead had a book and wished to create a cuboid, we begin with three choices.

A simple way to visualise the output of *spread*, is simply to consider the resulting arrangement of array values in memory. Often, Fortran array data is arranged in a contiguous linear arrangement, and so the output of:

spread( [1,2,3], 1, 3 )

is stored in memory as [ 1, 1, 1, 2, 2, 2, 3, 3, 3 ]. The output of:

spread( [1,2,3], 2, 3 )

however, is stored as [ 1, 2, 3, 1, 2, 3, 1, 2, 3 ]. Both are 2D arrays with a shape of [ 3, 3]; and the difference is clear. Here we take the second choice, and replicate the input vector (that is, a 1D array) across the second dimension.

We now use the *eoshift* intrinsic function to shift the individual rows of this 3×3 array to the right, filling the holes left behind with the value provided by *eoshift*‘s third (and final) argument, 1; aka. the multiplicative unit. The second argument is a vector indicating the degree of shifting required for each row. As the default shift direction is to the *left*, we are looking for a vector, [ -2, -1, 0 ]. It may be constructed using the following *implied-do* expression:

[(i,i=-(size([1,2,3])-1),0)]

So, obtaining the length of the input using the *size* intrinsic, we have built a new 1D array containing the arithmetic progression from *minus* one less than the length, up to zero.

Returning to *eoshift*, the following expression:

eoshift( spread( [1,2,3] ,2, size([1,2,3]) ), & [(i,i=-(size([1,2,3])-1),0)], & 1)

evaluates to a 2D array with shape [ 3, 3 ]; its contents in memory are: [ 1, 1, 1, 1, 1, 2, 1, 2, 3 ].

Fittingly, the *product* reduction intrinsic can at last be employed. As well as operating on a vector, *product* can be applied to all vectors aligned with the axis specified by its second (and optional) argument. In this case that is axis #1, and so final result can be obtained by applying the *product* intrinsic to the result of the *eoshift* expression from above:

product( eoshift( spread( [1,2,3] ,2, size([1,2,3]) ), & [(i,i=-(size([1,2,3])-1),0)], & 1), & 1)

## Wrapping up as a Function

Rather than crafting such an expression inline around each 1D array expression we wish the scanned product of, we can at least create a function; *scanl_product*. The *execution-part* of this function is basically the *product* expression above, with a 1D array variable abstraction *v* replacing the literal [ 1, 2, 3 ]:

function scanl_product(v) result(vout) integer, intent(in) :: v(:) integer :: i, vout(size(v)) vout = product( eoshift( spread( v, 2, size(v) ), & [(i,i=-(size(v)-1),0)], & 1), & 1) end function