Problem Set #8¶

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 the Stencil class 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¶

For this assignment (and for the rest of the quarter), we will be using the UW Hyak compute cluster.

As a reminder, the Hyak cluster is a shared compute resource, which sharing includes the compute resources themselves, as well as the file system and the front-end (login) node. 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 this assignment (no mining bitcoins or other activities for which a cluster might seem useful),

• Do not run anything compute-intensive on the head node (all compute-intensive jobs should be batch scheduled),

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

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

• Use the –time=5:00 argument to all jobs you submit.

Investigating the Cluster¶

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 -A niac --tasks=8 hostname

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

Based on what is printed as a result, are the tasks run on one host or multiple hosts?

Specifying the number of tasks is different from specifying the number of nodes:

% srun --nodes=8 hostname

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

Based on what is printed as a result, are the tasks run on one host or multiple hosts?

How many tasks are launched, and how many nodes are used, when you combine tasks and nodes?

% srun --nodes=8 --tasks=2 hostname

The nodes in Hyak have multiple CPU cores per nodes (how many?) and we can specify tasks per node to make best use of the cores (and other resources). Note also that we can still specify cpus per task for multi-threaded tasks. (You might try this with ompi_info.exe):

% srun --nodes=4 --tasks=2 --cpus-per-task=8 ./ompi_info.exe

Try running srun with hostname a few times. Experiment with various numbers of nodes, tasks, cpus.

$ srun -p ckpt --nodes=8 --tasks=4 --cpus-per-task=4 ./ompi_info.exe

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 --nodes=4 hello_script.bash

How many instances were run?

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

% srun --nodes=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 --nodes 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 --nodes=2 hello_script_p.bash
% sbatch --nodes=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 --nodes=2 hello_script_m.bash
% sbatch --nodes=4 hello_script_m.bash

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 “mpicxx” 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 mpicxx.

$ mpicxx hello.cpp

$ make hello.exe

$ srun --mpi=pmi2 --tasks=4 hello.exe
$ srun --mpi=pmi2 --nodes=4 hello.exe
$ srun --mpi=pmi2 --nodes=4 --tasks=2 hello.exe

$ sbatch --tasks=4 hello.bash
$ sbatch --nodes=4 hello.bash
$ sbatch --nodes=4 --tasks=2 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. Copy and paste relevant code lines that contain your edits to your report. Provide comments in the code near your edits to explain your approach.

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.)

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. mpiMath.hpp: contains the implementation of jacobi, which implements the Jacobi iteration. To parallelize this, we need to do two things. First, we need to correctly compute the global value of rho – which needs to be available on all ranks after its computation. Second, we need to properly exchange boundary values for the x Grid. This problem is the focus of this assignment.

2. mpiStencil.hpp: contains the implementation of the mult. This function also needs to exchange boundary values of x. You can use the same scheme as you used in jacobi. For extra credit, try a different scheme.

3. mpiMath.hpp contains an implementation of the mpi_dot() 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. mpiMath.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 mpiMath.hpp.

5. 583 only Repeat 4 but for the cg() algorithm also included in mpiMath.hpp.

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 --nodes=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.

Note that as with many of our benchmarking programs, we are assuming things divide nicely by powers of two (that is, sequential equivalence will only be guaranteed and/or expected when running on tasks and nodes that are powers of two).

$ ./grid.exe -j
# elapsed time [seqgrid]: 1 ms
# 128 iterations
# ||r|| = 0.0490717

$ sbatch grid.bash -j

Starting parallel grid
# elapsed time [mpigrid]: 5 ms
# 128 iterations
# ||r|| = 0.0490717
Done with parallel grid

$ sbatch --nodes=2 grid.bash -j

$ sbatch --nodes=2 grid.bash -j -d