Final Exam¶

Writing Jacobi as an Iteration with Matrix-Vector Product¶

In lecture we noted that we could express the “in-place” iterative update of the grid variables this way:

int jacobi(Grid& Phi0, size_t max_iters) {
  Grid Phi1(Phi0);
  for (size_t iter = 0; iter < max_iters; ++iter) {
    for (size_t i = 1; i < Phi0.numX()-1; ++i) {
      for (size_t j = 1; j < Phi.numY()-1; ++j) {
        Phi1(i, j) = (Phi0(i-1, j) + Phi0(i+1, j) + Phi0(i, j-1) + Phi0(i, j+1))/4.0;
      }
    }
    swap(Phi0, Phi1);
  }
  return max_iters;
}

An equivalent expression would be the following:

int jacobi(Grid& Phi, size_t max_iters) {
  Grid R(Phi);
  for (size_t iter = 0; iter < max_iters; ++iter) {
    for (size_t i = 1; i < Phi.numX()-1; ++i) {
      for (size_t j = 1; j < Phi.numY()-1; ++j) {
        R(i, j) = Phi(i,j) - (Phi(i-1, j) + Phi(i+1, j) + Phi(i, j-1) + Phi(i, j+1))/4.0;
      }
    }
    Phi = Phi - R;
  }
  return max_iters;
}

But notice what the two inner loops are: we are computing R Phi as an explicit stencil operation to effect \(R = A \times \phi\)

Let’s refactor and rewrite jacobi in a slightly more reusable and modular fashion. First, let’s define an operator between a stencil and a grid that applies the Laplacian.

Grid operator*(const Stencil& A, const Grid& x) {
  Grid y(x);

  for (size_t i = 1; i < x.numX()-1; ++i) {
    for (size_t j = 1; j < x.numY()-1; ++j) {
      y(i, j) = x(i,j) - (x(i-1, j) + x(i+1, j) + x(i, j-1) + x(i, j+1))/4.0;
    }
  }
  return y;
}

In general, we might want to carry this iteration with the stencil, and just dispatch to it from here (and use subtype or parametric polymorphism), but we are going to just use the Laplacian stencil here.

Now we can define a simple iterative solver (“ir” for iterative refinement):

size_t ir(const Stencil& A, Grid& x, const Grid& b, size_t max_iter) {
  for (size_t iter = 0; iter < max_iter; ++iter) {
    Grid r = b - A*x;
    x += r;
  }
  return max_iter;
}

Question for the reader. What do we need stencil to provide in this case? How would we define it? (Answer this question for yourself before you look into Stencil.hpp.

Second question for the reader. What functions do you need to define so that you can compile and run the code just as it is above – meaning, using matrix and vector notation? Again, answer this to yourself before looking in Grid.cpp.

(NB: Using the ir routine above is only equivalent to jacobi when we have unit diagonals in the matrix, otherwise we would need the additional step of doing that division (the general practice of which is known as preconditioning). But we’ve defined things here so that that they will work without needing to worry about preconditioning for now – even though preconditioning turns out to be a very important issue in real solvers.)

Compute cluster¶

The compute cluster for this assignment is set up in the Amazon Web Services (AWS) cloud just for this course and will be taken down once the course is completed. (NB: If you are interested in experimenting with parallel computing on your own or for a research project after this course, please let me know directly and I will try to see if we can make some arrangements to keep some of the clusters running a bit longer.)

The current IP address is: 34.209.177.35

The compute cluster is a shared compute resource, which sharing includes the compute resources themselves, as well as the file system. When you are logged in to the head node, other students will also be logged in at the same time, so be considerate of your usage. In particular,

• Use the cluster only for assignments in this course (no mining bitcoins or other activities for which a cluster might seem useful),

• Do not run anything compute-intensive on the head node,

• Do not run anything using mpirun on the head node,

• Do not put any really large files on the cluster, and

• Do not attempt to look at other students’ work.

Sanity Check(s)¶

Once you are connected to the head node, there are a few quick sanity checks you can do to check out (and learn more about) the cluster environment.

To find out more about the head node itself, try

% more /proc/cpuinfo

This will return copious information about the cores on the head node. Of interest are how many cores there are.

You should also double check the C++ compiler, the mpic++ compiler, and mpirun.

% c++ --version
% mpic++ --version

These should both print information about the installed version of clang. And, finally, issuing

% mpirun --version

should print out a message indicating that it is Open MPI.

Batch queue¶

We will be using the Slurm batch queuing system for this assignment, which we will discuss in more detail below. To get some information about it, you can issue the command

% sinfo

This will print a summary of the compute nodes along with some basic details about their capabilities.

If you pass the -N argument to sinfo you will get the info printed on a per-node basis. Passing the -l (ell) flag will provide additional information.

You may see a list that looks something like this

% sinfo -N -l
Sat Jun  8 05:04:34 2019
NODELIST         NODES PARTITION       STATE CPUS    S:C:T MEMORY TMP_DISK WEIGHT AVAIL_FE REASON
ip-10-11-10-44       1  compute*        idle    2    2:1:1   3711    16810      1   (null) none
ip-10-11-10-197      1  compute*        idle    2    2:1:1   3711    16810      1   (null) none
ip-10-11-10-201      1  compute*        idle    2    2:1:1   3711    16810      1   (null) none
ip-10-11-10-219      1  compute*        idle    2    2:1:1   3711    16810      1   (null) none

(“man sinfo” to access documentation to explain the different columns.)

The command “squeue” on the other hand will provide information about the jobs that are currently running (or waiting to be run) on the nodes.

% squeue
% squeue -arl

(“man squeue” to access information about the squeue command.)

Investigating the Cluster Queues¶

One program that is interesting to run across the cluster is the hostname command because it will do something different on different nodes (namely, print out the unique hostname of that machine). This can be very useful for quickly diagnosing issues with a cluster – and for just getting a feeling for how the queuing system works.

As an example, try the following:

% srun hostname

This will run the hostname command on one of the cluster nodes and return the results to your screen. You should see a different node name than the node you are logged into.

To run the same thing on multiple different nodes in the cluster, try the following:

% srun -n 8 hostname

The “-n 8” argument indicates that srun should run the indicated program – in this case, hostname, as 8 tasks.

This is different than using -N (capital N):

% srun -N 8 hostname

The “-N 8” argument indicates that srun should run the indicated program on 8 nodes.

The nodes in our cluster may have multiple CPU cores per nodes – so running 8 tasks may require only 2 or 4 nodes. However, specifying only four nodes will run 1 process per node (these processes can be multithreaded, however). Inspect the output from running the two version (-n 8 vs -N 8) and see if the output is consistent with two tasks per node.

--------------------------------------------------------------------------
[[56201,1],1]: A high-performance Open MPI point-to-point messaging module
was unable to find any relevant network interfaces:

Module: OpenFabrics (openib)
  Host: ip-10-11-15-142

Another transport will be used instead, although this may result in
lower performance.
--------------------------------------------------------------------------

Try running srun with hostname a few times. Experiment with various numbers of nodes. The fleet will currently only expand to 32 compute nodes, so if you request more nodes than that your job will not run (or may be rejected by srun).

sbatch¶

In the last assignment we looked at basic use of srun and sbatch and ran the following:

$ sbatch hello_script.bash

Look at the contents of the output file from this run. How many instances of the script were run (how many “hello” messages were printed)?

We can implement multiple instances with sbatch:

$ sbatch -N 4 hello_script.bash

How many instances were run?

This is different than what you might have expected. Run the script with srun

% srun -N 4 hello_script.bash

How many “hello” messages were printed?

The difference between the sbatch and the srun command is that sbatch just puts your job into the queuing system and requests resources to potentially run a parallel job (the -N option). But it doesn’t run the job in parallel. Running a job in parallel is the performed by srun (or, from a job in the queue, with mpirun).

Consider the following script (hello_script_p.bash), slightly modified from hello_script.bash:

#!/bin/bash
echo "Hello from Slurm"
srun hostname
echo "Goodbye"

Here, we run hostname with srun – note that we have not specified how many nodes we want. Launch a few trials of the second script

% sbatch hello_script_p.bash
% sbatch -N 2 hello_script_p.bash
% sbatch -N 4 hello_script_p.bash

Look at the .out files for each job. How many instances of hostname were run in each case? How many “hello” messages were printed?

We can also use mpirun to launch parallel jobs (and is what we will need to run our MPI jobs):

#!/bin/bash
echo "Hello from Slurm"
mpirun hostname
echo "Goodbye"

(This script is hello_script_m.bash)

Again, run this script a few different ways and see if the results are what you expect.

% sbatch hello_script_m.bash
% sbatch -N 2 hello_script_m.bash
% sbatch -N 4 hello_script_m.bash

sbatch and srun accept the same command-line arguments – so, again, tehre is a difference between -n and -N.

Autoscaling and Ganglia¶

A final note on using the AWS cluster. It is currently configured to “autoscale” – that is, when there is not a large demand for nodes it shuts them down (maintaining a minimum of 4 up at all times). If you submit a job with a value of -N that requires more nodes than are currently available, AWS will spin some more up. However, this will take a few minutes to complete. Your job will wait in the queue in the meantime.

If you want to see what the cluster is doing (how many nodes are up, how loaded they are, etc), you can visit the page http://IP/ganglia/ , where IP is the IP address of the head node. This shows a number of statistics from the cluster, monitored in real time.

Compiling and Running an MPI Program¶

As we have emphasized in lecture, MPI programs are not special programs in any way, other than that they make calls to functions in an MPI software library. An MPI program is a plain old sequential program, regardless of the calls to MPI functions. (This is in contrast to, say, an OpenMP program where there are special directives for compilation into a parallel program.)

Given that, we can compile an MPI program into an executable with the same C++ complier that we would use for a sequential program without MPI function calls. However, using a third-party library (not part of the standard library) requires specifying the location of include files, the location of library archives, which archives to link to, and so forth. Since MPI itself only specifies its API, implementations vary in terms of where their headers and archives are kept. Moreover, it is often the case that a parallel programmer may want to experiment with different MPI implementations with the same program.

It is of course possible (though not always straightforward) to determine the various directories and archives needed for compiling an MPI program. Most (if not all) implementations make this process transparent to the user by providing a “wrapper compiler” – a script (or perhaps a C/C++ program) that properly establishes the compilation environment for its associated MPI implementation. Using the wrapper compiler, as user does not have to specify the MPI environment. The wrapper compiler is simply invoked as any other compiler – and the MPI-associated bits are handled automagically.

The wrapper compilers associated with most MPI implementations are named along the lines of “mpicc” of “mpic++” for compiling C and C++, respectively. These take all of the same options as their underlying (wrapped) compiler. The default wrapped compiler is specified when the MPI implementation is built, but it can usually be over-ridden with appropriate command-line and/or environment settings. For MPI assignments in this course, the compiler is named mpic++ – the same name is used by both MPICH and Open MPI.

$ mpic++ hello.cpp

$ make hello.exe

$ srun -n 4 hello.exe

$ sbatch -n 4 hello.bash

if (argc >= 2) rounds = std::stol(argv[1]);

Ring¶

In the MPI ping-pong example we sent a token back and forth between two processes, which, while actually quite amazing, is still only limited to two processes.

Complete the code in ring.cpp so that instead of simply sending a token from process 0 to 1 and then back again, the token is sent from process 0 to 1, then from 1 to 2, then from 2 to 3 – etc until it comes back to 0. As with the ping pong example above, increment the value of the token every time it is passed.

Hint: The strategy is that you want to receive from myrank-1 and send to myrank+1. It is important to note that this will break down if myrank is zero or if it is mysize-1 because there is no rank -1 or mysize. These cases could be addressed explicitly. (Or, since with everything else in MPI, there may be capabilities in the library for handling this situation, since it is not uncommon.)

A completed ring.cpp program.

Norm¶

We couldn’t visit another high-performance computing paradigm without using it to compute the Euclidean norm of a vector. Computing the Euclidean norm with MPI has two parts, and the difficult part isn’t the one you would expect.

Solving Laplace’s Equation¶

The files that you will edit for this part of the assignment are in the grid subdirectory. In the include subdirectory you will find the following files:

• Grid.hpp implements a rectangular grid

• Stencil.hpp iplements a stencil operator (in particular, it implements the operator* between a stencil and a grid)

• ir.hpp that implements iterative refinement as described above

• jacobi.hpp that implements the Jacobi iteration to solve Laplace’s equation

• cg.hpp that implements the conjugate gradient algorithm

Your assignment is to parallelize – with MPI – these algorithms.

Your work will be carried out in the following files:

1. mpiStencil.hpp: contains implementation of the operator* function. This is basically the same as the sequential version except there is a call to update_halo(). You will need to fill in the body of this function. It takes a Grid x as an argument and should update the boundary values of x with the appropriate values from its neighbors. This problem is the focus of this assignment.

2. Final.hpp: contains an implementation of the ir(…) function. This is an overload of the ir(…) function in ir.cpp – it just takes an mpiStencil rather than a plain Stencil. In looking at the different functions called in the sequential version of ir() – what needs to change here to make it a parallel ir()? We have discussed everything you need for this in lecture. You may want to change the names of some functions that are called here but you probably do not need any MPI code in there directly. There is a Big Hint in the function that is defined at the top of Final.hpp (and that you will be finishing in the next question).

3. Final.hpp contains an implementation of the mpiDot() function.

   Think very carefully about this one. In this operation the distinction between real boundaries in your problem and the “pseudo-boundaries” realized by the ghost cells does become important.

4. 583 only Repeat 2 but for the cg() algorithm also included in Final.cpp.

There is a single driver program in the grid subdirectory to run the sequential and parallel versions of ir, jacobi, and cg.

This program takes a number of arguments:

• one of -j, -i, or -c to select jacobi, ir, or cg solver, respectively

• -x and -y specify the grid dimensions for the discretized Laplace’s equation

• -m specifies the maximum number of iterations to run the algorithm

• -t specifies the convergence tolerance for the iterations

• -d turns on debugging output

• -p specifies that the program should be run in parallel

You can run the program sequentially on the front end machine (the one you logged into). It will print the norm of the residual at each iteration of the algorithm.

The output from the sequential algorithm should match the corresponding parallel algorithm. E.g.., running parallel ir should give the same results as running sequential ir.

For running the parallel version of the program with sbatch, use the grid.bash script, e.g.,

$ sbatch --ntasks 8 -N 4 grid.bash -m 128

The grid.bash script automatically passes the -p option to grid.exe so you don’t need to specify that.