Lecture 5: Cilk Examples

Shankar Balachandran
CSE Department, IIT Madras
shankar@[Link]

The lecture is partly based on

Charles Leiserson’s Slides on Cilk
 The Master Method
The Master Method for solving recurrences
applies to recurrences of the form
T(n) = a T(n/b) + f (n) ,*
where a ¸ 1, b > 1, and f is asymptotically
positive.

IDEA: Compare nlogba with f (n) .

*The unstated base case is T(n) = (1) for

sufficiently small n.
 Master Method — CASE 1
T(n) = a T(n/b) + f (n)

nlogba À f (n)

Specifically, f (n) = O(nlogba – e) for some

constant e > 0.
Solution: T(n) = (nlogba) .
 Master Method — CASE 2
T(n) = a T(n/b) + f (n)

nlogba ¼ f (n)

Specifically, f (n) = (nlogba lgkn) for some

constant k ¸ 0.
Solution: T(n) = (nlogba lgk+1n) .
 Master Method — CASE 3
T(n) = a T(n/b) + f (n)

nlogba ¿ f (n)

Specifically, f (n) = (nlogba + e) for some

constant e > 0 and f (n) satisfies the
regularity condition that a f (n/b) · c f (n)
for some constant c < 1.
Solution: T(n) = (f (n)) .
 Master Method Summary
T(n) = a T(n/b) + f (n)

CASE 1: f (n) = O(nlogba – e), constant e > 0

 T(n) = (nlogba) .
CASE 2: f (n) = (nlogba lgkn), constant k  0
 T(n) = (nlogba lgk+1n) .
CASE 3: f (n) = (nlogba + e ), constant e > 0,
and regularity condition
 T(n) = ( f (n)) .
 Master Method Quiz
• T(n) = 4 T(n/2) + n
nlogba = n2 À n ) CASE 1: T(n) = (n2).
• T(n) = 4 T(n/2) + n2
nlogba = n2 = n2 lg0n ) CASE 2: T(n) = (n2lg n).
• T(n) = 4 T(n/2) + n3
nlogba = n2 ¿ n3 ) CASE 3: T(n) = (n3).
• T(n) = 4 T(n/2) + n2/ lg n
Master method does not apply!
 Outline
• Recurrences (Review)
• Matrix Multiplication
• Merge Sort
• Tableau Construction
• Conclusion
Square-Matrix Multiplication
c11 c12 L c1n a11 a12 L a1n b11 b12 L b1n
c21 c22 L c2n a21 a22 L a2n b21 b22 L b2n
M M O M
= M M O M
£ M M O M
cn1 cn2 L cnn an1 an2 L ann bn1 bn2 L bnn

C A B
n
cij = a ik bkj
k=1

Assume for simplicity that n = 2k.

Recursive Matrix Multiplication
Divide and conquer —
C11 C12 A11 A12 B11 B12
= £
C21 C22 A21 A22 B21 B22

A11B11 A11B12 A12B21 A12B22

= +
A21B11 A21B12 A22B21 A22B22

8 multiplications of (n/2) £ (n/2) matrices.

1 addition of n £ n matrices.
Matrix Multiply in Pseudo-Cilk
cilk void Mult(*C, *A, *B, n) {
float *T = Cilk_alloca(n*n*sizeof(float));
h base case & partition matrices i
spawn Mult(C11,A11,B11,n/2);
spawn Mult(C12,A11,B12,n/2);
spawn Mult(C22,A21,B12,n/2);
spawn Mult(C21,A21,B11,n/2);
spawn Mult(T11,A12,B21,n/2);
spawn Mult(T12,A12,B22,n/2);
spawn Mult(T22,A22,B22,n/2);
spawn Mult(T21,A22,B21,n/2);
sync;
spawn Add(C,T,n);
sync;
return;
}
Absence of type
C = A¢ B declarations.
Matrix Multiply in Pseudo-Cilk
cilk void Mult(*C, *A, *B, n) {
float *T = Cilk_alloca(n*n*sizeof(float));
h base case & partition matrices i
spawn Mult(C11,A11,B11,n/2);
spawn Mult(C12,A11,B12,n/2);
spawn Mult(C22,A21,B12,n/2);
spawn Mult(C21,A21,B11,n/2);
spawn Mult(T11,A12,B21,n/2);
spawn Mult(T12,A12,B22,n/2);
spawn Mult(T22,A22,B22,n/2);
spawn Mult(T21,A22,B21,n/2);
sync;
spawn Add(C,T,n);
sync;
return;
}
Coarsen base cases
C = A¢ B for efficiency.
Matrix Multiply in Pseudo-Cilk
cilk void Mult(*C, *A, *B, n) {
float *T = Cilk_alloca(n*n*sizeof(float));
h base case & partition matrices i
spawn Mult(C11,A11,B11,n/2);
spawn Mult(C12,A11,B12,n/2);
Also need a row-
spawn Mult(C22,A21,B12,n/2);
size argument for
spawn Mult(C21,A21,B11,n/2);
spawn Mult(T11,A12,B21,n/2);
array indexing.
spawn Mult(T12,A12,B22,n/2);
spawn Mult(T22,A22,B22,n/2);
spawn Mult(T21,A22,B21,n/2);
sync;
spawn Add(C,T,n);
sync; Submatrices are
return;
} produced by pointer
C = A¢ B calculation, not
copying of elements.
 Matrix Multiply in Pseudo-Cilk
cilk void Mult(*C, *A, *B, n) {
float *T = Cilk_alloca(n*n*sizeof(float));
h base case & partition matrices i
spawn Mult(C11,A11,B11,n/2);
spawn Mult(C12,A11,B12,n/2);
spawn Mult(C22,A21,B12,n/2);
spawn Mult(C21,A21,B11,n/2);
spawn Mult(T11,A12,B21,n/2);
spawn Mult(T12,A12,B22,n/2);
spawn Mult(T22,A22,B22,n/2);
spawn Mult(T21,A22,B21,n/2);
sync;
spawn Add(C,T,n); cilk void Add(*C, *T, n) {
sync; h base case & partition matrices i
return; spawn Add(C11,T11,n/2);
} spawn Add(C12,T12,n/2);
spawn Add(C21,T21,n/2);
spawn Add(C22,T22,n/2);
C = A¢ B sync;
return;
C=C+T }
Work of Matrix Addition
cilk void Add(*C, *T, n) {
h base case & partition matrices i
spawn Add(C11,T11,n/2);
spawn Add(C12,T12,n/2);
spawn Add(C21,T21,n/2);
spawn Add(C22,T22,n/2);
sync;
return;
}

Work: A1(n) = 4?A1(n/2) + (1)

= (n2) — CASE 1
nlogba = nlog24 = n2 À (1) .
 Span of Matrix Addition
cilk void Add(*C, *T, n) {
h base case & partition matrices i
spawn Add(C11,T11,n/2);
spawn Add(C12,T12,n/2);
maximum spawn Add(C21,T21,n/2);
spawn Add(C22,T22,n/2);
sync;
return;
}

Span: A1(n) = A?1(n/2) + (1)

= (lg n) — CASE 2
nlogba = nlog21 = 1 ) f (n) = (nlogba lg0n) .
Work of Matrix Multiplication
cilk void Mult(*C, *A, *B, n) {
float *T = Cilk_alloca(n*n*sizeof(float));
h base case & partition matrices i
spawn Mult(C11,A11,B11,n/2);
spawn Mult(C12,A11,B12,n/2);
8 M
spawn Mult(T21,A22,B21,n/2);
sync;
spawn Add(C,T,n);
sync;
return;
}

Work: M1(n) = 8?M1(n/2) +A1(n) + (1)

= 8 M1(n/2) + (n2)
= (n3) — CASE 1
nlogba = nlog28 = n3 À (n2) .
Span of Matrix Multiplication
cilk void Mult(*C, *A, *B, n) {
float *T = Cilk_alloca(n*n*sizeof(float));
h base case & partition matrices i
spawn Mult(C11,A11,B11,n/2);
spawn Mult(C12,A11,B12,n/2);
8 M
spawn Mult(T21,A22,B21,n/2);
sync;
spawn Add(C,T,n);
sync;
return;
}

Span: M1(n) = M?1(n/2) + A1(n) + (1)

= M1(n/2) + (lg n)
= (lg2 n) — CASE 2
nlogba = nlog21 = 1 ) f (n) = (nlogba lg1n) .
Parallelism of Matrix Multiply
Work: M1(n) = (n3)
Span: M1(n) = (lg2n)

M1(n)
Parallelism: = (n3/lg2n)
M1(n)

For 1000 £ 1000 matrices,

parallelism ¼ (103)3/102 = 107.
 Stack Temporaries
cilk void Mult(*C, *A, *B, n) {
float *T = Cilk_alloca(n*n*sizeof(float));
h base case & partition matrices i
spawn Mult(C11,A11,B11,n/2);
spawn Mult(C12,A11,B12,n/2);
M
spawn Mult(T21,A22,B21,n/2);
sync;
spawn Add(C,T,n);
sync;
return;
}

In hierarchical-memory machines (especially chip

multiprocessors), memory accesses are so expensive that
minimizing storage often yields higher performance.

IDEA: Trade off parallelism for less storage.

No-Temp Matrix Multiplication
cilk void MultA(*C, *A, *B, n) {
// C = C + A * B
h base case & partition matrices i
spawn MultA(C11,A11,B11,n/2);
spawn MultA(C12,A11,B12,n/2);
spawn MultA(C22,A21,B12,n/2);
spawn MultA(C21,A21,B11,n/2);
sync;
spawn MultA(C21,A22,B21,n/2);
spawn MultA(C22,A22,B22,n/2);
spawn MultA(C12,A12,B22,n/2);
spawn MultA(C11,A12,B21,n/2);
sync;
return;
}

Saves space, but at what expense?

Work of No-Temp Multiply
cilk void MultA(*C, *A, *B, n) {
// C = C + A * B
h base case & partition matrices i
spawn MultA(C11,A11,B11,n/2);
spawn MultA(C12,A11,B12,n/2);
spawn MultA(C22,A21,B12,n/2);
spawn MultA(C21,A21,B11,n/2);
sync;
spawn MultA(C21,A22,B21,n/2);
spawn MultA(C22,A22,B22,n/2);
spawn MultA(C12,A12,B22,n/2);
spawn MultA(C11,A12,B21,n/2);
sync;
return;
}

Work: M1(n) = 8?M1(n/2) + (1)

= (n3) — CASE 1
 Span of No-Temp Multiply
cilk void MultA(*C, *A, *B, n) {
// C = C + A * B
h base case & partition matrices i
spawn MultA(C11,A11,B11,n/2);
spawn MultA(C12,A11,B12,n/2);
maximum spawn MultA(C22,A21,B12,n/2);
spawn MultA(C21,A21,B11,n/2);
sync;
spawn MultA(C21,A22,B21,n/2);
spawn MultA(C22,A22,B22,n/2);
maximum spawn MultA(C12,A12,B22,n/2);
spawn MultA(C11,A12,B21,n/2);
sync;
return;
}

Span: M1(n) = 2?M1(n/2) + (1)

= (n) — CASE 1
Parallelism of No-Temp Multiply
Work: M1(n) = (n3)
Span: M1(n) = (n)

M1(n)
Parallelism: = (n2)
M1(n)
For 1000 £ 1000 matrices,
parallelism ¼ (103)3/103 = 106.
Faster in practice!
 Testing Synchronization
Cilk language feature: A programmer can
check whether a Cilk procedure is ―synched‖
(without actually performing a sync) by
testing the pseudovariable SYNCHED:
•SYNCHED = 0 ) some spawned children
might not have returned.
•SYNCHED = 1 ) all spawned children
have definitely returned.
 Best of Both Worlds
cilk void Mult1(*C, *A, *B, n) {// multiply & store
h base case & partition matrices i
spawn Mult1(C11,A11,B11,n/2); // multiply & store
spawn Mult1(C12,A11,B12,n/2);
spawn Mult1(C22,A21,B12,n/2);
spawn Mult1(C21,A21,B11,n/2);
if (SYNCHED) {
spawn MultA1(C11,A12,B21,n/2); // multiply & add
spawn MultA1(C12,A12,B22,n/2);
spawn MultA1(C22,A22,B22,n/2);
spawn MultA1(C21,A22,B21,n/2);
} else {
float *T = Cilk_alloca(n*n*sizeof(float));
spawn Mult1(T11,A12,B21,n/2); // multiply & store
This code is just as parallel
spawn Mult1(T12,A12,B22,n/2);
spawn Mult1(T22,A22,B22,n/2);
as the original, but it only
spawn Mult1(T21,A22,B21,n/2);
sync;
}
uses more space if runtime
spawn Add(C,T,n); // C = C + T
sync;
return;
parallelism actually exists.
}
Ordinary Matrix Multiplication
n
cij = a ik bkj
k=1

IDEA: Spawn n2 inner Work: (n3)

products in parallel.
Compute each inner Span: (lg n)
product in parallel. Parallelism: (n3/lg n)

BUT, this algorithm exhibits

poor locality and does not
exploit the cache hierarchy
of modern microprocessors,
especially CMP’s.
 Outline
• Recurrences (Review)
• Matrix Multiplication
• Merge Sort
• Tableau Construction
• Conclusion
 Merging Two Sorted Arrays
void Merge(int *C, int *A, int *B, int na, int nb) {
while (na>0 && nb>0) {
if (*A <= *B) {
*C++ = *A++; na--;
} else {
*C++ = *B++; nb--;
}
} Time to merge n
while (na>0) {
*C++ = *A++; na--; elements = (n).
?
}
while (nb>0) {
*C++ = *B++; nb--;
}
}

3 12 19 46

4 14 21 23
 Merge Sort
cilk void MergeSort(int *B, int *A, int n) {
if (n==1) {
B[0] = A[0];
} else {
int *C;
C = (int*) Cilk_alloca(n*sizeof(int));
spawn MergeSort(C, A, n/2);
spawn MergeSort(C+n/2, A+n/2, n-n/2);
sync;
Merge(B, C, C+n/2, n/2, n-n/2);
}
}
3 4 12 14 19 21 33 46
merge
3 12 19 46 4 14 21 33
merge
3 19 12 46 4 33 14 21
merge
19 3 12 46 33 4 21 14
 Work of Merge Sort
cilk void MergeSort(int *B, int *A, int n) {
if (n==1) {
B[0] = A[0];
} else {
int *C;
C = (int*) Cilk_alloca(n*sizeof(int));
spawn MergeSort(C, A, n/2);
spawn MergeSort(C+n/2, A+n/2, n-n/2);
sync;
Merge(B, C, C+n/2, n/2, n-n/2);
}
}

Work: T1(n) = 2 ?T1(n/2) + (n)

= (n lg n) — CASE 2
nlogba = nlog22 = n ) f (n) = (nlogba lg0n) .
 Span of Merge Sort
cilk void MergeSort(int *B, int *A, int n) {
if (n==1) {
B[0] = A[0];
} else {
int *C;
C = (int*) Cilk_alloca(n*sizeof(int));
spawn MergeSort(C, A, n/2);
spawn MergeSort(C+n/2, A+n/2, n-n/2);
sync;
Merge(B, C, C+n/2, n/2, n-n/2);
}
}

? (n/2) + (n)
Span: T1(n) = T1
= (n) — CASE 3
nlogba = nlog21 = 1 ¿ (n) .
Parallelism of Merge Sort
Work: T1(n) = (n lg n)
Span: T1(n) = (n)

T1(n)
Parallelism: = (lg n)
T1(n)

We need to parallelize the merge!

 Parallel Merge
0 na/2 na
A · A[na/2] ¸ A[na/2]

Recursive Binary search Recursive

merge merge

B · A[na/2] ¸ A[na/2] na ¸ nb
0 j j+1 nb

KEY IDEA: If the total number of elements to be

merged in the two arrays is n = na + nb, the total
number of elements in the larger of the two
recursive merges is at most (3/4)
? n.
 Parallel Merge
cilk void P_Merge(int *C, int *A, int *B,
int na, int nb) {
if (na < nb) {
spawn P_Merge(C, B, A, nb, na);
} else if (na==1) {
if (nb == 0) {
C[0] = A[0];
} else {
C[0] = (A[0]<B[0]) ? A[0] : B[0]; /* minimum */
C[1] = (A[0]<B[0]) ? B[0] : A[0]; /* maximum */
}
} else {
int ma = na/2;
int mb = BinarySearch(A[ma], B, nb);
spawn P_Merge(C, A, B, ma, mb);
spawn P_Merge(C+ma+mb, A+ma, B+mb, na-ma, nb-mb);
sync;
}
}

Coarsen base cases for efficiency.

 Span of P_Merge
cilk void P_Merge(int *C, int *A, int *B,
int na, int nb) {
if (na < nb) {
M
} else {
int ma = na/2;
int mb = BinarySearch(A[ma], B, nb);
spawn P_Merge(C, A, B, ma, mb);
spawn P_Merge(C+ma+mb, A+ma, B+mb, na-ma, nb-mb);
sync;
}
}

? (3n/4) + (lg n)
Span: T1(n) = T1
= (lg2n) — CASE 2
nlogba = nlog4/31 = 1 ) f (n) = (nlogba lg1n) .
 Work of P_Merge
cilk void P_Merge(int *C, int *A, int *B,
int na, int nb) {
if (na < nb) {
M
} else {
int ma = na/2;
int mb = BinarySearch(A[ma], B, nb);
spawn P_Merge(C, A, B, ma, mb);
spawn P_Merge(C+ma+mb, A+ma, B+mb, na-ma, nb-mb);
sync;
}
}

Work: T1(n) = T?1(n) + T1((1–)n) + (lg n),

where 1/4 ·  · 3/4 .
CLAIM: T1(n) = (n) .
Parallelism of P_Merge

Work: T1(n) = (n)

Span: T1(n) = (lg2n)

T1(n)
Parallelism: = (n/lg2n)
T1(n)
 Parallel Merge Sort
cilk void P_MergeSort(int *B, int *A, int n) {
if (n==1) {
B[0] = A[0];
} else {
int *C;
C = (int*) Cilk_alloca(n*sizeof(int));
spawn P_MergeSort(C, A, n/2);
spawn P_MergeSort(C+n/2, A+n/2, n-n/2);
sync;
spawn P_Merge(B, C, C+n/2, n/2, n-n/2);
}
}
Work of Parallel Merge Sort
cilk void P_MergeSort(int *B, int *A, int n) {
if (n==1) {
B[0] = A[0];
} else {
int *C;
C = (int*) Cilk_alloca(n*sizeof(int));
spawn P_MergeSort(C, A, n/2);
spawn P_MergeSort(C+n/2, A+n/2, n-n/2);
sync;
spawn P_Merge(B, C, C+n/2, n/2, n-n/2);
}
}

Work: T1(n) = 2 T1(n/2) + (n)

= (n lg n) — CASE 2
Span of Parallel Merge Sort
cilk void P_MergeSort(int *B, int *A, int n) {
if (n==1) {
B[0] = A[0];
} else {
int *C;
C = (int*) Cilk_alloca(n*sizeof(int));
spawn P_MergeSort(C, A, n/2);
spawn P_MergeSort(C+n/2, A+n/2, n-n/2);
sync;
spawn P_Merge(B, C, C+n/2, n/2, n-n/2);
}
}

? (n/2) + (lg2n)
Span: T1(n) = T1
= (lg3n) — CASE 2
nlogba = nlog21 = 1 ) f (n) = (nlogba lg2n) .
Parallelism of Merge Sort

Work: T1(n) = (n lg n)

Span: T1(n) = (lg3n)

T1(n)
Parallelism: = (n/lg2n)
T1(n)
 Outline
• Recurrences (Review)
• Matrix Multiplication
• Merge Sort
• Tableau Construction
• Conclusion
 Tableau Construction
Problem: Fill in an n £ n tableau A, where
A[i, j] = f ( A[i, j–1], A[i–1, j], A[i–1, j–1] ).
00 01 02 03 04 05 06 07 Dynamic
10 11 12 13 14 15 16 17 programming
20 21 22 23 24 25 26 27 • Longest common
subsequence
30 31 32 33 34 35 36 37
• Edit distance
40 41 42 43 44 45 46 47
• Time warping
50 51 52 53 54 55 56 57
60 61 62 63 64 65 66 67
Work: (n2).
70 71 72 73 74 75 76 77
 Recursive Construction
n
Cilk code
I II spawn I;
sync;
n spawn II;
spawn III;
sync;
III IV spawn IV;
sync;
 Recursive Construction
n
Cilk code
I II spawn I;
sync;
n spawn II;
spawn III;
sync;
III IV spawn IV;
sync;

? 1(n/2) + (1)
Work: T1(n) = 4T
= (n2) — CASE 1
 Recursive Construction
n
Cilk code
I II spawn I;
sync;
n spawn II;
spawn III;
sync;
III IV spawn IV;
sync;

? 1(n/2) + (1)
Span: T1(n) = 3T
= (nlg3) — CASE 1
Analysis of Tableau Construction
Work: T1(n) = (n2)
Span: T1(n) = (nlg3)
¼ (n1.58)

T1(n)
Parallelism: ¼ (n0.42)
T1(n)
A More-Parallel Construction
n spawn I;
sync;
spawn II;
I II IV spawn III;
sync;
spawn IV;
spawn V;
n III V VII spawn VI
sync;
spawn VII;
spawn VIII;
VI VIII IX sync;
spawn IX;
sync;
A More-Parallel Construction
n spawn I;
sync;
spawn II;
I II IV spawn III;
sync;
spawn IV;
spawn V;
n III V VII spawn VI
sync;
spawn VII;
spawn VIII;
VI VIII IX sync;
spawn IX;
sync;

? 1(n/3) + (1)
Work: T1(n) = 9T
= (n2) — CASE 1
A More-Parallel Construction
n spawn I;
sync;
spawn II;
I II IV spawn III;
sync;
spawn IV;
spawn V;
n III V VII spawn VI
sync;
spawn VII;
spawn VIII;
VI VIII IX sync;
spawn IX;
sync;

? 1(n/3) + (1)
Span: T1(n) = 5T
= (nlog35) — CASE 1
Analysis of Revised Construction
Work: T1(n) = (n2)
Span: T1(n) = (nlog35)
¼ (n1.46)

T1(n)
Parallelism: ¼ (n0.54)
T1(n)
More parallel by a factor of
(n0.54)/(n0.42) = (n0.12) .
 Puzzle
What is the largest parallelism that
can be obtained for the tableau-
construction problem using Cilk?

• You may only use basic Cilk control

constructs (spawn, sync) for
synchronization.
• No locks, synchronizing through
memory, etc.