CSE 463: Parallel and Distributed Systems Fall 2024 Final Exam Solutions

CSE 463: Parallel and Distributed // Step 1: Forward wave propagation

Systems for it = 2 to nt -1:
for ix = 2 to nx -1:
Fall 2024 Final Examination  Model Solutions for iy = 2 to ny -1:
Alexandria University  Computers and Systems for iz = 2 to nz -1:
s[iz ,iy ,ix , it ] = WaveEq (s ,v ,
Engineering it ,ix ,iy , iz )

// Step 2: Backward wave propagation

Question 1: Multithreading and Message for it = nt -1 downto 1:
Passing [42] for ix = 2 to nx -1:
for iy = 2 to ny -1:
1.a) Find Largest and Smallest in List [E][15] for iz = 2 to nz -1:
g[iz ,iy ,ix , it ] = WaveEq (g ,v ,
Algorithm: Divide the array among p threads. Each it ,ix ,iy , iz )
thread nds local min/max, then combine results.
// Step 3: Image generation
// Shared variables for ix = 0 to nx -1:
shared : A[n], global_max , global_min for iy = 0 to ny -1:
shared : lock for iz = 0 to nz -1:
I[iz ,iy , ix ] = 0
// Thread function for thread_id = 0 to p -1 for it = 0 to nt -1:
function FindMinMax ( thread_id ): I[iz ,iy , ix ] += s[iz ,iy , ix , it
start = thread_id * (n/p ) ]
end = ( thread_id + 1) * (n/p ) - 1 * g[iz ,iy ,ix , it
]
local_max = A[ start ] return I
local_min = A[ start ]

for i = start +1 to end : Time Complexity:

if A[i ] > local_max :
local_max = A[i ] ˆ Wave propagation: O(nt · nx · ny · nz) each
if A[i ] < local_min :
local_min = A[i ] ˆ Image generation: O(nx · ny · nz · nt)
ˆ Total: O(nt · nx · ny · nz) = O(N · T ) where N =
// Critical section nx · ny · nz
acquire ( lock )
if local_max > global_max :
global_max = local_max Space Complexity:
if local_min < global_min :
global_min = local_min ˆ Waveelds s, g : O(nx · ny · nz · nt) each
release ( lock ) ˆ Image I : O(nx · ny · nz)
// Main ˆ Total: O(N · T )
global_max = - infinity
global_min = + infinity 1.c.2) Multithreaded Parallel Algorithm [D][3]
parallel_for tid = 0 to p -1: function ParallelSeismicImaging () :
FindMinMax ( tid ) // Parallelize spatial loops ( ix ,iy , iz )
// Time loop must remain sequential ( data
dependency )
1.b) Time Complexity & Cost Optimality [I][10]
Time Complexity Analysis: // Step 1: Forward propagation
for it = 2 to nt -1:
ˆ Each thread processes n/p elements: O(n/p) # pragma omp parallel for collapse (3)
for ix = 2 to nx -1:
ˆ Combining results: O(p) in worst case (sequential lock) for iy = 2 to ny -1:
ˆ Total parallel time: Tp = O(n/p + p) for iz = 2 to nz -1:
s[iz ,iy ,ix , it ] = WaveEq (s ,v
Serial Time: T1 = O(n) ,...)
barrier () // Sync before next timestep
Speedup: S = T1
Tp = n
n/p+p
// Step 2: Backward propagation
Cost: C = p · Tp = p · (n/p + p) = n + p2 for it = nt -1 downto 1:
# pragma omp parallel for collapse (3)
Cost Optimality Condition: For cost-optimal
√ algo- for ix , iy , iz in spatial domain :
rithm, C = O(n), so we need p2 ≤ n, thus p = O( n). g[iz ,iy ,ix , it ] = WaveEq (g ,v ,...)
barrier ()
Making it Cost-Optimal:
// Step 3: Image - fully parallel
√ # pragma omp parallel for collapse (3)
1. Use p = O( n) processors
for ix = 0 to nx -1:
2. Use tree-based reduction for combining: reduces O(p) for iy = 0 to ny -1:
to O(log p) for iz = 0 to nz -1:
3. With reduction: Tp = O(n/p + log p), Cost = O(n + sum = 0
for it = 0 to nt -1:
p log p) sum += s [...] * g [...]
4. With p = O(n/ log n): Cost = O(n)  cost-optimal I[iz ,iy , ix ] = sum

1.c.1) Serial Seismic Algorithm [I][10]

// Serial algorithm for seismic imaging
1.c.3) Complexity Analysis [D][2]
function SeismicImaging () : Time Complexity (with p threads):

1
CSE 463: Parallel and Distributed Systems Fall 2024 Final Exam Solutions

ˆ Wave propagation: O(nt · N

p)  time loop sequential Solution: Use shared memory to cache tiles of B .
ˆ Image generation: O( N p·nt ) __global__ kernel MatVecOptimized (A , B , C , m , n)
  :
ˆ Total: Tp = O Np·T __shared__ float Bs [ TILE_SIZE ]

Space Complexity: O(N ·T ) (unchanged  shared mem- i = blockIdx .x * blockDim .x + threadIdx .x

sum = 0
ory)
// Process B in tiles
Speedup: S = T1
Tp = N ·T
N ·T /p =p for tile = 0 to (n / TILE_SIZE - 1) :
// Collaboratively load B tile
Speedup approaches p for large spatial domains. Limited if threadIdx .x < TILE_SIZE :
by time-step barriers. Bs [ threadIdx .x] = B[ tile * TILE_SIZE
+ threadIdx .x]
1.c.4) Hybrid MPI + Multithreading [I][2] __syncthreads ()
// Domain decomposition across nodes ( MPI )
// Thread parallelism within each node // Compute partial sum using shared mem
if i < m:
function HybridSeismic () : for j = 0 to TILE_SIZE -1:
MPI_Init () col = tile * TILE_SIZE + j
rank = MPI_Comm_rank () sum += A[i ][ col ] * Bs [j]
__syncthreads ()
// Decompose spatial domain across nodes
my_x_start , my_x_end = decompose (nx , rank ) if i < m:
C[i ] = sum
for it = 2 to nt -1:
// Exchange ghost cells with neighbors
MPI_SendRecv ( boundaries to neighbors ) Computation-to-Memory Ratio Analysis:
// Local computation with threads ˆ Original: Each thread loads n elements of B . Total
# pragma omp parallel for collapse (2)
for ix = my_x_start to my_x_end : loads for B : m · n
for iy , iz in local domain : ˆ Optimized: B loaded once per block. Loads for B :
s [...] = WaveEq (...) n · ⌈m/blockSize⌉
MPI_Barrier () ˆ Improvement factor: ≈ blockSize (e.g., 256×)

  Time Complexity:
Time Complexity: T = O N ·T
p·M + O(T · comm)
ˆ Serial: T1 = O(m · n)
where M = number of MPI processes, p = threads per
ˆ Parallel (with p threads): Tp = O(m · n/p)
process
Space: O(N · T /M ) per node (distributed) Speedup: S = p (ideal speedup)

Question 2: Data Parallelism [58] Space: O(m · n) for matrix + O(n + m) for vectors
2.a) Matrix-Vector Multiplication [E][20] 2.c.1) Out-of-Core Matrix-Vector Multiply [I][5]
Matrix Am×n multiplied by vector Bn gives result Cm : When matrix doesn't t in memory, process in blocks:
// Serial out -of - core matrix - vector multiply
n−1
X // Matrix A is 2^ n x 2^n , stored on disk
Ci = Ai,j · Bj
j=0 function OutOfCoreMatVec ( A_file , B , C , N):
// N = 2^ n
block_size = available_memory / N
// CUDA / Data - parallel pseudo - code
// Each thread computes one element of C
// Initialize result
for i = 0 to N -1:
__global__ kernel MatVec (A , B , C , m , n ):
C[i ] = 0
i = blockIdx . x * blockDim .x + threadIdx .x
// Process matrix in row blocks
if i < m:
for block = 0 to ( N/ block_size - 1) :
sum = 0
row_start = block * block_size
for j = 0 to n -1:
row_end = ( block + 1) * block_size
sum += A[i ][ j] * B[ j]
C[i ] = sum
// Load block of rows from disk
A_block = load_from_disk ( A_file ,
// Launch with m threads
row_start , row_end )
MatVec < < <( m +255) /256 , 256 > > >(A , B , C , m , n)
// Compute partial results
Each thread: Reads row i of A (n elements) + vector B for i = 0 to block_size - 1:
for j = 0 to N - 1:
(n elements), performs n multiply-adds. C[ row_start + i] +=
2.b) Improved Algorithm with Shared Memory A_block [ i ][ j] * B [j]
[I][20] free ( A_block )
Problem: In 2.a, every thread reads entire vector B 
redundant global memory accesses. return C

2
CSE 463: Parallel and Distributed Systems Fall 2024 Final Exam Solutions

Execution Time: SpMV_CSR < < <( m +255) /256 , 256 > > >(...)

ˆ I/O: O(N 2 /bandwidth)  read entire matrix 2.e) Load Balancing for SpMV [H][3]
ˆ Computation: O(N 2 ) multiply-adds 1. Load Balance Problem:
ˆ Total: T = O(N 2 /B + N 2 ) = O(N 2 ) dominated by
I/O In CSR with one thread per row:

2.c.2) Multi-GPU Extension [D][3] ˆ Rows have dierent numbers of non-zeros

// Multi - GPU matrix - vector multiplication ˆ Thread with dense row takes much longer
function MultiGPU_MatVec ( A_file , B , C , N , ˆ Sparse rows nish quickly and wait (idle)
num_gpus ):
rows_per_gpu = N / num_gpus ˆ SIMT execution causes warp divergence

// Broadcast B to all GPUs 2. Load-Balanced Algorithm:

for gpu = 0 to num_gpus - 1:
// Approach : One thread per non - zero element
cudaSetDevice ( gpu )
// Use segmented reduction for row sums
cudaMemcpy ( d_B [ gpu ], B , N* sizeof ( float ))
__global__ kernel SpMV_Balanced ( values , col_idx ,
// Each GPU processes its row block
row_idx , B , C , nnz , m):
parallel for gpu = 0 to num_gpus - 1:
// row_idx [ k] = which row element k belongs
cudaSetDevice ( gpu )
to
row_start = gpu * rows_per_gpu
tid = blockIdx . x * blockDim .x + threadIdx .x
// Stream rows from disk to GPU
if tid < nnz :
for block in gpu 's row range :
// Each thread handles one non - zero
load A_block from disk
row = row_idx [ tid ]
cudaMemcpyAsync ( d_A , A_block , ...)
col = col_idx [ tid ]
MatVecKernel < < <... > > >( d_A , d_B [ gpu ] ,
product = values [ tid ] * B[ col ]
d_C [ gpu ], ...)
// Atomic add to result ( or use
// Copy results back
reduction )
cudaMemcpy ( C + row_start , d_C [ gpu ], ...)
atomicAdd (& C[ row ] , product )
// Results already in correct positions
// Alternative : Use warp - level primitives
return C
// Assign work units of fixed size , not rows

Speedup: ≈ num_gpus (limited by disk I/O bandwidth) Better Approach  CSR-Vector:

2.d.1) Sparse Matrix Data Structure [I][5]
ˆ Assign one warp (32 threads) per row
Compressed Sparse Row (CSR) Format:
ˆ Threads in warp cooperate on row using reduction
ˆ values[]: Non-zero values (size = nnz) ˆ More uniform work distribution
ˆ col_idx[]: Column index for each value (size = nnz) ˆ Handles variable row lengths within warp
ˆ row_ptr[]: Start index in values for each row (size = // CSR - Vector : One warp per row
m+1) __global__ kernel SpMV_Vector ( values , col_idx ,
row_ptr , B , C , m):
Example: int warp_id = ( blockIdx .x * blockDim .x
+ threadIdx .x) / 32
 
1 0 2
int lane = threadIdx .x % 32
A = 0 3 0
4 0 5 if warp_id < m:
float sum = 0
int start = row_ptr [ warp_id ]
values = [1, 2, 3, 4, 5] int end = row_ptr [ warp_id + 1]
col_idx = [0, 2, 1, 0, 2]
row_ptr = [0, 2, 3, 5] // Each lane processes subset of row
for ( int j = start + lane ; j < end ; j +=
Space: O(nnz + m) instead of O(m · n) 32)
sum += values [j ] * B[ col_idx [j ]]
2.d.2) Sparse Matrix-Vector (SpMV) [D][2]
// CSR SpMV - one thread per row // Warp - level reduction
__global__ kernel SpMV_CSR ( values , col_idx , for ( int offset = 16; offset > 0; offset
row_ptr , B , C , m): /= 2)
row = blockIdx . x * blockDim .x + threadIdx .x sum += __shfl_down_sync (0 xffffffff ,
sum , offset )
if row < m :
sum = 0 if ( lane == 0)
start = row_ptr [ row ] C[ warp_id ] = sum
end = row_ptr [ row + 1]

for idx = start to end - 1:

col = col_idx [ idx ] End of Solutions
sum += values [ idx ] * B[ col ]

C[ row ] = sum

// Launch