Fortran, Uncategorized

Diagonal Selection

A good example of the expressivity of the forall statement from Fortran 95 and HPF, is the selection of the elements on the diagonal of an array. We begin such a thing with the declarations of a rank 2 array a, and a vector da to hold the diagonal:

integer :: i, a(4,2) = reshape([(i,i=1,size(a))],shape(a))
integer :: da(minval(shape(a)))

Note that regardless of the rank of the original array, the length of the sought diagonal will always be equal to the mininum element of the array’s shape vector. In our case the shape vector of a is [ 4, 2 ], and so the diagonal has a length of 2. Of course this is due to our definition of the diagonal of an array. Informally, our notion of diagonal corresponds to the set of array elements, indexed by a set of equal values; and ordered by those values. So, for a, this will be the elements indexed by a(1,1) and a(2,2): [ 1, 6 ]. So then, let’s look at the forall statement to achieve this:

forall (i=1:size(da)) da(i) = a(i,i)

How would this look as an array expression? As discussed in my first blog post, such array expressions hold the potential to handle input of arbitrary dimensionality; unlike the forall statement; though certainly more verbose.

Diagonal Selection

If an array of rank n is flattened, say by the Fortran pack intrinsic, it becomes clear that the number of elements separating each adjacent pair of diagonal members, is constant. Conveniently, the set of recursive equations which can express this constant may also be encoded using the scan_product function from the previous blog post.

The array a has shape [ 4, 2 ]. The product scan of [ 4, 2 ] is [ 4, 8 ]. We’re looking for a logical mask vector to select the diagonal elements. ( True and False will also be denoted as 1 or  0 as required. ) Some examples: the earlier [ 4, 2 ] shape array will require a mask of [ 1, 0, 0, 0, 0, 1, 0, 0 ]; while an array of shape [ 3, 4 ] requires a mask of [ 1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0 ]. The pattern to notice is a repeating string of False values, delimited by True values, then succeeded by False values. The repeating string of False values has a length of four in the first example, and a length of three in the second. How many False values in general? With an input array a, the following code will assign the answer to the zps (zeros_per_step) integer variable:

zps = sum(eoshift(scan_product(shape(a)),-1))

This value be used to size a logical vector (zs) of False values as shown. A True value should be cons’d to the start of this vector, and the result repeated mi times, where mi is one less than the length of the final diagonal. The implicit do loop can handle such concatenative actions concisely, as shown in the code below:

zs  = spread(.false.,1,zps)
mi  = minval(shape(a))-1
tzs = [(.true.,zs,j=1,mi)]

The final, padding, string of False values, the vector zpd below, is longest in arrays having a relatively large final dimension; for example, shapes such as [ 2, 10 ], or [ 4, 5, 100 ].

nz  = product(shape(a)) - (mi * (zps+1)) - 1
zpd = spread(.false.,1,nz)

To construct the final 1D boolean selection mask, msk, the tzs vector has both a single True value appended; and also the zpd vector:

msk  = [(tzs,.true.,zpd,i=1,1)]

To apply this mask vector, we need to flatten the input array a using the pack intrinsic. The same pack intrinsic can then be used to apply the mask we have constructed above:

pack(pack(a,.true.),msk)

If the operation is to be encapsulated in a function, we are once again forced to choose a single, static rank for the result. The code below is then arbitrarily for an array of rank 3; again, an interface statement, could help construct a verbose generic representation. I’ve also used the associate construct from Fortran 2003, which operates like a non-recursive let expression. Hopefully it can aid readability; at the least one can easily print out the contents of each binding.

function diag3(a) result(r)
  integer, intent(in) :: a(:,:,:)
  integer :: i,j,r(minval(shape(a)))

  associate (sa  => shape(a)         )
  associate (sc  => scanl_product(sa))
  associate (eo  => eoshift(sc,-1)   )
  associate (zps => sum(eo)          )
  associate (ps  => product(sa)      )
  associate (mi  => minval(sa)-1     )
  associate (nz  => ps - (mi * (zps+1))-1)
  associate (zpd => spread(.false.,1,nz))
  associate (zs  => spread(.false.,1,zps))
  associate (tzs =>    [(.true.,zs,j=1,mi)])
  associate (msk => [(tzs,.true.,zpd,i=1,1)])
    r = pack(pack(a,.true.),msk)
  end associate
  end associate
  end associate
  end associate
  end associate
  end associate
  end associate
  end associate
  end associate
  end associate
  end associate
end function
Fortran, Uncategorized

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 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
coz = seventeen1(z)
Fortran, Uncategorized

Fragrant Fortran

Fortran is a programming language like C, with more capital letters, and less pointers. Since around 1990, the Fortran standard has also supported arrays as first class language elements. Such array features, provided by Fortran 90, Fortran 95, and Fortran 2003, are reminiscent of Kenneth Iverson’s seminal APL (A Programming Language), and J; though less developed.

Fortran 9x includes many useful and essential procedures for manipulating arrays. For example, the function shape takes one array argument (of any dimensionality) and return a one-dimensional array (often referred to as a vector) which lists the extent/size of each dimension. The length of this vector is therefore equal to the dimensionality of the argument to the shape function.

Tens of such array functions are specified by the Fortran standard, and included in essentially every Fortran compiler. Along with the more numerous scalar functions, these are collectively referred to as the Fortran intrinsics.

Fortran array intrinsics are often defined for arrays of arbitrary dimensionality (rank). However, the end user is made to work hard to replicate such genericity. A user-defined function accepting or returning arrays must specify a fixed dimensionality. Consequently, Fortran’s expression language is massively more concise than its procedure language. Hence,

m = sum(a) / size(a)

will assign the mean (excepting zero-sized arrays, and other desiderata) of input array a, to m, regardless of the rank of the a. A function to perform a similar action on a one-dimensional array could look like this:

function mean1(a) result(m)
  real, intent(in) :: a(:)
  real :: m

  m = sum(a) / size(a)
end function

which could be called with:

m = mean1(a)

We would need another for an array of rank 2:

function mean2(a) result(m)
  real, intent(in) :: a(:,:)     ! Extra colon !
  real :: m

  m = sum(a) / size(a)
end function

The two versions of mean differ only by the rank in array a‘s declaration. Also, apart from such “boilerplate” declarations, the core of our function is lexically the same as the earlier expression; “m = sum(a) / size(a)”. At least Fortran arrays have a maximum rank of seven, so we would “only” need seven versions of this mean function. Tragically however, if a function takes n array arguments, we will need to write 7^n functions.

One ray of light comes in the form of Fortran’s interface statement. The interface statement allows a user to access multiple procedures via a single name; similar to C++ overloading. For example, the following interface statement allows a call to mean to be instantiated appropriately, determined by the types of its arguments.

interface mean
  module procedure mean1, mean2
end interface mean

It’s worth mentioning that Professor Martin Erwig’s Parametric Fortran project provides one solution to this convoluted aspect of Fortran.

Finally, note also that alongside array expressions, 1990 also eschewed ye olde brutal and capitalist “FORTRAN” in favour of the effervescent stateswoman we know and love as “Fortran” today.