Multithreading

Parallel Algorithms

David Mayerich

Scalable Tissue Imaging and Modeling (STIM) Laboratory

Department of Electrical and Computer Engineering

Cullen College of Engineering

University of Houston

David Mayerich

STIM Laboratory, University of Houston

Memory Allocation

  • stack vs. heap

  • A stack is a data structure that is very fast at removing stuff added last

  • Basic allocations happen on the stack:

David Mayerich

STIM Laboratory, University of Houston

int a = 12;
int b = 300;
float x = 3.14;
float y[3];

float* z = (float*) malloc(sizeof(float) * 3);

previous allocations

a = 12
b = 300
x = 3.14
y[0]
y[1]
y[2]
z*
z[0]
z[1]
z[2]

heap

Practical Stack vs. Heap Allocation

  • The stack is fast - resets with scope:

David Mayerich

STIM Laboratory, University of Houston

int N = 100;
float z = 0.0;
for(int i = 0; i < N; i++) {
	float x = sin( a[i] );
    float y = cos( b[i] );
    z = x * y;
}

  • The stack has a limited size (stack overflow)

  • The heap is slower to allocate/free

previous allocations

N = 100
z = 0.0
i = 0
x
y

Dynamic Allocation

  • DO NOT use variables in static allocations:



     

    • Not supported in recent versions of the C standard

    • Can result in a stack overflow

David Mayerich

STIM Laboratory, University of Houston

  • DO use dynamic allocation:

int n;
scanf("size of array: %d", n);
float A[n];
float* A = (float*) malloc( n * sizeof(float));

OR (using C++):

float* A = new float[n];

Segmentation Faults

David Mayerich

STIM Laboratory, University of Houston

int N = 2000;
//.....
double* x;
x = (double*) malloc(N * sizeof(double));
x[2000] = 3.14159;		// possible segmentation fault

Heap Memory

Case 1 - allocation is contiguous in memory

Case 2 - allocation is not contiguous in memory

A[6] overwrites B[0]
(no error reported)

segmentation fault

A[0]

B[0]

A[0]

B[0]

Multidimensional Arrays

  • Stack allocation for 2D arrays:
     

David Mayerich

STIM Laboratory, University of Houston

  • Heap allocation for indirect addressing:

double A[20][20];
A[7][3] = 7.0 * 3.0;
int M = 20;		int N = 20;
//.....
double** A;
A = (double**) malloc(M * sizeof(double*));
for(int i = 0; i < M; i++) {
	A[i] = (double*) malloc(N * sizeof(double));
}
A[7][3] = 7.0 * 3.0;
A[N][0]
A[N][1]
A[N][N-1]
A[1][0]
A[1][1]
A[1][N-1]
double*
double**

A

double
\vdots

Multidimensional Arrays

  • Indirect addressing is slower

    • requires following a pointer in the heap

    • difficult to cache


       

  • Heap allocation for direct addressing:

David Mayerich

STIM Laboratory, University of Houston

int M = 20;		int N = 20;
//.....
double* A;
A = (double*) malloc(M * N * sizeof(double*));
A[7 * N + 3] = 7.0 * 3.0;

integer math is fast

A

A[0]

A[1*M+0]

A[1*M+2]

A[(N-1)*M+0]

Storing Matrices

  • Matrices are just 2D arrays, but they're very common and have their own terminology:

    • row major order format stores elements such that row values are contiguous

    • column major order format stores elements such that column values are contiguous

David Mayerich

STIM Laboratory, University of Houston

\begin{bmatrix} m_{1,1} & m_{1,2} & m_{1,3}\\ m_{2,1} & m_{2,2} & m_{2,3}\\ m_{3,1} & m_{3,2} & m_{3,3}\\ \end{bmatrix}

Row-Major Order

Column-Major Order

m_{1,1}
m_{2,1}
m_{3,1}
m_{1,2}
m_{2,2}
m_{3,2}
m_{1,3}
m_{2,3}
m_{3,3}
m_{1,1}
m_{2,1}
m_{3,1}
m_{1,2}
m_{2,2}
m_{3,2}
m_{1,3}
m_{2,3}
m_{3,3}
\begin{bmatrix} m_{1,1} & m_{1,2} & m_{1,3}\\ m_{2,1} & m_{2,2} & m_{2,3}\\ m_{3,1} & m_{3,2} & m_{3,3}\\ \end{bmatrix}
\begin{bmatrix} m_{1,1} & m_{1,2} & m_{1,3}\\ m_{2,1} & m_{2,2} & m_{2,3}\\ m_{3,1} & m_{3,2} & m_{3,3}\\ \end{bmatrix}

Higher Dimensions

  • Three dimensions \((x, y, z)\):

David Mayerich

STIM Laboratory, University of Houston

int X = 20;		int Y = 30;		int Z = 10;

float* A = (float*) malloc(X * Y * Z * sizeof(float));

// access (7, 8, 9):
A[9 * X * Y + 8 * X + 7] = 7.0f * 3.0f;

Function Pointers

  • Functions, like data, have locations in memory

  • If we define two functions with similar signature:

David Mayerich

STIM Laboratory, University of Houston

double cosine_expansion(double x) {
  const double x2 = x * x;
  double term = 1.0;
  double sum = 1.0;
  for (int k = 1; k < max_terms; ++k) {
    term *= -x2 / ((2.0 * k - 1.0) * (2.0 * k));
    if (sum + term == sum) break;
    sum += term;
  }
  return sum;
}
double sine_expansion(double x) {
  const double x2 = x * x;
  double term = x;
  double sum = x;
  for (int k = 1; k < max_terms; ++k) {
    term *= -x2 / ((2.0 * k) * (2.0 * k + 1.0));
    if (sum + term == sum) break;
    sum += term;
  }
  return sum;
}

  • These functions can be selectively called:

double (*trig)(double);			// create a function pointer "trig"
if(user_input == "cos") 		// if the user wants to call cosine
	trig = cosine_expansion;	// assign that function to trig 
else trig = sine_expansion;		// otherwise assign the sine function
double y = trig(x);				// call the user-specified trig function

Callback Functions

  • These are useful for creating callback functions: user-specified functions called by non-user-specified functions

David Mayerich

STIM Laboratory, University of Houston

// user-defined function that prints the
// mouse coordinates
void user_mouse_click(int px, int py){
	printf("Mouse coordinate: %d, %d", px, py);
}

Callback Functions

  • These are useful for creating callback functions: user-specified functions called by non-user-specified functions

David Mayerich

STIM Laboratory, University of Houston

// user-defined function that prints the
// mouse coordinates
void user_mouse_click(int px, int py){
	printf("Mouse coordinate: %d, %d", px, py);
}

// tell the GUI (someone else's code)
// what function to call when the mouse
// is clicked
GUI.set_mouse_func(user_mouse_click);

Threads and Concurrency

  • A thread is a lightweight process that manages an independent set of instructions

  • Each thread has its own call stack, including registers

  • The main thread is the execution pipeline for main() and used to launch and control other concurrent threads:

David Mayerich

STIM Laboratory, University of Houston

execution time

shared state

heap

libraries
global functions

main()

registers

call stack:
inst 2
launch t0
inst 3
launch t1
inst 4
...
 

thread0()

registers

call stack:

inst a1
inst a2
...

thread1()

...

C++ Thread Library

  • C++11 standard supports threads

  • Launch independent threads by assigning a function pointer to the constructor:

David Mayerich

STIM Laboratory, University of Houston

#include <thread>

void func(float p1, int p2, ...) {
	// do independent work using
    // p1, p2, etc.
}
// ...
// create a thread by passing:
//	1) the function to execute
//	2) all parameters required by the function
std::thread t(func, p1, p2, ...);

Matrix-Vector Multiplication (Example)

  • Consider the following operation:

David Mayerich

STIM Laboratory, University of Houston

\mathbf{y}=\mathbf{Mx} = \begin{bmatrix} m_{1,1} & \cdots & m_{1,N} \\ \vdots & \ddots & \vdots \\ m_{N,1} & \cdots & m_{N,N} \end{bmatrix} \begin{bmatrix} x_1 \\ \vdots \\ x_N \end{bmatrix}

  • This can be computed as a set of \(N\) dot (inner) vector products:

\mathbf{y}= \mathbf{Mx} = \begin{bmatrix} y_1 \\ \vdots \\ y_N \end{bmatrix} \quad \text{where} \quad y_i = \begin{bmatrix} m_{i,1} & \cdots & m_{i,N} \end{bmatrix} \begin{bmatrix} x_{1} \\ \cdots \\ x_N \end{bmatrix}

Inner Product Implementation

David Mayerich

STIM Laboratory, University of Houston

// Calculate the dot (inner) product of two vectors v1 and v2
// of length N and store the result in c
void dot(float* v1, float* v2, float* c, float N) {

}

Inner Product Implementation

David Mayerich

STIM Laboratory, University of Houston

// Calculate the dot (inner) product of two vectors v1 and v2
// of length N and store the result in c
void dot(float* v1, float* v2, float* c, float N) {
	// create a local register to accumulate the result
	float r = 0;
}

Inner Product Implementation

David Mayerich

STIM Laboratory, University of Houston

// Calculate the dot (inner) product of two vectors v1 and v2
// of length N and store the result in c
void dot(float* v1, float* v2, float* c, float N) {
	// create a local register to accumulate the result
	float r = 0;
    // for each element in both vectors
    for(int i = 0; i < N; i++)
}

Inner Product Implementation

David Mayerich

STIM Laboratory, University of Houston

// Calculate the dot (inner) product of two vectors v1 and v2
// of length N and store the result in c
void dot(float* v1, float* v2, float* c, float N) {
	// create a local register to accumulate the result
	float r = 0;
    // for each element in both vectors
    for(int i = 0; i < N; i++)
    	r += v1[i] * v2[i];
}

Inner Product Implementation

David Mayerich

STIM Laboratory, University of Houston

// Calculate the dot (inner) product of two vectors v1 and v2
// of length N and store the result in c
void dot(float* v1, float* v2, float* c, float N) {
	// create a local register to accumulate the result
	float r = 0;
    // for each element in both vectors
    for(int i = 0; i < N; i++)
    	r += v1[i] * v2[i];
    *c = r;
}

Implement Single-Process Multiplication

  • Define the dot product function (previous slide):

David Mayerich

STIM Laboratory, University of Houston

void dot(float* v1, float* v2, float* c, float N);

  • Single-process implementation of \(\mathbf{y}=\mathbf{Mx}\):

// allocate memory
float* x = new float[N];
float* y = new float[N];
float* M = new float[N*N];
// fill values x, M ...

Implement Single-Process Multiplication

  • Define the dot product function (previous slide):

David Mayerich

STIM Laboratory, University of Houston

void dot(float* v1, float* v2, float* c, float N);

  • Single-process implementation of \(\mathbf{y}=\mathbf{Mx}\):

// allocate memory
float* x = new float[N];
float* y = new float[N];
float* M = new float[N*N];
// fill values x, M ...

// iterate across all matrix rows
for(int ri = 0; ri < N; ri++)
	dot(x, &M[ri * N], &y[ri], N);

Implement Multithreading

  • Define the dot product function (previous slide):

David Mayerich

STIM Laboratory, University of Houston

void dot(float* v1, float* v2, float* c, float N);

  • Single-process implementation of \(\mathbf{y}=\mathbf{Mx}\):

// allocate memory
float* x = new float[N];
float* y = new float[N];
float* M = new float[N*N];
// fill values x, M ...
// create an array of threads
std::thread* t = new std::thread[N];
//launch one thread for each row
for(int ri = 0; ri < N; ri++)
	t[ri] = std::thread(dot, y, &M[ri * N], &x[ri], N);

Synchronization

  • This loop will return immediately

  • HOWEVER, y will be invalid until all threads complete

  • Implement a synchronization barrier with join():

David Mayerich

STIM Laboratory, University of Houston

for(int ri = 0; ri < N; ri++)
	t[ri] = std::thread(dot, y, &M[ri * N], &x[ri], N);
for(int ri = 0; ri < N; ri++)
	t[ri] = std::thread(dot, y, &M[ri * N], &x[ri], N);

// synchronization barrier -----------
for(int ri = 0; ri < N; ri++)
	// join() returns when the thread is finished
	t[ri].join();

Summary

  • By default, a C/C++ program is operating along a single main thread

  • Child threads can be spawned using an std::thread object

    • Create the std::thread

    • Pass the work to the thread's constructor as a function pointer (*func)

  • Child threads have their own state space: parameters passed to func through the constructor are unique to that thread

  • Child threads can access global memory via pointers

  • When a thread is launched, the main thread returns immediately: you can't expect the thread to be finished without checking!

  • Threads can be synchronized with main using the join() function

David Mayerich

STIM Laboratory, University of Houston