Why parallelization?¶

Top 500 computers in the world: www.top500.org

top500.jpg

Fast computers have several million cores, which need to be used efficiently & simultaneously.

Remember:

Prefix Factor
Kilo $10^3$
Mega $10^6$
Giga $10^9$
Tera $10^{12}$
Peta $10^{15}$
Exa $10^{18}$

What is theoretical peak performance¶

A rough theoretical peak multiplies the number of cores, the clock frequency, and the number of floating-point operations completed per cycle. This is an upper bound for performance, not a prediction of application performance. The actual performance is obtained by running LINPACK benchmark.

In [8]:
# my laptop
cores = 12
clock_GHz = 2.4
flops_per_cycle = 8
number_of_nodes = 10**7

peak_GFLOPS = cores * clock_GHz * flops_per_cycle * number_of_nodes
print(f'Theoretical peak: {peak_GFLOPS:.1f} GFLOP/s')
print(f'                  {peak_GFLOPS/10**6:.4f} PFLOP/s')

Theoretical peak: 2304000000.0 GFLOP/s
                  2304.0000 PFLOP/s

Why parallelization: growth in core counts

Why parallelization: processor performance

Moore's law

Why parallelization: memory latency

OpenMP and multicore execution

OpenMP Mandelbrot example in C++

OpenMP Mandelbrot example in Fortran

Measured Mandelbrot scaling

Python example: serial and parallel Mandelbrot calculations¶

Each point in the complex plane is independent. We therefore parallelize the outer image loop with Numba's prange. The inner loop remains an ordinary serial loop within each worker.

In [9]:
import time
import numpy as np
import matplotlib.pyplot as plt
from numba import njit, prange, get_num_threads, set_num_threads

%matplotlib inline
In [10]:
@njit
def mandelbrot_serial(data, ext, max_steps):
    Nx, Ny = data.shape
    for i in range(Nx):
        x = ext[0] + (ext[1] - ext[0])*i/(Nx - 1.0)
        for j in range(Ny):
            y = ext[2] + (ext[3] - ext[2])*j/(Ny - 1.0)
            c = x + 1j*y
            z = 0j
            iteration = max_steps
            for k in range(max_steps):
                z = z*z + c
                if z.real*z.real + z.imag*z.imag > 4.0:
                    iteration = k
                    break
            data[i, j] = iteration

@njit(parallel=True)
def mandelbrot_parallel(data, ext, max_steps):
    Nx, Ny = data.shape
    for i in prange(Nx):
        x = ext[0] + (ext[1] - ext[0])*i/(Nx - 1.0)
        for j in range(Ny):
            y = ext[2] + (ext[3] - ext[2])*j/(Ny - 1.0)
            c = x + 1j*y
            z = 0j
            iteration = max_steps
            for k in range(max_steps):
                z = z*z + c
                if z.real*z.real + z.imag*z.imag > 4.0:
                    iteration = k
                    break
            data[i, j] = iteration
In [11]:
ext = np.array([-2.0, 1.0, -1.0, 1.0])
Nxy = 1000
max_steps = 1000
serial_data = np.empty((Nxy, Nxy), dtype=np.int32)
parallel_data = np.empty_like(serial_data)

# Warm-up calls compile both functions. Do not include them in the timing.
mandelbrot_serial(serial_data[:10, :10], ext, 10)
mandelbrot_parallel(parallel_data[:10, :10], ext, 10)

mandelbrot_serial(serial_data, ext, max_steps)
mandelbrot_parallel(parallel_data, ext, max_steps)
print('The implementations agree:', np.array_equal(serial_data, parallel_data))

OMP: Info #276: omp_set_nested routine deprecated, please use omp_set_max_active_levels instead.

The implementations agree: True
In [12]:
fig, ax = plt.subplots(figsize=(9, 6))
im = ax.imshow(parallel_data.T, extent=ext, origin='lower', cmap='magma')
ax.set_xlabel(r'$\operatorname{Re}(c)$')
ax.set_ylabel(r'$\operatorname{Im}(c)$')
ax.set_title('Mandelbrot escape iterations')
fig.colorbar(im, ax=ax, label='iteration')
fig.tight_layout()
plt.show()
No description has been provided for this image
In [13]:
get_num_threads()  # gives number of threads on your computer
Out[13]:
12
In [14]:
thread_counts = range(1,get_num_threads()+1)

times = []
for p in thread_counts:
    set_num_threads(p)          # using only limited number of cores after that line
    start = time.perf_counter()
    mandelbrot_parallel(parallel_data, ext, max_steps)
    elapsed = time.perf_counter() - start
    times.append(elapsed)
    print(f'{p:3d} threads: {elapsed:.4f} s')

set_num_threads(get_num_threads())
times = np.array(times)
speedup = times[0]/times

fig, ax = plt.subplots(figsize=(7, 5))
ax.plot(thread_counts, speedup, 'o-', label='measured')
ax.plot(thread_counts, thread_counts, '--', color='0.5', label='ideal')
ax.set_xlabel('number of threads')
ax.set_ylabel('speedup')
ax.grid()
ax.legend()
fig.tight_layout()
plt.show()

  1 threads: 0.6866 s
  2 threads: 0.4794 s
  3 threads: 0.4373 s
  4 threads: 0.4328 s
  5 threads: 0.3276 s
  6 threads: 0.3141 s
  7 threads: 0.2772 s
  8 threads: 0.2285 s
  9 threads: 0.2315 s
 10 threads: 0.2066 s
 11 threads: 0.1901 s
 12 threads: 0.1897 s
No description has been provided for this image

OpenMP reduction example

OpenMP critical section example

In [15]:
%%file pi_examp.cc

#include <iostream>
#include <ctime>
#include <cmath>
#include <omp.h>
using namespace std;

double f(double x){
  return 4.0/(1.0+x*x);
}
double calcPi(int n)
{
  const double dx = 1.0/n;
  double fSum = 0.0;
  #pragma omp parallel for reduction(+:fSum)
  for (int i=0; i<n; ++i){
    double x = (i+0.5)*dx;
    fSum += f(x);
  }
  return fSum*dx;
}

double calcPi_bad(int n)
{
  const double dx = 1.0/n;
  double fSum = 0.0;
  #pragma omp parallel for
  for (int i=0; i<n; ++i){
    double x = (i+0.5)*dx;
    double df = f(x);
    #pragma omp critical
    fSum += df;
  }
  return fSum*dx;
}
int main()
{
  int n=1000000;

  clock_t startTimec = clock();
  double start = omp_get_wtime();
  double fpi;
  for (int j=0; j<100; ++j)
    fpi = calcPi(n);
    //fpi = calcPi_bad(n);
  
  clock_t endTimec = clock();
  double diffc = double(endTimec-startTimec)/CLOCKS_PER_SEC;
  double diff = omp_get_wtime()-start;
  clog<<"clock time : "<<diffc<<"s"<<" with wall time="<<diff<<"s "<<endl;

  std::cout << "At n=" << n << " approximation for pi= " << fpi << " with error= " << std::abs(fpi-M_PI) << std::endl; //<< " " << gpi << std::endl;
}

Overwriting pi_examp.cc
In [17]:
from pathlib import Path

# Homebrew uses versioned names for GNU C++, for example g++-15.
candidates = [p for root in ('/opt/homebrew/bin', '/usr/local/bin')
              for p in Path(root).glob('g++-[0-9]*')]

if not candidates:
    raise FileNotFoundError('GNU C++ was not found. On macOS, install it with: brew install gcc')
    
gxx = str(max(candidates, key=lambda p: int(p.name.rsplit('-', 1)[1])))
print('Using', gxx)

cmd=f'"{gxx}" -fopenmp -O3 -o pi_examp pi_examp.cc'

!{cmd}
!{export OMP_NUM_THREADS=1; pi_examp}
!{export OMP_NUM_THREADS=6; pi_examp}
!{export OMP_NUM_THREADS=12; pi_examp}

Using /opt/homebrew/bin/g++-15
clock time : 0.082733s with wall time=0.082731s 
At n=1000000 approximation for pi= 3.14159 with error= 2.88658e-14
clock time : 0.073281s with wall time=0.014555s 
At n=1000000 approximation for pi= 3.14159 with error= 8.9706e-14
clock time : 0.11407s with wall time=0.0168s 
At n=1000000 approximation for pi= 3.14159 with error= 8.30447e-14

Python example: a parallel reduction¶

Numba recognizes total += value as a reduction inside a prange loop. Each thread accumulates a private partial sum, and Numba combines the partial results after the loop.

In [18]:
@njit(parallel=True)
def midpoint_pi(N):
    dx = 1.0/N
    total = 0.0
    for i in prange(N):
        x = (i + 0.5)*dx
        total += 4.0/(1.0 + x*x)
    return total*dx

midpoint_pi(10)  # compile before examining the timing
for N in (1_000, 100_000, 10_000_000):
    estimate = midpoint_pi(N)
    error = abs(estimate - np.pi)
    print(f'N={N:>10,d}: pi={estimate:.15f}, error={error:.3e}')

N=     1,000: pi=3.141592736923126, error=8.333e-08
N=   100,000: pi=3.141592653598126, error=8.333e-12
N=10,000,000: pi=3.141592653589811, error=1.776e-14
In [19]:
set_num_threads(12)
In [21]:
N = 100_000_000
calls_per_trial = 5

original_threads = get_num_threads()
thread_counts = np.arange(1, original_threads + 1)

thread_counts = range(1,get_num_threads()+1)

times = []
for p in thread_counts:
    set_num_threads(p)          # using only limited number of cores after that line
    start = time.perf_counter()
    for _ in range(calls_per_trial):
        estimate = midpoint_pi(N)
    elapsed = time.perf_counter() - start
    times.append(elapsed)
    print(f'{p:3d} threads: {elapsed:.4f} s')

set_num_threads(get_num_threads())
times = np.array(times)
speedup = times[0]/times

fig, ax = plt.subplots(figsize=(7, 5))
ax.plot(thread_counts, speedup, 'o-', label='measured')
ax.plot(thread_counts, thread_counts, '--', color='0.5', label='ideal')
ax.set_xlabel('number of threads')
ax.set_ylabel('speedup')
ax.grid()
ax.legend()
fig.tight_layout()

  1 threads: 0.3852 s
  2 threads: 0.1788 s
  3 threads: 0.1302 s
  4 threads: 0.0974 s
  5 threads: 0.0790 s
  6 threads: 0.0665 s
  7 threads: 0.0703 s
  8 threads: 0.0634 s
  9 threads: 0.0579 s
 10 threads: 0.0546 s
 11 threads: 0.0493 s
 12 threads: 0.0502 s
No description has been provided for this image

Memory latency and bandwidth

Multicore cache architecture

Modern CPU design and hardware threads

Memory management recommendations

Cache sequence 1

Cache sequence 2

Cache sequence 3

Cache sequence 4

Cache sequence 5

Cache sequence 6

Cache sequence 7

Cache sequence 8

Row-major and column-major memory access

Python example: contiguous and strided access¶

NumPy arrays use C order by default, so the last index is contiguous. These functions perform identical arithmetic but visit memory in different orders.

In [22]:
@njit
def transform_contiguous(a, result):
    for i in range(a.shape[0]):
        for j in range(a.shape[1]):
            result[i, j] = 2.0*a[i, j] + 1.0


@njit
def transform_strided(a, result):
    for j in range(a.shape[1]):
        for i in range(a.shape[0]):
            result[i, j] = 2.0*a[i, j] + 1.0
In [23]:
def typical_time(function, *args, repeat=7, number=5):
    """
    Return the median execution time per function call.
    """
    measurements = []

    for _ in range(repeat):
        start = time.perf_counter()

        for _ in range(number):
            function(*args)

        elapsed = time.perf_counter() - start
        measurements.append(elapsed/number)

    return np.median(measurements)
In [24]:
a = np.random.default_rng(509).random((3000, 3000))
result = np.empty_like(a)

# Compile both functions.
transform_contiguous(a, result)
transform_strided(a, result)

contiguous_time = typical_time(transform_contiguous, a, result)

strided_time = typical_time(transform_strided, a, result)

print(f'contiguous access: {contiguous_time:.6f} s')
print(f'strided access:    {strided_time:.6f} s')
print(
    f'slowdown:          '
    f'{strided_time/contiguous_time:.2f} times'
)

contiguous access: 0.001399 s
strided access:    0.035034 s
slowdown:          25.04 times

MPI and multi-node parallelization