Introduction

We have been talking about multithreaded algorithms.
These were algorithms which run on a Dynamic Multithreading concurrency platform, a concurrency platform with a load balancing scheduler and supported nested parallelism (spawning) and parallel loops (parallel-for).
We discussed the performance measures for our algorithms of work, `T_1`, total amount of computation needed to be done or the amount of time on one processor; and span, `T_(infty)`, the amount of time if we had as many processors as desired.
We gave the work and span laws as lower bounds on how long a parallel algorithm takes to run on p processors, `T_P`.
We also proved Brent's Theorem which gave an upper bound.
We analyzed a parallel version of our Fibonacci algorithm using span-and-sync, and computed that its parallelism, `T_1/T_(infty)`, was exponential.
Thus far, we haven't given an example with pseudo-code showing the use of parallel for, so we begin today with an example.

Parallel Loops

Many algorithms contain loops all of whose iterations can operate in parallel.
Although we could parallelize these using the spawn and sync keywords, it is more convenient to instead introduce a new keyword parallel which we'll put before for statements whose iterations we want to run in parallel.
For example let `A = (a_(ij))` be an `n times n` matrix and let `vec x= (x_j)` be a column vector.

The following is a parallel algorithm for computing `A vec x`, the multiplication of `A` with `x`.

MAT-VEC(A, x)
1 n = A.rows
2 let y be a new vector of length n
3 parallel for i = 1 to n:
4     y[i] = 0
5 parallel for i = 1 to n:
6     for j = 1 to n:
7          y[i] = y[i] + a[i][j] * x[j]
8 return y

Simulating Parallel For using Spawn and Sync

A compiler can implement each parallel for loop as a divide-and-conquer subroutine using nested parallelism.

For example, lines 5-7 can be implement with the call MAT-VEC-MAIN-LOOP(A, x, y, n, 1, n) where MAT-VEC-MAIN-LOOP is defined as:

MAT-VEC-MAIN-LOOP (A, x, y, n, i, i')
1 if i == i':
2     for j = 1 to n
3         y[i] = y[i] + a[i][j] * x[j]
4 else: 
      mid = floor((i+ i')/2)
5     spawn MAT-VEC-MAIN-LOOP(A, x, y, n, i, mid)
6     MAT-VEC-MAIN-LOOP(A, x, y, n, mid + 1, i')
7     sync

Analyzing MAT-VEC(A, x)

To calculate `T_1(n)` for MAT-VEC we just remove the parallel keywords and analyze the resulting code.
Because of the doubly nested loops this gives, `T_1(n) = Theta(n^2)`.
As we are converting parallel-for loops to spawn sync calls in a balanced fashion, the span to execute these calls will be the depth of the tree of recursive calls and so logarithmic in the number of tree leaves, in this case `n`.
More precisely, for a parallel loop with `n` iterations in which the `i`th iteration has span `iter_(infty)(i)`, the span is
`T_(infty)(n) = Theta(log n) + max_(1 le i le n) iter_(infty)(i)`.
In the case of MAT-VEC, the first parallel for loop has span `Theta(log n)` as each iteration is `Theta(1)`. In the case of the second parallel for loop, we have an inner serial for loop so get the span to be `Theta(n)`.
Hence, the span of the whole algorithm will be `Theta(n)`, as it will be dominated by the second parallel for loop.
The parallelism of this whole algorithm is thus `(Theta(n^2))/(Theta(n)) = Theta(n)`.

In-Class Exercise

Rewrite MAT-VEC(A,x) using spawn and sync so that the inner for loop on j is done is a parallel fashion.
What is `T_1/T_infty` in this case?
Would the added complexity be worth it if we had a fixed number of processors? Why or why not?
Post your solution to the Feb 16 In-Class Exercise Thread.

Race Conditions

A multithreaded algorithm is deterministic if it always does the same thing on the same input, no matter how the instructions are scheduled on the multicore computer; otherwise, it is said to be nondeterministic.
A determinancy race occurs when two logically parallel instructions access the same memory location and at least one of the instructions performs a write.
```
RACE-EXAMPLE
1 x = 0
2 parallel for i = 1 to 2
3     x = x + 1
4 print x
```
After line 1 two parallel strands try to increment x in line 3, and print x in line 4.
This might print 1 or it might print 2.
If both strands read x before either increments it, then 1 will be output.
If one strand increments x before the other reads it then the output will be 2.

Preventing Race Conditions

There are many different ways used to try to prevent race conditions: exclusion locks, synchronization blocks, etc.
To keep things simple, we shall insist that strands that operate in parallel are independent: they have no determinancy races among them.
So in a parallel for construct, all the iterations should be independent.
Between a spawn and the corresponding sync, the code of the spawned child should be independent of the code of the parent, including code executed by additional spawned or called children.

Another Race Example

Here is another example of code which has a race condition:

MAT-VEC-WRONG(A, x)
1 n = A.rows
2 let y be a new vector of length n
3 parallel for i = 1 to n
4    y[i] = 0
5 parallel for i = 1 to n
6    parallel for j = 1 to n
7        y[i] = y[i] + a[i][j] * x[j]
8 return y

This code is like our code from earlier, except the second inner for loop has been parallelized, giving it a span of `Theta(log n)` rather than `Theta(n)`.
Unfortunately, the updates in line 7 are not all independent and so we have a determinancy race.

Chess Program Socrates

The chess program Socrates was programmed using the Cilk programming language and competed in matches in late 1990s, early 2000s.
It was developed on a 32 processor system but ran in matches on a 512 processor system.
During development the programmers came up with an optimization that reduced `T_(32) = 65` to `T_(32)' = 40` seconds on one of the testing benchmarks.
The original code was eventually favored over the optimization because it was determined the original would be faster on 512 processors.
Here is the analysis:
- The original code had `T_1 =2048 sec` and `T_(infty) = 1 sec`.
- Suppose Brent's Theorem was exactly true rather than an inequality. Then `T_P = T_1/P + T_(infty)` and the run time on 32 processors would be
  `T_(32) = 2048/32 + 1 = 65`.
- The optimization had `T_1' =1024 sec` and `T_(infty)' = 8 sec`. Using the same approximation we get
  `T_(32)' = 1024/32 + 8 = 40`.
- For 512 processors, though, `T_(512) = 2048/512 + 1 = 5` and `T_(512)' = 1024/512 + 8 = 10`
From this, we see that work and span together provide a better means of extrapolating performance than just measuring running times.

Multithreaded Matrix Multiplication

If we parallelize the standard matrix multiply algorithm, we get the following:

P-SQUARE-MATRIX-MULTIPLY(A, B)
1 n = A.rows
2 let C be a new n x n matrix
3 parallel for i = 1 to n
4     parallel for j = 1 to n
5         c[i][j] = 0
6         for k = 1 to n
7             c[i][j] = c[i][j] + a[i][k] * b[k][j]
8 return C

The serialization of this has three nested for loops up to n and so `T_1(n) = Theta(n^3)`.
The span is `T_(infty)(n) = Theta(n)` because we follow a `Theta(log n)` path through the tree of recursion for the outer parallel for, followed by another path of `Theta(log n)` path through the tree of recursion for the inner parallel for and then we serial execute the innermost for `n` times.
So we have `Theta(log n) + Theta(log n) + Theta(n) = Theta(n)`.
This gives a parallelism of `(Theta(n^3))/(Theta(n)) = Theta(n^2)`.

Multithreaded Matrix Multiplication - Take 2

There are several ways we can improve the parallel runtime of matrix multiplication.
One could start with a faster single processor algorithm such as Strassen's algorithm -- the book shows the analysis of such an algorithm.
Or we can try to improve the use of parallelism in the existing algorithm
We take this approach and give a first pass at reducing `T_(infty)(n)` to `Theta(log^2 n)`. This can actually be improved to `Theta(log n)`.
To start, let's assume `n` is a power of 2. We can split the matrices `A` and `B` into four quadrants `A_(11)`, `A_(12)`, `A_(21)`, `A_(22)` and `B_(11)`, `B_(12)`, `B_(21)`, `B_(22)` each with `n/2 times n/2` entries. We can then multiply these quadrants to get `C` using matrix multiply for `2 times 2` matrices:
`A = [[A_(11),A_(12)],[A_(21),A_(22)]]`, `B = [[B_(11),B_(12)],[B_(21),B_(22)]]`, `C = [[C_(11),C_(12)],[C_(21),C_(22)]]`
where `C = A cdot B` is given by:
`[[C_(11),C_(12)],[C_(21),C_(22)]] = [[A_(11),A_(12)],[A_(21),A_(22)]][[B_(11),B_(12)],[B_(21),B_(22)]]`
`= [[A_(11)B_(11) + A_(12)B_(21),A_(11)B_(12) + A_(12)B_(22)],[A_(21)B_(11) + A_(22)B_(21),A_(21)B_(12) + A_(22)B_(22)]]`
We use the above to give a parallel divide-and-conquer algorithm for matrix multiply.

Race Conditions, Chess, Matrix Multiplication

Outline