Programming Lab #4
Cardiac Electrophysiology Simulation with SSE Vector Intrinsics

Changelog

Contents

• Cardiac Electrophysiology Simulation
• The starter code
• The Assignment
• Turning in your assignment
• Grading

Overview of the Assignment

In this assignment you'll implement the heart electrophysiology simulator discussed in class using the SSE intrinsics and run on the Bang cluster. Using your implementation, you'll conduct experiments on the Bang cluster to assess performance. The assignment includes a brief report, and the requirements for your report are described at the end of this writeup.

‣ ‣ ‣ Required solution files: solve.cpp

This means that the files that contain your solution to the assignment must be turned in by the deadline. The turnin procedure is described below in the section "Turning in your assignment".

The procedure for obtaining the starter code is similar to the previous assignments, and can be found by visiting the following web page:

A4 Git Instructions

In doing this assignment, it's your responsibility to understand the course rules for integrity of scholarship.

The assignment includes a report, and the requirements for this report are described below.

Cardiac Electrophysiology Simulation

Simulations play an important role in science, medicine, and engineering. A cardiac electrophysiology simulator like the one we are using in this assignment can be used for clinical diagnostic and therapeutic purposes. To reduce computational costs, simplifying assumptions are often made to make the simulations feasible. In our clinical example, time may be of the essence in determining an appropriate life-saving treatment.

Cell simulators entail solving a coupled set of a equations: a system of Ordinary Differential Equations (ODEs) together with a Partial Differential Equation (PDE). If you are interested in learning more about cardiac modeling peruse www.labheart.org or thevirtualheart.org or visit the web page of the UCSD Cardiac Mechanics Research Group.

Our simulator models the propagation of electrical signals in the heart, and it incorporates a cell model describing the kinetics of the membrane of a single cell. There can be different cell models, of varying degrees of complexity. We will use a model known as the Aliev-Panfilov model, that maintains 2 state variables, and solves one PDE (The PDE couples multiple cells into a system). This simple model can account for complex behavior such as how spiral waves break up and form elaborate patterns. Spiral waves can lead to life threatening situations such as ventricular fibrillation.

Our simulator models electrical signals in an idealized biological system in which voltages vary at discrete points in time, called timesteps, at the discrete positions of a mesh. The simulator is an iterative procedure that updates the simulation state variables (on the mesh) many times over a sequence of discrete timesteps. A simulation will run for some amount of simulated time, which is a real number. (For our purposes, the units of time are unimportant.) The timesteps subdivide the simulation time line that starts at t=0 and ends at some specified time t=T. The unit of subdivision is called Δt, or the time discretization. It is important to keep in mind that when we use the term "time," we need to specify whether we are talking about real-valued simulation time, or the discrete timestep. We will use the fully qualified name to avoid ambiguity. By convention, we number the discrete timesteps with integers. Thus, if we say that we are at timestep k of the simulation, then the current simulation time is kΔt.

To discretize space, we'll use a uniformly spaced two-dimensional mesh, shown in the left side of the figure (Irregular meshes are also possible, but are more difficult to parallelize).

At each time step, the simulator updates the voltage according to nearest neighboring positions in space and time. This is done first in space, and next in time. Nearest neighbors in space are defined on the 4 Manhattan directions only: north, south, east and west, as shown on the right side of the above figure.

In general, the finer the mesh — and hence the more numerous the points — the more accurate the solution, but at an increased computational cost. In a similar fashion, the smaller the timestep, the more accurate the solution, but at a higher computational cost. Thus, if h is the mesh spacing, Δt is the time discretization, then the running time is proportional to (Δt×h²)^-1. However, to simplify our performance studies, we will run our simulations for a given number of iterations rather than a given amount of simulated time.

The simulator algorithm

The simulator maintains two state variables characterizing the cell electrophysiology model: excitation (E) and recovery (R). Since we are using the method of finite differences to solve the problem, we discretize the variables E and R by considering the values only at a regularly spaced set of discrete points. We represent the variables as 2D arrays: E[][] and R[][], respectively. The PDE solver also maintains the voltage at the previous timestep, E_prev.

Here is the formula for solving the PDE, where the superscripts t and t-1 refer to the current and previous timestep, respectively, and the constant α is defined in the simulator (see apf.cpp and solve.cpp):

Et[i,j] =  Et-1[i,j] +  α*(Et-1[i+1,j] + Et-1[i-1,j] + Et-1[i,j+1] + Et-1[i,j-1] - 4*Et-1[i,j])

Here are the formulas for solving the ODE:

E = E - Δt*(k*E * (E-a) * (E-1)+E * R);
R = R + Δt*(ε+ μ1*R / ( E+μ2)) * (-R-k*E * (E-b-1));

The variables k, a, b, ε, μ₁ and μ₂ are constant floating point scalars, and are defined in the simulator. References to E and R correspond to whole arrays. Expressions involving whole arrays (that is, array names without subscripts) are pointwise operations on corresponding array elements.

We interpret a pointwise whole array operation as follows. If A, B and C are 2D arrays of size N × N, then A = B + C is equivalent to the following loop nest:

for (i = 0; i < N; i++)
    for (j = 0; j < N; j++)
        A[i][j] = B[i][j] + C[i][j];

When an array is combined with a constant, e.g. (E-a), the constant is extended to cover each element of the array, i.e. we compute the E[i][j] - a for all i and j.

To organize our work, we observe a number of bookkeeping conventions. We divide each axis into N-1 subintervals using N points positioned uniformly along the line at positions x=i×h, i=0, 1, 2, ..., N-1. We also pad the boundaries with a cell on each side, in order to properly handle updates next to a boundary. (This padding will simplify the code and you'll see comments about it in solve.cpp.) This is shown in below, for 1 axis of the mesh drawn above, where the mesh spacing h = 1/N. The points labeled 'o' correspond to the physical boundary, the points labeled 'P' are padding.

i=0   i=1    i=2                    i=N-1   i=N   i=N+1
P------o------|------| ...... |------|------|------o------P
       ^                                           ^
       |                                           |
      x=0.0  x=h    x=2h                 x=1.0-h   x=1.0

Consider the case where we set N=64. We allocate an array of size 67 × 67, but we solve the problem on a central 65 × 65 submesh. Thus, our loops range from 1 to N+1, that is, 1 to 65. Points 0 and 66 correspond to buffer points.

It follows that when we subdivide the data to partition it for parallel execution, we are actually subdividing N+1 rather than N.

The provided code

In addition to the starter code, we also provide a simple example of how to use SSE intrinsics, which is found on Bang in $PUB/Examples/SSE/Vec

All of your changes will made to just one file:

solve.cpp

To build your code, use the provided Makefile, as it sets up the appropriate flags for the compiler via the "arch" file located in $(PUB)/Arch/arch.gnu-c++11.generic. The Makefile produces a target called apf. We will build your code from the provided Makefile (with vec=1 asserted) so don't depend on changes to this file. In particular, don't change the name of make target nor add additional code files. Moreover, to ensure that your program is correctly graded, do not change any of the following files:

apf.cpp  apf.h  cmdLine.cpp  Plotting.h   Plotting.cpp   Report.cpp   Timer.cpp

The code and Makefile have been set up to enable you to compile with or without the flags needed to support SSE vector intrinsics. The default is to compile without the intrinsics. To compile with intrinsics, specify following make command line argument to Make:

vec=1 

This setting defines the SSE_VEC macro which you can use to cause certain code to be conditionally compiled. To take advantage of this feature in any code you write, follow the convention used in the provided starter code. (There is a message printed out that tell you when you are building code with SSE intrinsics).

The simulator has various options, controlled via command line arguments (parsed in in cmdLine.cpp), with default values set in apf.cpp:

Thus, the following command line will conduct the simulation on a 255 × 255 box, run for 1000 iterations, plot the excitation variable every 20 iterations, and print out summary values of the excitation variable every 10 iterations.

    ./apf -n 255 -i 1000 -s 10 -p 20

If we leave off the -s and -p flags, then the output is as follows. The last line of output, beginning with the string "@ " outputs performance data along with relevant parameters, that is suitable for automated extraction from a sequence of runs.

[Simulation begins] Sat May  9 09:02:54 2015

dt= 0.040555, # iters = 1000
m x n = 255 x 255

[ *** Simulation ends at] Sat May  9 09:02:56 2015

 iteration 1000
Max norm: 9.928313e-01, L2norm: 7.007398e-01
Running Time: 0.933977 sec.
GFlop rate: 1.96

    M x N   [#iter]   {T_p, Gflops}   Linf, L2
@> 255 255     1000   0.934 1.965      9.92831e-01   7.00740e-01

The plotting capability, controlled via the -p flag, uses gnuplot, and may be useful in debugging. Do not enable plotting when you are measuring running times for your performance studies. The plotter capability only works on a single core and is discussed in more detail later in this writeup.

You'll notice that the provided code uses a different naming convention for the bounds of the computational box. In particular, while we use the single letter N in this document, the code uses different bounds, for the x and y axes, in terms of variables m and n. This convention enables the code to accommodate non-square local problem sizes that arise when the data has been partitioned for parallel execution. The computational boxes have been padded by one cell in each dimension. This change is needed to handle boundary conditions, which are described in solve()

Also note the conventions used to index the 2D solution arrays: the indices are swapped, with the i index being the most rapidly varying index, that is, indexing along the last dimension of the array. Thus, array elements in the same row have different i indices and elements i-1, i and i+1 are adjacent in memory.

The Assignment

Vectorizing the code

You'll use the SSE intrinsics discussed in class to nearly double the speed of your code. Note that, the compiler will attempt to vectorize certain code and will, as the optimization is enabled by -O3. However, the compiler will not succeed in vectorizing the most computationally intensive parts of the application and this is where your hand optimizations will make a significant difference in performance.

Not all code is amenable to optimization and code that consumes little of the application's running time won't make a significant difference even if it can take advantage of vectorization. The code used to to set up the "ghost cells" is case an point. Your efforts will be much more rewarding if you optimize the PDE and ODE code. Vectorize these two parts on at a time. We have set up a simple vectorization example to help get you started. However, in the example code the stencils are one dimensional and there are no alignment issues. In the PDE and ODE loops, you'll have two dimensional stencils and array strides are not integral powers of two in length.

Assignment Checklist

1. Vectorize your code using the SSE intrinsics. to obtain a performance improvement of about 60% for an input size of N=255. As you increase N, this performance improvement will shrink. We will test this case to assess, in part, performance of your code.
   Assume that N is an odd number.


2. Try to make the code go as fast as you can. There may be other ways of improving performance of your code. We will award 10% extra credit if you can run 15% or more faster than the provided reference binary.


3. Plot the performance benefit of your code compared to the provided starter code, over the range of N=127 to N=1023. Adjust the number of iterations so that your runs last just a few seconds. For N=1023, a reasonable number of iterations (-i) is about 20, for N=127, several thousand. The work scales quadratically with N).
   
4. Report any performance anomalies or outstanding issues in your code. We will provide some partial credit for test cases your code doesn't pass, if we receive a carefully documented description. Negative results can be a valuable means of learning, so we are giving you an incentive to explore this venue.

In the spirit of friendly competition, post your results to the A4 Piazza Forum.

Testing your code

Your code will be tested for correctness and performance against various combinations of the N and niter. N will always be an odd number; the provided reference will get a Segmentation Fault on even numbers.

Test your code by comparing against the provided implementation. You can catch obvious errors using the graphical output (see below), but a more precise test is to compare two summary quantities reported by the simulator: the L₂ norm and the L_∞ norm, the maximum value (magnitude) of the excitation variable. The former value is obtained by summing the squares of all solution values on the mesh, dividing by the number of points in the mesh (called normalization), and taking the square root of the result. For the latter measure, we take the "absolute max," that is the maximum of the absolute value of the solution, taken over all points. While there may be slight variations in the last digits of the reported quantities, variations in the first few digits indicate an error.

Verify that results are independent of the number of whether or not the code was vectorized, for various (legal) values of N.

Plotting capability

As mentioned previously, the code includes with a plotting capability to observe the simulation dynamics which can be helpful when debugging. The plotter (Plotting.cpp) interfaces with gnuplot, which has been installed on bang (and is freeware). Note that two entries are provided, one to the provided multidimensional array implementation, and another to support flattened arrays.

Remote visual display can be slow, so either work from a campus machine, or increase the frequency to at least 20 (-p 20). The summary statistics may be more useful than plotting in some cases. If a display window does not appear on your screen, then your ssh display settings are not set correctly. If you are connecting to the remote machine from a Linux platform using ssh, use the -X option:

% ssh -X yourlogin@ccom-bang-login.ucsd.edu

If you are connecting from a Windows platform, and are using a 'secure shell client' AND 'exceed' (or xwin32, or cygwin) to log into a remote machine you must ensure that 'X11 tunneling' is set. To do this, you must have installed Cygwin /X with X11.

The simulations will produce a spiral pattern that evolves with time. However, if you set the number of iterations (-i) to a sufficiently large value, the organized spiral pattern will break up. This indicates non-physical behavior, and the length of time that the simulator produces physically meaningful results depends on the mesh size N.

Be sure that, while collecting performance data, that you have not enabled plotting (-p) or statistics reporting (-s).

Turning in your assignment

1. Submission via GitLab

When you have completed the assignment, you must submit it. You will do this by pushing your changes to GitLab. Follow these steps:

1. Run a ’git status’ to see what files have been changed.
2. Use ’git add <files>’ to stage certain files for a commit, use ’git add .’ to add all files in the working directory to the commit (hopefully you will not be doing this at this late stage—all files should already be in your repository!)
3. Use ’git commit -m "FINAL SUBMISSION"’ to commit the changes to your repository and indicate that you are done and this is your final submission.
4. To push your changes onto your repository on GitLab run ’git push origin master’ Don’t forget this push! You must do this to get credit.
5. As an added safety check, you should log on to GitLab on the web and verify that your push went through.

Some of the finer points about this process:

1. Watch out for (and resolve) merge conflicts before your final push!!
   If you have been developing on multiple computers and not syncing with GitLab between switching computers, you will likely get merge conflicts. This happens because your code from one repo is out of sync with your code from another and GitLab doesn't know how to merge the files. You will get these conflicts "at the last minute," because git must sync the code in the repo you're trying to push with the code in the repo you're trying to push to. MERGE CONFLICTS WILL PREVENT YOUR CODE FROM COMPILING. So, even if your code was working perfectly just before you attempted your final push, conflicts will break your code.
   
   It is your responsibility to recognize when you have conflicts (it's not hard--you get a big warning), and to fix them. Please leave yourself enough time to do this.


2. If you are submitting with a partner, just push to the joint repository on behalf of both partners. That’s your joint submission.


3. Use the commit message FINAL SUBMISSION for your final submission in only one of your individual and partner repositories. You shouldn't be using both anyway, but if you happen to have pushed changes to both, just make sure only one contains the FINAL SUBMISSION commit message.


4. If you decide you want to resubmit your assignment before the deadline, just do another commit with the message FINAL SUBMISSION (and don’t forget to push to GitLab!!). We will grade the last commit before the deadline that is tagged as the FINAL SUBMISSION.


5. If you submit the assignment within 24 hours after the deadline (i.e. if there is a commit tagged "FINAL SUBMISSION") you will be charged a 10% lateness penalty. . If it is more than 24 but within 48 hours, you will be charged a 20% penalty. After 48 hours, we’ll roll back to the last commit tagged as "FINAL SUBMISSION" that was submitted before the deadline.

2. Required files in the submission

1. A PDF document that presents all requested performance data. This file must have the following format. The stem part of the file name (before the .pdf) should begin with "A4-report" and end with you repository name, e.g. A4-report-baden-powell.pdf or, if working alone: A4-report-baden.pdf.

   You may provide data in tabular form only. If you provide a plot, you must still provide the tabular data as sometimes plots can be hard to read. Be sure to provide a caption for each plot or table that describes the experiment you conducted, e.g. "strong scaling."

   Plot your data carefully, plotting only what's necessary to get your points across, and choosing log and linear plots carefully. You must also include any plotted data in tabular form.

   Negative results are also valid and especially valued if you can document them, with the underlying causes, or best guesses if the causes are not clear. So if results weren't as you expected, let us know why and we may be able to give you some partial credit!

   Cite any written works or software as appropriate in the text of the report. Citations should include any software you used that was written by someone else, or that was written by you for purposes other than this class.

2. Source code. Even though we will be substituting fresh copies of all but one source module, include all code files your worked with, in case we need to unravel a problem with your turnin. For this reason, also include your Makefile.
3. 2 Required text documents

   Your turnin must also fill in two important documents, which accompany the provided code.

   1. A completed self-evaluation file named teameval.txt discussing your work; if in a team this form address the division of labor and other aspects of how you worked together as a team.
   2. A file named DECLARATION, that states that the work was your own, and that you did not violate the Academic Integrity policies in effect for the course. Modify the DECLARATION by filling in your correct login ID(s) and the date only. Don't change anything else.
      

   These documents are required and we will not report a grade for your assignment without them. A 10% penalty will be assessed in cases of missing forms, empty forms, or forms not properly filled out.

Grading

Your grade will depend on 3 factors: performance, correctness and the depth and your explanations of observed performance and how the design of your code influenced performance. In determining your grade for the assignment, explanations and insight will be emphasized more heavily than performance itself, provided you were able to get your code to scale.

We will test your code on bang. Be sure to do your final testing on that machine, to avoid any surprises at turnin time. Our test script will build your code against fresh copies of the provided source files mentioned above. Then, it will run your code against various inputs to measure your code's performance, and test for correctness.

The grade for your comments will be based on clear explanations in your source code and GitLab log entries. When you commit files with GitLab, be sure to add brief, but descriptive comments, suitable for understanding your development process. If you aren't sure, be sure to ask us.

Frequently asked question, where classmates answer questions

SSE Documentation

• Intrinsics Guide This is a handy "cheat sheet," a complete reference to all SIMD intrinsics on Intel Architectures which you can filter by type, SSE version number. Bang supports SSE, SSE2 and SSE3 instruction sets but most of what you need is in SSE2.
• How to Write Fast Code SIMD Vectorization, Markus Püschel, CMU (2008). This presentation gives you a nice overview.
• SSE2 Preview, and MMX Primer, by Stefano Tommesani. (Many writeups on SSE accessible from this page.
• Intel also publishes complete reference guides and manuals. Save a tree and preview these electronically.
  • Intel 64 and IA-32 Architectures Software Developer Manuals
  • Intel Architecture Instruction Set Extensions Programming Reference, October 2014. This manual specifies the details behavior of machine instructions.
• GCC documentation - We are using GCC 4.8 on Bang, but this page documents older versions too. Look at the section entitled "Intel 386 and AMD x86-64 Options," especially the optimizations concerning SSE.
• An Introduction to GCC Compiler Intrinsics in Vector Processing, by George Koharchik and Kathy Jones, Linux Journal, Sep 21, 2012.

[Sat May 9 10:17:55 PDT 2015]