HPJava subarrays

<h1>Subarrays</h1>

A <em>subarray</em> is a section of a array, typically selected by some
combination of scalar and triplet subscripts.

<p>

A subarray is an object in its own right - its type is that of a suitable
multi-dimensional array.  It describes some subset of the elements of
a parent array.  The syntax for creating a subarray is modelled on the
Fortran syntax for creating regular sections, but the subscripts are
enclosed in double brackets.  This distinguishes subarray construction
clearly from the simple subscripting operation introduced earlier.  In
<font color=Blue>
<pre>
  int [[]] e = a [[2,  2 : ]] ;

  foo(b [[ : , 0, 1 : 10 : 2]]) ;
</pre>
</font>
assume <tt>a</tt> is a two-dimensional array and <tt>b</tt> is a
three-dimensional array.  <tt>e</tt> becomes an alias for the 3rd row
of elements of <tt>a</tt>.  The method <tt>foo</tt> expects a
two-dimensional array as argument.  It can read or write to the set of
elements of <tt>b</tt> selected by the section.  Note that the
subscripts of <tt>e</tt>, like any other array, start at 0, although
the first element is identified with <tt>a[2, 2]</tt>.

<p>

Subarray construction is a relatively complex operation.  Besides
scalar integers and triplets, allowed subscripts include locations and
ranges.  A location subscript must be a location in the range
associated with the relevant dimension of the array, and a range
subscript must be a subrange of the appropriate range of the array.  By
definition, integer subscript <tt>n</tt> in dimension <tt>r</tt> of
array <tt>a</tt> is equivalent to location subscript
<tt>a.rng(r)[n]</tt>, and triplet subscript <tt>l:u:s</tt> is
equivalent to range subscript <tt>a.rng(r)[l:u:s]</tt>.

<p>

The mapping (distribution group and ranges) of a subarray 
is determined as follows.  Suppose all integer and triplet subscripts
are replaced by their equivalent location or range subscripts.  If the
location subscripts are <tt>i, j, ...</tt> the distribution group of
the subarray is
<font color=Blue>
<pre>
  a.grp() / i / j / ...
</pre>
</font>
where <tt>a</tt> is the parent array.  The <tt>s</tt>'th range of the
subarray is equal to the <tt>s</tt>'th range-valued subscript.

<p>

We can use subranges and subarrays to give a parallel implementation
Cholesky decomposition.  In pseudocode the algorithm is
<font color=Blue>
<pre>
  for k=1 to n-1
    l<sub>kk</sub> = a<sub>kk</sub><sup>1/2</sup>
    for s=k+1 to n
      l<sub>sk</sub>=a<sub>sk</sub>/l<sub>kk</sub>
    for j=k+1 to n
      for i=j to n
        a<sub>ij</sub> = a<sub>ij</sub> - l<sub>ik</sub> l<sub>jk</sub>
  l<sub>nn</sub> = a<sub>nn</sub><sup>1/2</sup>
</pre>
</font>
The HPJava parallel version is
<font color=Blue>
<pre>
  Procs1 p = new Procs1(P) ;
  on(p) {
    Range x = new CyclicRange(N, p.dim(0));

    float [[*,]] a = new float [[N, x]] ;

    float [[*]]  b = new float [[N]] ;

    ... some code to initialise `a' ...

    for(int k = 0 ; k < N - 1 ; k++) {

      at(l = x[k]) {
        float d =  Math.sqrt(a[k, l]) ;

        a[k, l] = d ;
        for(int s = k + 1 ; s < N ; s++)
          a[s, l] /= d ;
      }

      Adlib.remap(b[[k + 1 : ]], a[[k + 1 :, k]]) ;

      overall(j = x[k + 1 : ])
        for(int i = x.idx(j) ; i < N ; i++)
          a[i, j] -= b[i] * b[x.idx(j)] ;
    }

    at(l = x [N - 1])
      a[N - 1, l] = Math.sqrt(a[N - 1, l]) ;
  }
</pre>
</font>
A cyclic range is used to allocate the original array on
multiprocessors.  The result lower triangular matrix overwrites part of
its storage.  During the computation, the collective communication
function <tt>remap</tt> is used for broadcasting updated columns.

<p>

The example also illustrates the use of the member
<font color=Blue>
<pre>
  int idx(Location)
</pre>
</font>
of <tt>Range</tt>, which recovers a global subscript from a location. 

<!--
<b>Next:</b>
<a href="subarrays.html">Subarrays</a>
-->


