3.3. Element-wise Array Implementation

The simplest way to use the POOMA Toolkit is to use the POOMA Array class instead of C arrays. Arrays automatically handle memory allocation and deallocation, support a wider variety of assignments, and can be used in expressions. Example 3-2 implements Doof2d using Arrays and element-wise accesses. Since the same algorithm is used as Example 3-1, we will concentrate on the differences.

Example 3-2. Element-wise Array Implementation of Doof2d


#include <iostream>		// has std::cout, ...
#include <stdlib.h>		// has EXIT_SUCCESS
#include "Pooma/Arrays.h"
    // has POOMA's Array declarations  (1)

// Doof2d: POOMA Arrays, element-wise implementation

int main(int argc, char *argv[])
{
  // Prepare the POOMA library for execution.  (2)
  Pooma::initialize(argc,argv);
  
  // Ask the user for the number of averagings.
  long nuAveragings, nuIterations;
  std::cout < < "Please enter the number of averagings: ";
  std::cin >  > nuAveragings;
  nuIterations = (nuAveragings+1)/2;
    // Each iteration performs two averagings.

  // Ask the user for the number n of values along
  // one dimension of the grid.
  long n;
  std::cout < < "Please enter the array size: ";
  std::cin >  > n;

  // Specify the arrays' domains [0,n) x [0,n).  (3)
  Interval<1> N(0, n-1);
  Interval<2> vertDomain(N, N);

  // Create the arrays.  (4)
  // The Array template parameters indicate
  // 2 dimensions, a 'double' value
  // type, and ordinary 'Brick' storage.
  Array<2, double, Brick> a(vertDomain);
  Array<2, double, Brick> b(vertDomain);

  // Set up the initial conditions.
  // All grid values should be zero except for the
  // central value.  (5)
  for (int j = 1; j < n-1; j++)
    for (int i = 1; i < n-1; i++)
      a(i,j) = b(i,j) = 0.0;
  b(n/2,n/2) = 1000.0;

  // In the average, weight elements with this value.
  const double weight = 1.0/9.0;

  // Perform the simulation.
  for (int k = 0; k < nuIterations; ++k) {
    // Read from b.  Write to a.
    for (int j = 1; j < n-1; j++)
      for (int i = 1; i < n-1; i++)
        a(i,j) = weight *  (6)
          (b(i+1,j+1) + b(i+1,j  ) + b(i+1,j-1) +
           b(i  ,j+1) + b(i  ,j  ) + b(i  ,j-1) +
           b(i-1,j+1) + b(i-1,j  ) + b(i-1,j-1));

    // Read from a.  Write to b.
    for (int j = 1; j < n-1; j++)
      for (int i = 1; i < n-1; i++)
        b(i,j) = weight *
          (a(i+1,j+1) + a(i+1,j  ) + a(i+1,j-1) +
           a(i  ,j+1) + a(i  ,j  ) + a(i  ,j-1) +
           a(i-1,j+1) + a(i-1,j  ) + a(i-1,j-1));
  }

  // Print out the final central value.
  Pooma::blockAndEvaluate();
    // Ensure all computation has finished.
  std::cout < <
     (nuAveragings % 2 ? a(n/2,n/2) : b(n/2,n/2))
     < < std::endl;

  // The arrays are automatically deallocated.  (7)

  // Tell the POOMA library execution finished.  (8)
  Pooma::finalize();
  return EXIT_SUCCESS;
}
(1)
To use POOMA Arrays, the Pooma/Arrays.h must be included.
(2)
The POOMA Toolkit structures must be constructed before their use.
(3)
Before creating an Array, its domain must be specified. The N Interval represents the one-dimensional integral set {0, 1, 2, …, n-1}. The Interval<2> vertDomain object represents the entire two-dimensional index domain.
(4)
An Array's template parameters indicate its dimension, its value type, and how the values will be stored or computed. The Brick Engine type indicates values will be directly stored. It is responsible for allocating and deallocating storage so new and delete statements are not necessary. The vertDomain specifies the array index domain.
(5)
The first loop initializes all Array values to the same scalar value. The second statement illustrates assigning one Array value. Indices, separated by commas, are surrounded by parentheses rather than surrounded by square brackets ([]).
(6)
Array element access uses parentheses, rather than square brackets.
(7)
The Arrays deallocate any memory they require, eliminating memory leaks.
(8)
The POOMA Toolkit structures must be destructed after their use.

We describe the use of Array and the POOMA Toolkit in Example 3-2. Arrays, declared in the Pooma/Arrays.h, are first-class objects. They "know" their index domain, can be used in expressions, can be assigned scalar and array values, and handle their own memory allocation and deallocation.

The creation of the a and b Arrays requires an object specifying their index domains. Since these are two-dimensional arrays, their index domains are also two-dimensional. The two-dimensional Interval<2> object is the Cartesian product of two one-dimensional Interval<1> objects, each specifying the integral set {0, 1, 2, …, n-1}.

An Array's template parameters indicate its dimension, the type of its values, and how the values are stored. Both a and b are two-dimension arrays storing doubles so their dimension is 2 and their value type is double. An Engine stores an Array's values. For example, a Brick Engine explicitly stores all values. A CompressibleBrick Engine also explicitly stores values if more than one value is present, but, if all values are the same, storage for just that value is required. Since an engine can store its values any way it desires, it might instead compute its values using a function or compute using values stored in separate engines. In practice, most explicitly specified Engines are either Brick or CompressibleBrick.

Arrays support both element-wise access and scalar assignment. Element-wise access uses parentheses, not square brackets. For example, b(n/2,n/2) specifies the central element. The scalar assignment b = 0.0 assigns the same 0.0 value to all array elements. This is possible because the array knows the extent of its domain. We illustrate these data-parallel statements in the next section.

Any program using the POOMA Toolkit must initialize the toolkit's data structures using Pooma::initialize(argc,argv). This extracts POOMA-specific command-line options from the program's command-line arguments and initializes the interprocessor communication and other data structures. When finished, Pooma::finalize() ensures all computation and communication has finished and the data structures are destructed.