To: J3                                                     J3/18-247
From:    Van Snyder
Subject: Supporting rank genericity
Date: 2018-September-13
References: 04-195

Introduction
============

MAXLOC(A) returns a rank-1 array whose extent is equal to the rank of A.
It cannot be used directly as a subscript for A.  One must assign it to
a variable of the appropriate rank, and use elements of that array as
subscripts.  For example

  integer :: S(3)
  real :: A(10,10,10)
  ...
  s = maxloc(a)
  print *, 'Maximum value in A is ', a(s(1),s(2),s(3))

UBOUND(A) returns a rank-1 array whose extent is equal to the rank of A.
It cannot be used directly as upper bounds for an automatic variable.
For example

  real, intent(in) :: A(:,:,:)
  real :: X(ubound(a,1),ubound(a,2),ubound(a,3))

or

  real, allocatable :: X(:,:,:)
  ...
  allocate ( x(ubound(a,1),ubound(a,2),ubound(a,3)) )

or, if you need a halo around all dimensions of X

  allocate ( x(lbound(a,1)-1:ubound(a,1)+1, &
           &   lbound(a,2)-1:ubound(a,2)+1, &
           &   lbound(a,3)-1:ubound(a,3)+1) )

These are tedious if you know the rank, and impossible if we develop
facilities for rank-generic programming.

Proposal
========

1. Array bounds
---------------

Allow a single array B of rank one and extent R to specify the upper or
bounds of a rank R array A.  The extent of B shall be a constant, except
maybe for the A(..) case.  For example

  real, intent(in) :: A(:,:,:)
  real :: X(ubound(a))

or

  real, allocatable :: X(:,:,:)
  ...
  allocate ( x(ubound(a)) )

or (for the halo case)

  allocate ( x(lbound(a)-1:ubound(a)+1) )

The case of one bound being a scalar and the other not is obvious:

  allocate ( x(0:ubound(a)) )

is the same as

  allocate ( x([0,0,0]:ubound(a)) )

Indeed, the default case (above) of lower bounds being all 1

  allocate ( x(ubound(a)) )

is the same as

  allocate ( x(1:ubound(a)) )

2. Subscripts
-------------

Allow a single array of rank one and extent R as the subscript of a
rank R array.  For example

  real :: A(10,10,10)
  ...
  print *, 'Maximum value in A is ', a(maxloc(a))

Then A(S) is the same as A(3,4,5), not A(3:3,4:4,5:5).

Let's not do kludges and half measures, which is how we got nonpointer
optional arguments corresponding to disassociated pointers being
absent.  The original solution proposed for the combinatorial-explosion
problem was a special case of the conditional expression proposal, with
an empty ELSE part meaning "if the logical expression is true, the
argument is the consequent (not its value), else the argument is
absent."

The obvious generalization here is to allow a single subscript S of rank
K+1 and extents (R,N1, ..., NK) as a subscript for a rank-R array A.
The elements of the rank-one sections in the first dimension of S are
used consecutively as subscripts for A, as in the rank one case,
resulting in a rank K array of extents N1, ..., NK.

This provides a more general scatter/gather facility than the present
vector subscript facility.  This is not the same as using the elements
of the rank-one sections in S as vector subscripts for A, which would
result in a rectangular section of shape [n,n] (in the case A is
of rank 2 and S is of shape [2,n]).

Example:

Suppose we have arrays A3 with dimensions (10,10,10) and S3 with
dimensions (3,2).  If we assume

                                           [ 3 6 ]
S3 = reshape( [3, 4, 5, 6, 7, ], [3,2] ) = [ 4 7 ]
                                           [ 5 8 ]

then A3(S3) is a rank-1 extent-2 array that can appear in a variable-
definition context (except as an actual argument associated with a dummy
argument having INTENT(OUT) or INTENT(INOUT)); it specifies the same
array as [ A3(3,4,5), A3(6,7,8) ], which cannot appear in a
variable-definition context.  This is different from
A3(S3(1,:),S3(2,:),S3(3,:)), which can appear in a variable-definition
context, but is an object with shape [2,2,2], not [2].  The former is an
arbitrary collection of elements of A3, while the latter is a
rectangular section of A3.

The degenerate case (same as above) is

  real :: A(10,10,10)
  integer :: S(3) = [ 3, 4, 5 ]

Then A(S) is the same as A(3,4,5), not A(3:3,4:4,5:5).  This illustrates
that a rank K+1 subscript with first extent R for a rank R array
produces a rank K result.  In this example, K = 0, so A(S) is a scalar,
not a rank-3 array with shape [1,1,1].  I.e., S is NOT a vector
subscript.  But many of the same restrictions for vector subscripts
remain, e.g. not being an actual argument corresponding to an
INTENT(OUT) dummy argument.

This can probably be made to work with DIMENSION(..) arrays.  A run-time
check is probably required to make sure that in A(S), the first extent
of S is the same as the (dynamic) rank of A.  This might be especially
useful to create automatic variables with the same shape as
DIMENSION(..) arrays.