Parallelization

In computer programming parallelization is a way to speed a program up by spreading the calculations over more compute resources. The idea itself is rather old. For example, Amdahl published his paper suggesting his famous law about the most benefit you can get from parallelization in 1967[1], and Kuck et al. published a paper analyzing Fortran IV codes regarding the potential for parallel execution in 1972[2]. Nevertheless, Moore’s law[3], from 1965, meant that for a considerable amount of time one could just wait for computers to become faster rather than spend the effort to parallelize programs.

Still in the 1980s efforts towards parallel programming started that resulted in the Parallel Virtual Machine (PVM)[4] and Message Passing Interface (MPI)[5] approaches in the early 1990s. Today MPI still is the leading technology for distributed data parallel computing.

Another realization early on was that compilers can do parallelization for you. This realization led to efforts such as High Performance Fortran[6]. In principle compilers can even generate distributed data parallel codes. In practice having compilers write good MPI codes is highly non-trivial, and even Virtual Shared Memory models tend to have performance issues due to communication latencies[7]. Other problems include that static code analysis cannot resolve which loops are good targets for parallelization, and data dependencies can be difficult to analyze[8]. In addition successful porting efforts to parallel platforms typically require significant code rewrites to express the parallelism on a much higher level in the code design. The need to bring programmer knowledge into the code to direct the compiler’s parallelization efforts led to directive based approaches such as OpenMP[9]. The latest developments with OpenMP focus on offloading computation to accelerator devices such as GPUs[10].

Introducing the matrix-matrix multiplication

Given the more than 50 year long history of parallel computing it is a valid question to ask “where to start?”. This page illustrates currently leading approaches based on a matrix-matrix multiplication example. Matrix-matrix multiplication is a good example because:

  • The basic algorithm is simple and easily understood

  • It is an important step in many engineering applications

  • It requires a substantial amount of computation so parallelization is a meaningful thing to do

  • The algorithm is not completely trivial so some generally relevant issues will come up

The implementations shown here are available as mini apps[11]. The codes are written in Fortran, mainly because Fortran supports convenient ways of using multi-dimensional arrays. There is no fundamental reason why these examples could not be expressed in C or C++ (Python is a different story because the interpreter model imposes additional restrictions, but MPI is available in through mpi4py[12].)

Finally, the codes provided here are truly examples. In a real application you should call an optimized library routine for matrix-matrix multiplication. The BLAS[13] collection of routines is a good place to find suitable routines. Here the matrix-matrix multiplication is explicitly coded so that the implementation can be modified.

In terms of matrices the calculation can be expressed as C := C + A * B. The corresponding code is

do jj = 1, nsize
  do kk = 1, nsize
    do ii = 1, nsize
      c(ii,jj) = c(ii,jj) + a(ii,kk)*b(kk,jj)
    enddo
  enddo
enddo

There are a few things to note here:

  1. First of all, the three matrices are each \(O(N^2)\) in size, i.e. if the dimension of the matrix is nsize then you need nsize * nsize numbers to store it.

  2. Second, the algorithm has a computational complexity of \(O(N^3)\). I.e. it takes nsize * nsize * nsize operations to perform this calculation.

  3. Third , the loops are put in a particular order because of the way Fortran lays a matrix out in memory. In Fortran the distance in memory between c(ii,jj) and c(ii+1,jj) is 1, whereas the distance between c(ii,jj) and c(ii,jj+1) is nsize. The loops as shown here ensure the algorithm takes the smallest steps in memory which is typically beneficial for performance as it maximizes cache reuse. Note that in C/C++ the memory layout is the other way around, and you would want the inner loop to run over jj and the outer loop over ii.

Distributed Data Parallelization with MPI

In the collection of Mini Apps the case called mpi-mxm demonstrates a distributed data parallelization of the matrix-matrix multiplication using MPI. Distributed data refers to the spreading of data across all available compute resources. This way if you have \(P\) processors then each processor only needs to hold about \(1/P\) of the data of the problem. At the same time you can bring the compute power of \(P\) processors to bear on the problem.

In order to split the problem over a potentially large number of processors it is not sufficient to partition the matrices over just 1 index. Doing that would limit the amount of processors the compute can be shared with to the number of elements along one dimension, i.e. nsize in our case example. If the matrices are split into blocks then, in principle, the compute could be shared over nsize * nsize processors. To keep the algorithm simple it is assumed that the number of processors is the square of an integer \(p\). If it isn’t a square then the remaining processors will simply not be used. In addition, if \(P = p^2\) we also assume that \(\mathrm{nsize} = n p\) where \(n\) is some whole number. This means that all matrices can be broken down into \(p^2\) blocks that are each \(n\)-by-\(n\) elements large.

Thus the matrices are partitioned as shown below. This example shows the distribution on a 4-by-4 processor grid (i.e. \(p=4\))

        block-beta
  columns 3
  block
    columns 1
    block
      C11["C<sub>0,0</sub>"] C21["C<sub>0,1</sub>"] C31["C<sub>0,2</sub>"] C41["C<sub>0,3</sub>"]
    end
    block
      C12["C<sub>1,0</sub>"] C22["C<sub>1,1</sub>"] C32["C<sub>1,2</sub>"] C42["C<sub>1,3</sub>"]
    end
    block
      C13["C<sub>2,0</sub>"] C23["C<sub>2,1</sub>"] C33["C<sub>2,2</sub>"] C43["C<sub>2,3</sub>"]
    end
    block
      C14["C<sub>3,0</sub>"] C24["C<sub>3,1</sub>"] C34["C<sub>3,2</sub>"] C44["C<sub>3,3</sub>"]
    end
  end

  block
    columns 1
    block
      A11["A<sub>0,0</sub>"] A21["A<sub>0,1</sub>"] A31["A<sub>0,2</sub>"] A41["A<sub>0,3</sub>"]
    end
    block
      A12["A<sub>1,0</sub>"] A22["A<sub>1,1</sub>"] A32["A<sub>1,2</sub>"] A42["A<sub>1,3</sub>"]
    end
    block
      A13["A<sub>2,0</sub>"] A23["A<sub>2,1</sub>"] A33["A<sub>2,2</sub>"] A43["A<sub>2,3</sub>"]
    end
    block
      A14["A<sub>3,0</sub>"] A24["A<sub>3,1</sub>"] A34["A<sub>3,2</sub>"] A44["A<sub>3,3</sub>"]
    end
  end

  block
    columns 1
    block
      B11["B<sub>0,0</sub>"] B21["B<sub>0,1</sub>"] B31["B<sub>0,2</sub>"] B41["B<sub>0,3</sub>"]
    end
    block
      B12["B<sub>1,0</sub>"] B22["B<sub>1,1</sub>"] B32["B<sub>1,2</sub>"] B42["B<sub>1,3</sub>"]
    end
    block
      B13["B<sub>2,0</sub>"] B23["B<sub>2,1</sub>"] B33["B<sub>2,2</sub>"] B43["B<sub>2,3</sub>"]
    end
    block
      B14["B<sub>3,0</sub>"] B24["B<sub>3,1</sub>"] B34["B<sub>3,2</sub>"] B44["B<sub>3,3</sub>"]
    end
  end

In algorithm we need to choose how to distribute the compute. For this I have chosen to do all the compute associated with a block of matrix C on the processor that holds that block. You could make different choices here, for example, do everything associated with the local block of A on the processor that holds it. To compute the local block of C we need pairs of blocks from rows of A and columns of B as shown below.

        block-beta
  columns 3
  block
    columns 1
    block
      C11["C<sub>0,0</sub>"] C21["C<sub>0,1</sub>"] C31["C<sub>0,2</sub>"] C41["C<sub>0,3</sub>"]
    end
    block
      C12["C<sub>1,0</sub>"] C22["C<sub>1,1</sub>"] C32["C<sub>1,2</sub>"] C42["C<sub>1,3</sub>"]
    end
    block
      C13["C<sub>2,0</sub>"] C23["C<sub>2,1</sub>"] C33["C<sub>2,2</sub>"] C43["C<sub>2,3</sub>"]
    end
    block
      C14["C<sub>3,0</sub>"] C24["C<sub>3,1</sub>"] C34["C<sub>3,2</sub>"] C44["C<sub>3,3</sub>"]
    end
  end

  block
    columns 1
    block
      A11["A<sub>0,0</sub>"] A21["A<sub>0,1</sub>"] A31["A<sub>0,2</sub>"] A41["A<sub>0,3</sub>"]
    end
    block
      A12["A<sub>1,0</sub>"] A22["A<sub>1,1</sub>"] A32["A<sub>1,2</sub>"] A42["A<sub>1,3</sub>"]
    end
    block
      A13["A<sub>2,0</sub>"] A23["A<sub>2,1</sub>"] A33["A<sub>2,2</sub>"] A43["A<sub>2,3</sub>"]
    end
    block
      A14["A<sub>3,0</sub>"] A24["A<sub>3,1</sub>"] A34["A<sub>3,2</sub>"] A44["A<sub>3,3</sub>"]
    end
  end

  block
    columns 1
    block
      B11["B<sub>0,0</sub>"] B21["B<sub>0,1</sub>"] B31["B<sub>0,2</sub>"] B41["B<sub>0,3</sub>"]
    end
    block
      B12["B<sub>1,0</sub>"] B22["B<sub>1,1</sub>"] B32["B<sub>1,2</sub>"] B42["B<sub>1,3</sub>"]
    end
    block
      B13["B<sub>2,0</sub>"] B23["B<sub>2,1</sub>"] B33["B<sub>2,2</sub>"] B43["B<sub>2,3</sub>"]
    end
    block
      B14["B<sub>3,0</sub>"] B24["B<sub>3,1</sub>"] B34["B<sub>3,2</sub>"] B44["B<sub>3,3</sub>"]
    end
  end

  style C23 stroke-width:4px
  style A13 stroke-width:4px
  style A23 stroke-width:4px
  style A33 stroke-width:4px
  style A43 stroke-width:4px
  style B21 stroke-width:4px
  style B22 stroke-width:4px
  style B23 stroke-width:4px
  style B24 stroke-width:4px

So the pseudo code becomes for a processor at coordinates \(p_{row}\), \(p_{col}\) on the processor grid with dimension \(p\)

MPI parallel matrix-matrix multiplication

allocate local A-block buffer A_buf
allocate local B-block buffer B_buf
do pp = 0, p-1
  if (pp == p_col) then
    A_buf = Ap_row,p_col
  end if
  MPI_Bcast(A_buf,communicator(p_row))
  if (pp == p_row) then
    B_buf = Bp_row,p_col
  end if
  MPI_Bcast(B_buf,communicator(p_col))
  do jj = 1, nsize_block
    do kk = 1, nsize_block
      do ii = 1, nsize_block
        C-block(ii,jj) = C-block(ii,jj) + A_buf(ii,kk) * B_buf(kk,jj)
      end do
    end do
  end do
end do
deallocate A_buf, B_buf

Where communicator(p_row) is the MPI communicator that includes all processors with the same value for p_row. I.e. it is a communicator for all processors in the row p_row. Similarly communicator(p_col) is the communicator for the column p_col.

Shared Memory Parallelization with OpenMP

Above one way was shown to parallelize matrix-matrix multiplication with MPI. While for some applications it is essential to be able to use compute resources that exceed the capabilities of a single node, for many applications that may not be required. If you want to speed your code up but don’t need more than a single node you can parallelize your code with multithreading through OpenMP[10]. A nice feature of OpenMP is that you can get the performance from multithreading without explicitly having to deal with threads. Instead OpenMP relies on directives that you put into the code that guide the compiler in generating parallel code for you. The directives are inserted in comments so that if you compile the code without OpenMP you will get a working serial program. So you don’t lose your original program by inserting OpenMP directives.

An additional advantage that follows from this is that you can parallelize your code piecewise. You can pick the most expensive part of your code, parallelize it, and then see what the next most expensive part is. In this context, mind that if you parallelize your code you might want to run larger problems, for which different parts may have a high cost compared to small problems.

One key consideration is that although the compiler generates the parallel code for you, you are still responsible that the code functions correctly. The most common problem is that OpenMP relies on shared memory. That means that by default all threads use the same variables. Of course if two threads change the same variable in an uncoordinated fashion this leads to unpredictable results. This is referred to as a race condition, i.e. running your program is like a race between processors and the results depend on which processor gets somewhere first.

OpenMP has many ways to manage access to variables to avoid race conditions. For example, a variable can be declared thread local so that only one thread can access it. Also all threads are numbered and you might use the thread number to access only some elements of an array. Alternatively, you can define critical regions, which are pieces of code that only one thread at a time can execute. We won’t go into all these options here but mind that there are multiple options and it pays to see what is available and what is suitable for your code.

In the introduction the matrix-matrix multiplication was presented with a three-fold nested loop ii, kk, jj where the ii-loop is the inner most loop. As written we can parallelize the code straightforwardly over the jj-loop.

!$omp parallel do private(kk,ii)
do jj = 1, nsize
  do kk = 1, nsize
    do ii = 1, nsize
      c(ii,jj) = c(ii,jj) + a(ii,kk)*b(kk,jj)
    enddo
  enddo
enddo
!$omp end parallel do

Note that in OpenMP the loop counter that is part of a work sharing construct is automatically made private. So we don’t need to say that jj is private. But kk and ii are not part of a worksharing construct so we need to explicitly say they are private.

Depending on the dimension of the matrix this provides rather limited opportunities. For example, if the dimension of the matrix is 100, then I can parallelize only over 100 threads. On amplitUDE that wouldn’t even exploit all the 112 cores of a single node. As in this case matrix \(c\) would be 100-by-100 in size, and every element can be calculated independently we should be able to do much better than that. OpenMP has a directive to parallelize over a combination of loops to help with this.

!$omp parallel do collapse(2), private(ii)
do jj = 1, nsize
  do kk = 1, nsize
    do ii = 1, nsize
      c(ii,jj) = c(ii,jj) + a(ii,kk)*b(kk,jj)
    enddo
  enddo
enddo
!$omp end parallel do

The collapse(2) directive tells the compiler to combine the jj- and kk-loop into one. I.e. the compiler will actually do:

!$omp parallel do private(jj,kk)
do jk = 1, nsize*nsize
  jj = (jk-1)/nsize + 1
  kk = mod(jk-1,nsize) + 1
  ...

and then parallelize the jk-loop. Note that in C/C++ this will be slightly different as there you count starting at 0 instead of 1. Also note that now jj and kk are no longer loop counters and so we need to make them private explicitly.

Unfortunately, the resulting code is wrong! When the kk-loop is split over multiple threads these threads change the same pieces of matrix \(c\) and a race condition has been introduced. In this case the problem is easily fixed by interchanging the ii- and kk-loops to get:

!$omp parallel do collapse(2), private(kk)
do jj = 1, nsize
  do ii = 1, nsize
    do kk = 1, nsize
      c(ii,jj) = c(ii,jj) + a(ii,kk)*b(kk,jj)
    enddo
  enddo
enddo
!$omp end parallel do

This code is available in openmp-mxm. Note that, as was stated in the introduction, this loop structure leads to less efficient memory access. Nevertheless, this loop structure is the most straightforward way to ensuring correctness. Therefore, paying a little loss in memory access efficiency is a worthwhile price for correctly getting more parallelism.

MPI + OpenMP Parallelization

We saw that we can use MPI to parallelize across nodes. A remaining issue is that in an MPI code the number of messages we need to send increases with the number MPI processes. For example, if I run the matrix-matrix multiplication MPI code on a single CPU then I don’t need to do any communication, but I have to do \(\mathrm{nsize}^3\) multiplications. If I ran the code on \(\mathrm{nsize}^2\) CPUs then every CPU would do \(\mathrm{nsize}\) multiplications, but also \(2*\mathrm{nsize}\) MPI broad casts. So, for a matrix-matrix multiplication, if I distribute larger tasks I have less of a communication overhead, but at the same time each task involves more compute and therefore takes longer to complete.

One way to resolve this issue is instead of running a single MPI process per CPU to run a single MPI process per node, for example. The time it takes to do the resulting large tasks can be reduced by using OpenMP to accelerate the MPI processes. As MPI and OpenMP parallelization are independent of eachother all we need to do is to put the two approaches together as shown below (see also mpi-openmp-mxm)

MPI+OpenMP parallel matrix-matrix multiplication

allocate local A-block buffer A_buf
allocate local B-block buffer B_buf
do pp = 0, p-1
  if (pp == p_col) then
    A_buf = Ap_row,p_col
  end if
  MPI_Bcast(A_buf,communicator(p_row))
  if (pp == p_row) then
    B_buf = Bp_row,p_col
  end if
  MPI_Bcast(B_buf,communicator(p_col))
  !$omp parallel do collapse(2) private(kk)
  do jj = 1, nsize_block
    do ii = 1, nsize_block
      do kk = 1, nsize_block
        C-block(ii,jj) = C-block(ii,jj) + A_buf(ii,kk) * B_buf(kk,jj)
      end do
    end do
  end do
  $!omp end parallel do
end do
deallocate A_buf, B_buf

One note worthy restriction on this approach is that MPI communicates between processes, it cannot communicate between threads. Therefore, you cannot place MPI calls inside an OpenMP parallel region! You will need to keep that in mind in your algorithm design.

MPI + OpenMP GPU off-loading

Sofar parallelization of matrix-matrix multiplications over CPUs has been discussed in various options. On the lastest high performance computing systems the CPUs unfortunately are the slowest resources available. Instead, the most compute power comes from accelerators such as GPUs. Hence it would be very beneficial if we could move most of the compute to the accelerator devices. One aspect to keep in mind though is that with current architectures the bandwidth between the host’s main memory and the device memory is relatively low. So moving data from the host to the device or from the device to the host is typically slow. Over time this is likely to become less important, in particular as improved device-to-device communications tend to reduce the need to communicate data with the host. At present, nevertheless, it is not worthwhile to move small tasks to the device. Conversely, when you have all the required data on the device already, it may also not be worthwhile to move small amounts of compute to the host, even if the compute is not suited to GPUs.

In recent times off loading compute to devices has been greatly facilitated by additions to the OpenMP standard. In order to understand these additions it helps to know a little bit more about the architecture of GPUs. Like any multi-core processor a GPU consists of a number of compute cores. On regular multi-core processors every core is a fully fletched CPU. The down side of this is that a lot of silicon and hence a lot of energy is spent on the control infrastructure that manages the scheduling and execution of instructions. On GPUs, this has been addressed by having multiple cores share their control infrastructure. On Nvidia GPUs 32 cores share their control infrastructure and this unit is referred to as a warp. On AMD GPUs there are 64 cores to a warp. Intel GPUs use wavefronts in a similar manner, except that you can control the number of cores per wavefront to some extent. The important restriction that comes from sharing control infrastructure is that all cores in a warp or wavefront must execute the same instructions. Therefore, the matrix-matrix multiplication example discussed here is an excellent candidate because for every pair of ii and jj the same dot-product (the loop over kk) is executed.

Algorithms with if-statements are less well suited. Consider the following code snippet:

do ii = 1, huge_number
  if (d(ii) < 0.0) then
    c(ii) = a(ii) * b(ii)
  else
    c(ii) = a(ii) + b(ii)
  end if
end do

Now dependent on the values of d some cores in a warp may have to perform a multiplication and others an addition. But because all cores in a warp have to do the same operation this situation has to be resolved in a special way. The way this done is that all cores will do both instructions. For example, all cores in the warp do the multiplication first and then the result on all cores where d(ii) >= 0.0 is thrown away. Then all cores do the addition and the result on all cores where d(ii) < 0.0 is thrown away. Clearly, this means that if your code contains if-statements that require that different cores to do different computations (also referred to as branch-divergence) this becomes very costly. The larger the number of branches, and the more compute there is to do in each branch, the higher the impact on performance on GPUs. This is something to keep in mind when selecting code sections to off-load to a GPU.

When it comes to off-loading compute to GPUs with OpenMP the directives do refer in some extent to the GPU architecture. Instead of the usual parallel do we have to indicate that we are targeting GPUs for the execution. In OpenMP the target construct is used for this purpose and its clauses can be used to manage the data migration to and from the device. In addition there is the teams construct to group threads into “teams” where each team executes on a warp. Put together using OpenMP to off-load a matrix-matrix multiplication to a single GPU looks like

!$omp target map(tofrom: c) map(to: a, b)
!$omp teams distribute parallel do collapse(2) private(kk)
do jj = 1, nsize
  do ii = 1, nsize
    do kk = 1, nsize
      c(ii,jj) = c(ii,jj) + a(ii,kk) * b(kk,jj)
    end do
  end do
end do
!$omp end teams distribute parallel do
!$omp end target

Here the map clauses are the main new thing. The tofrom clause indicates that c has to be moved from the host to the device, then it can be changed on the device, and the result needs to be moved back to the host. The other two matrices only need to be moved to the device, as indicated with the to clause but as they don’t change there is no need to move those matrices back.

If you need an intermediate data structure on the device you can create it on the device with alloc and destroy it when you’re done with delete. For example, a matrix transformation \(C = B^T A B\) can be rewritten as \(D = A B; C = B^T D\) to lower the computational complexity and off-loaded to the GPU as (see also openmp-offload-mxmxm

!$omp target enter data map(to: a, b) map(alloc: c, d)
!$omp target map(to: a, b, d)
!$omp teams distribute parallel do collapse(2) private(kk)
do jj = 1, nsize
  do ii = 1, nsize
    d(ii,jj) = 0.0
    do kk = 1, nsize
      d(ii,jj) = d(ii,jj) + a(ii,kk) * b(kk,jj)
    end do
  end do
end do
!$omp end teams distribute parallel do
!$omp end target
!$omp target map(to: b, d) map(from: c)
!$omp teams distribute parallel do collapse(2) private(kk)
do jj = 1, nsize
  do ii = 1, nsize
    c(ii,jj) = 0.0
    do kk = 1, nsize
      c(ii,jj) = c(ii,jj) + b(kk,ii) * d(kk,jj)
    end do
  end do
end do
!$omp end teams distribute parallel do
!$omp end target
!$omp target exit data map(delete: c, d)

Here target enter data and target exit data constructs have been added. These constructs allow to move data to and from the device outside of OpenMP parallel regions. For example, it is possible to have code to be executed on the host between target enter data and target, such as MPI calls to move data between compute nodes.

In the target enter data construct there is an map alloc clause that arranges memory on the GPU for matrices c and d. The following target construct has a map to clause that includes d. Clearly, that is strange. Why do we need to tell OpenMP to move d to the device when it is actually first calculated on the device? This is a consequence of the design fundamentals of OpenMP. OpenMP uses compiler directives embedded in comments. An OpenMP capable compiler can translate these directives and generate the code we want, but a compiler that is unaware of OpenMP should be able to ignore these directives and still generate a correctly working code. This requirement that a compiler that is not OpenMP aware or only supports part of an OpenMP standard still produces working code requires special rules. One such rule is that there has to be a map clause for every variable that is used in a device kernel. If no map clause is given then this is resolved by the compiler implicitly inserting a map tofrom clause. So if we would have written

!$omp target enter data map(to: a, b) map(alloc: c, d)
!$omp target map(to: a, b)

then the compiler would read this as if we had written

!$omp target enter data map(to: a, b) map(alloc: c, d)
!$omp target map(to: a, b) map(tofrom: d)

because tofrom is the safe choice with respect to program correctness if the compiler ignores part of OpenMP. But this would mean that the program would move the data from d back to the host at the end of the parallel do construct. As d is not needed on the host we don’t want the code to do that. Therefore we are forced to explicitly provide a map clause that says something about d. When we say map(to: d) then the code will check whether d is already on the device, which it is because we have just created it there, and then decide that no data transfer is needed. At the end of the parallel do no data transfer will happen either because we have said explicitly that we don’t want d to be moved to the host (there is no from part in the map clause related to d). So now, the data in d won’t be moved, which is what we wanted even if we had to say this in a bit of a strange way.

Putting the OpenMP off-loading and the MPI framework together we get (see also mpi-openmp-offload-mxm)

MPI+OpenMP parallel matrix-matrix multiplication

allocate local A-block buffer A_buf
allocate local B-block buffer B_buf
!$omp target enter data map(to: C-block)
do pp = 0, p-1
  if (pp == p_col) then
    A_buf = Ap_row,p_col
  end if
  MPI_Bcast(A_buf,communicator(p_row))
  if (pp == p_row) then
    B_buf = Bp_row,p_col
  end if
  MPI_Bcast(B_buf,communicator(p_col))
  !$omp target map(to: A_buf, B_buf)
  !$omp teams distribute parallel do collapse(2) default(shared) &
  !$omp private(jj,ii,kk)
  do jj = 1, nsize_block
    do ii = 1, nsize_block
      do kk = 1, nsize_block
        C-block(ii,jj) = C-block(ii,jj) + A_buf(ii,kk) * B_buf(kk,jj)
      end do
    end do
  end do
  !$omp end teams distribute parallel do
  !$omp end target
end do
!$omp target exit data map(from: C-block)
deallocate A_buf, B_buf

Note the use of target enter data and target exit data to move C-block to the device before the loop over processors starts, keep it on the device throughout the processor loop, and bring it back to the host after the loop has finished. Without this C-block would be moved to and from the device for every off-load region.

OpenMP off-loading to multiple GPUs

The common approach with MPI+OpenMP off-loading is to run one MPI-process per device. This process then uses OpenMP to off-load compute to the device. What if we want to use OpenMP on a node and create a thread per device and off-load from those threads to the devices? That way we should be able to use a single (MPI) process and use OpenMP to off-load to all the devices in a node. Naively, one might expect this to work fully analogously to MPI+OpenMP off-loading, in that you simply use OpenMP threads instead of MPI processes to off-load from. In practice this is slightly more complicated. The reason is that MPI processes are separate processes that are independent from eachother except for the MPI communication between them. Within OpenMP the threads are created in the context of an overarching process, and the threads share a lot of the process state. As a result if one device changes data in a variable for which a map tofrom clause was specified this triggers a synchronization between the off-loaded kernels. At that moment the benefits of using multiple GPUs is lost.

To avoid the problem described above we need to explicitly suppress synchronization by off-loading the kernels asynchronously by using the nowait clause. The resulting code is (see also openmp-offload-mxm)

      num_devices = 1
#ifdef _OPENMP
      num_devices = omp_get_num_devices()
#endif
      num_per_thread = (nsize-1+num_devices)/num_devices
      !$omp parallel do private(ijlo, ijhi)
      do k_dev = 0, num_devices-1
        ijlo = k_dev*num_per_thread+1
        ijhi = min((k_dev+1)*num_per_thread,nsize)
!
!       Nowait on the next line is essential to suppress synchronization between the threads that drive the GPUs
!$omp   target device(k_dev) map(tofrom: c(:,ijlo:ijhi)) map(to: a, b(:,ijlo:ijhi), ijlo, ijhi, nsize) nowait
!$omp   teams distribute parallel do collapse(2) default(shared) firstprivate(ijlo, ijhi, nsize) private(jj,ii,kk)
        do jj = ijlo, ijhi
          do ii = 1, nsize
            do kk = 1, nsize
              c(ii,jj) = c(ii,jj) + a(ii,kk)* &
                                    b(kk,jj)
            enddo
          enddo
        enddo
!$omp   end teams distribute parallel do
!$omp   end target
      end do
      !$omp end parallel do
!$omp taskwait

Note that we have to use the omp_get_num_devices() OpenMP API call to get the number of devices available in the node. If we are not using an OpenMP capable compiler this function will not be available. The #ifdef statement checks for the presence of the _OPENMP macro that is set by OpenMP capable compilers. As a result this function will only be called when the compiler supports it. To ensure code correctness for compilers that do not support OpenMP we must initialize num_devices explicitly to 1.

In the parallelization of the loop over devices we break matrices c and b into blocks by the second index. That way we off-load contiguous blocks of a data. Nevertheless we are parallelizing only over one index, whereas before in the MPI case we had argued for parallelizing over two indeces. The reason both approaches are reasonable is that with MPI you might create thousands of MPI processes and therefore you want collectively enough iterations to parallelize over that many processes. With OpenMP you parallelize just over the devices within a node. At present the number of devices per node tends to be limited to about 8 (this is anno 2025, if Moore’s law applies to the number of GPUs per node then there might be 256 GPUs per node by 2035, which is not as crazy as it seems given that amplitUDE has 112 cores per node and one could envisage an architecture where you have one GPU per CPU core). For most matrices partitioning them reasonably into 8 pieces is certainly feasible.

Another new clause here is the device clause. Every process has a default device. Therefore, unless we state otherwise all OpenMP threads would connect to the same device. With the device clause we can tell every thread to connect to its own device.

Finally, the nowait clause asynchronously off-loads the compute to the device. Nowait explicitly removes all synchronization so that the device kernels run entirely independently from eachother. As a side effect the code would immediately continue beyond the parallel do construct. Hence we need the taskwait directive to tell OpenMP to wait until all off-loaded kernels have finished.

With this we have a single proces using OpenMP off-loading to multiple GPUs.

MPI + OpenMP off-loading to multiple GPUs

As a final step we might combine MPI with OpenMP off-loading but now using OpenMP to off-load to multiple GPUs. The reason to want to do such a thing is the same that we have seen before in section MPI+OpenMP. When we break a job up into fewer but larger MPI tasks then we need less communication. Larger tasks tend to take longer to complete which tends to make load-balancing worse. We can mitigate the length of larger tasks by parallelizing those tasks over multiple threads. The key difference with respect to OpenMP + Offloading to multiple devices is that in this scenario we want to keep the parts of c on the device between MPI communications and successive off-loads.

The resulting algorithm can be sketched as (see also mpi-openmp-multiple-offload-mxm)

MPI+OpenMP parallel matrix-matrix multiplication

allocate local A-block buffer A_buf
allocate local B-block buffer B_buf
num_per_thread = (nsize_block-1+num_devices)/num_devices
!$omp parallel do private(ijlo, ijhi)
do k_dev = 0, num_devices-1
  ijlo = k_dev*num_per_thread+1
  ijhi = min((k_dev+1)*num_per_thread,nsize_block)
!$omp target enter data device(k_dev) map(to: c(:,ijlo:ijhi)) nowait
enddo
!$omp end parallel do
!$omp taskwait
do pp = 0, p-1
  if (pp == p_col) then
    A_buf = Ap_row,p_col
  end if
  MPI_Bcast(A_buf,communicator(p_row))
  if (pp == p_row) then
    B_buf = Bp_row,p_col
  end if
  MPI_Bcast(B_buf,communicator(p_col))
  !$omp parallel do private(ijlo, ijhi)
  do k_dev = 0, num_devices-1
    ijlo = k_dev*num_per_thread+1
    ijhi = min((k_dev+1)*num_per_thread,nsize_block)
    !$omp target device(k_dev) map(to: A_buf, B_buf(:,ijlo:ijhi))
    !$omp teams distribute parallel do collapse(2) default(shared) &
    !$omp private(jj,ii,kk)
    do jj = ijlo, ijhi
      do ii = 1, nsize_block
        do kk = 1, nsize_block
          C-block(ii,jj) = C-block(ii,jj) + A_buf(ii,kk) * B_buf(kk,jj)
        end do
      end do
    end do
    !$omp end teams distribute parallel do
    !$omp end target
  end do
  !$omp end parallel do
  !$taskwait
end do
!$omp parallel do private(ijlo, ijhi)
do k_dev = 0, num_devices-1
  ijlo = k_dev*num_per_thread+1
  ijhi = min((k_dev+1)*num_per_thread,nsize_block)
!$omp target exit data device(k_dev) map(from: c(:,ijlo:ijhi)) nowait
enddo
!$omp end parallel do
!$omp taskwait
deallocate A_buf, B_buf

With this we have explored every combination of MPI and OpenMP, apart from a mixed model. In a mixed model one could imagine using OpenMP threads to off-load to multiple devices and combine that with further OpenMP threads to do the compute on CPUs. Such a model would need a way to address the difference in the compute power between devices and CPUs. Addressing this issue might prove non-trivial. Furthermore the compute power of modern devices is so much larger than that of CPUs that this model might not make a lot of sense regardless.