Array Handling

by Michael Metcalf / CERN CN-AS

Array handling is included in Fortran 90 for two main reasons:

the notational convenience it provides, bringing he code closer to the underlying mathematical form;
for the additional optimization opportunities it gives compilers (although there are plenty of opportunities for degrading optimization too!).

At the same time, major extensions of the functionality in this area have been added. We have already met whole arrays in Parts 2 and 3 of this series - here we develop the theme.

Zero-sized arrays

A zero-sized array is handled by Fortran 90 as a legitimate object, without special coding by the programmer. Thus, in

     DO i = 1,n
        x(i) = b(i) / a(i, i)
        b(i+1:n) = b(i+1:n) - a(i+1:n, i) * x(i)
     END DO

no special code is required for the final iteration where i = n. We note that a zero-sized array is regarded as being defined; however, an array of shape (0,2) is not conformable with one of shape (0,3), whereas

    x(1:0) = 3

is a valid 'do nothing' statement.

Assumed-shape arrays

These are an extension and replacement for assumed-size arrays. Given an actual argument like:

     REAL, DIMENSION(0:10, 0:20) :: a
        :
     CALL sub(a)

the corresponding dummy argument specification defines only the type and rank of the array, not its size. This information has to be made available by an explicit interface, often using an interface block (see Part 5). Thus we write just

  SUBROUTINE sub(da)
     REAL, DIMENSION(:, :) :: da

and this is as if da were dimensioned (11,21). However, we can specify any lower bound and the array maps accordingly. The shape, not bounds, is passed, where the default lower bound is 1 and the default upper bound is the corresponding extent.

Automatic arrays

A partial replacement for the uses to which EQUIVALENCE is put is provided by this facility, useful for local, temporary arrays, as in

     SUBROUTINE swap(a, b)
        REAL, DIMENSION(:)       :: a, b
        REAL, DIMENSION(SIZE(a)) :: work
        work = a
        a = b
        b = work
     END SUBROUTINE swap

The actual storage is maintained on a stack.

ALLOCATABLE and ALLOCATE

Fortran 90 provides dynamic allocation of storage; it relies on a heap storage mechanism (and replaces another use of EQUIVALENCE). An example, for establishing a work array for a whole program, is

     MODULE work_array
        INTEGER n
        REAL, DIMENSION(:,:,:), ALLOCATABLE :: work
     END MODULE
     PROGRAM main
        USE work_array
        READ (input, *) n
        ALLOCATE(work(n, 2*n, 3*n), STAT=status)
        :
        DEALLOCATE (work)

The work array can be propagated through the whole program via a USE statement in each program unit. We may specify an explicit lower bound and allocate several entities in one statement. To free dead storage we write, for instance,

     DEALLOCATE(a, b)

We will meet this later, in the context of pointers.

Elemental operations and assignments

We have already met whole array assignments and operations:

   REAL, DIMENSION(10) :: a, b
   a = 0.          ! scalar broadcast; elemental assignment
   b = sqrt(a)     ! intrinsic function result as array object

In the second assignment, an intrinsic function returns an array-valued result for an array-valued argument. We can write array-valued functions ourselves (they require an explicit interface):

   PROGRAM test
      REAL, DIMENSION(3) :: a = (/ 1., 2., 3./),       &
                            b = (/ 2., 2., 2. /),  r
      r = f(a, b)
      PRINT *, r
   CONTAINS
      FUNCTION f(c, d)
      REAL, DIMENSION(:) :: c, d
      REAL, DIMENSION(SIZE(c)) :: f
      f = c*d        ! (or some more useful functio of c and d)
      END FUNCTION f
   END PROGRAM test

WHERE

Often, we need to mask an assignment. This we can do using the WHERE, either as a statement:

     WHERE (a /= 0.0) a = 1.0/a  ! avoid division by 0

(note: test is element-by-element, not on whole array), or as a construct:

     WHERE (a /= 0.0)
        a = 1.0/a
        b = a             ! all arrays same shape
     END WHERE

     WHERE (a /= 0.0)
        a = 1.0/a
     ELSEWHERE
        a = HUGE(a)
     END WHERE

Array elements

Simple case: given

         REAL, DIMENSION(100, 100) :: a

we can reference a single element as, for instance, a(1, 1). For a derived-data type like

     TYPE triplet
        REAL                  u
        REAL, DIMENSION(3) :: du
     END TYPE triplet

we can declare an array of that type:

     TYPE(triplet), DIMENSION(10, 20) :: tar

and a reference like

                     tar(n, 2)

is an element (a scalar!) of type triplet, but

                     tar(n, 2)%du

is an array of type real, and

                     tar(n, 2)%du(2)

is an element of it. The basic rule to remember is that an array element always has a subscript or subscripts qualifying at least the last name.

Array subobjects (sections)

The general form of subscript for an array section is

 
       [lower] : [upper] [:stride]

as in

       REAL a(10, 10)
       a(i, 1:n)                ! part of one row
       a(1:m, j)                ! part of one column
       a(i, : )                 ! whole row
       a(i, 1:n:3)              ! every third element of row
       a(i, 10:1:-1)            ! row in reverse order
       a( (/ 1, 7, 3, 2 /), 1)  ! vector subscript
       a(1, 2:11:2)             ! 11 is legal as not referenced
       a(:, 1:7)                ! rank two section

Note that a vector subscript with duplicate values cannot appear on the left-hand side of an assignment as it would be ambiguous. Thus,

       b( (/ 1, 7, 3, 7 /) ) = (/ 1, 2, 3, 4 /)

is illegal. Also, a section with a vector subscript must not be supplied as an actual argument to an OUT or INOUT dummy argument. Arrays of arrays are not allowed:

       tar%du             ! illegal

We note that a given value in an array can be referenced both as an element and as a section:

       a(1, 1)            !  scalar (rank zero)
       a(1:1, 1)          !  array section (rank one)

depending on the circumstances or requirements. By qualifying objects of derived type, we obtain elements or sections depending on the rule stated earlier:

       tar%u              !  array section (structure component)
       tar(1, 1)%u        !  component of an array element

Arrays intrinsic functions

Vector and matrix multiply

      DOT_PRODUCT        Dot product of 2 rank-one arrays
      MATMUL             Matrix multiplication

Array reduction

      ALL                True if all values are true
      ANY                True if any value is true. Example:
                             IF (ANY( a > b)) THEN
      COUNT              Number of true elements in array
      MAXVAL             Maximum value in an array
      MINVAL             Minimum value in an array
      PRODUCT            Product of array elements
      SUM                Sum of array elements

Array inquiry

      ALLOCATED          Array allocation status
      LBOUND             Lower dimension bounds of an array
      SHAPE              Shape of an array (or scalar)
      SIZE               Total number of elements in an array
      UBOUND             Upper dimension bounds of an array

Array construction

      MERGE              Merge under mask
      PACK               Pack an array into an array of rank
      SPREAD             Replicate array by adding a dimension
      UNPACK             Unpack an array of rank one into an array under mask

Array reshape

      RESHAPE            Reshape an array

Array manipulation

      CSHIFT             Circular shift
      EOSHIFT            End-off shift
      TRANSPOSE          Transpose of an array of rank two

Array location

      MAXLOC             Location of first maximum value in an array
      MINLOC             Location of first minimum value in an array