Scanners!
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 3x3 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