To: J3 J3/18-248
From: Van Snyder
Subject: Supporting rank variability
Date: 2018-September-13
References: 18-245, 18-247
Introduction
============
It is occasionally necessary to compute element position from subscripts
and bounds, or to compute subscripts from element position and bounds.
Here's an example.
A geophysical quantity might be represented in one, two, or three
physical dimensions. It might or might not depend upon frequency. It
might or might not have more than one component (e.g. magnetic field).
If it is necessary to interpolate from that representation to another,
for example a line of sight, interpolation coefficients are required.
Interpolation is a matrix-vector product. The interpolation
coefficients are best represented by a sparse matrix. Off-the-shelf
sparse-matrix software has only a single column index for each element,
not a tuple of up to five subscripts. In order to compute one element
of the matrix-vector product, an iterator (see 18-245) is used to
traverse a row of the interpolation matrix. Each element it produces
consists of the value, the row index, and the column index. To compute
the product, the column index of the matrix element needs to be
converted to subscripts of the geophysical state.
Proposal
========
Provide two intrinsic functions.
One, tentatively named ELEMENT_POSITION here, has SUBSCRIPTS and UBOUNDS
arguments, and an optional LBOUND argument (values assumed all 1 if it
is absent). The result value is the scalar element position of an array
element at position SUBSCRIPT in an array having shape
( [ LBOUNDS : ] UBOUNDS ) (see 18-247).
The other, tentatively named SUBSCRIPTS here, has ELEMENT_POSITION and
UBOUNDS arguments, and an optional LBOUND argument (values assumed all 1
if it is absent). The result value is a vector of subscripts that can
be used to access the element of an array having shape
( [ LBOUNDS : ] UBOUNDS ) (see 18-247) at the specified element
position.
These are relatively trivial calculations but they are more expensive
than the element multiplication. If implemented by users as module
procedures, it is unlikely processors will inline them (and users might
need several attempts to get the calculations done correctly). If they
are provided by the processor, which already knows how to do at least
the ELEMENT_POSITION calculation, as intrinsic functions, one might hope
that at least some processors will inline them.
Caveat
======
One might believe that the problem can be solved by representing the
geophysical quantity as a one-dimensional (or at least contiguous)
object, use pointer rank remapping to access it as a geophysical
quantity, and use the one-dimensional representation for interpolation.
But it is sometimes necessary to interpolate only in a subset of
discontiguous dimensions. For example, a quantity that depends upon
frequency and two geophysical dimensions might be represented as a
three-dimensional array ( frequency, vertical position, horizontal
position ). In some cases, it is necessary to interpolate in
all three dimensions, in which case rank remapping might be usable. In
other cases, it is only necessary to interpolate in geophysical
position, the same way for every frequency. The column indices of the
elements of the interpolation coefficient matrix in the second case
would depend only upon geophysical coordinates.
One might be tempted to suggest that the column indices in the elements
of the interpolation coefficient matrix ought to depend upon all three
dimensions of the geophysical quantity. In many cases, however, there
are several geophysical quantities, all having the same geophysical
locations, some depending upon frequency and others not depending upon
frequency. For geophysical interpolation alone, the same interpolation
coefficients can be used for all quantities. A "solution" that requires
another matrix of elements all having the same values and row indices,
but different column indices, is not helpful.
One might also be tempted to recommend that "all ya gotta do is" permute
the order of dimensions. In a legacy code of hundreds of thousands or
millions of lines, this is not a helpful suggestion. And it doesn't
always solve the problem, as for some calculations one order allows to
consider a contiguous (i.e., leading) subset of dimensions, while for
other calculations using the same object, a different subset of
dimensions -- discontiguous ones -- are used.
And, of course, there is the possibility that the geophysical object is
represented by an assumed-shape dummy argument that is not contiguous,
even in its leading dimensions. Working with it as if it were
contiguous would require copy-in or VALUE argument association.