# CS 5220

## Applications of Parallel Computers

### Stationary iterations

Prof David Bindel

Please click the play button below.

### Fixed Point Iteration

$$x_{k+1} = f(x_k) \rightarrow x_* = f(x_*)$$

### Iterative Idea

• $$f$$ is a contraction if $$\|f(x)-f(y)\| < \|x-y\|$$, $$x \neq y$$.
• $$f$$ has a unique fixed point $$x_* = f(x_*)$$.
• For $$x_{k+1} = f(x_k)$$, $$x_k \rightarrow x_*$$.
• If $$\|f(x)-f(y)\| < \alpha \|x-y\|$$, $$\alpha < 1$$, for all $$x, y$$, then $\|x_k-x_*\| < \alpha^k \|x-x_*\|$
• Looks good if $$\alpha$$ not too near 1...

### Stationary Iterations

Write $$Ax = b$$ as $$A = M-K$$; get fixed point of $M x_{k+1} = K x_k + b$ or $x_{k+1} = (M^{-1} K) x_k + M^{-1} b.$

• Convergence if $$\rho(M^{-1} K) < 1$$
• Best case for convergence: $$M = A$$
• Cheapest case: $$M = I$$

### Stationary Iterations

• Realistic: choose something between

 Jacobi $$M = \operatorname{diag}(A)$$ Gauss-Seidel $$M = \operatorname{tril}(A)$$

### Reminder: Discretized 2D Poisson Problem

$$(Lu)_{i,j} = h^{-2} \left( 4u_{i,j}-u_{i-1,j}-u_{i+1,j}-u_{i,j-1}-u_{i,j+1} \right)$$

### Jacobi on 2D Poisson

Assuming homogeneous Dirichlet boundary conditions

for step = 1:nsteps

for i = 2:n-1
for j = 2:n-1
u_next(i,j) = ...
( u(i,j+1) + u(i,j-1) + ...
u(i-1,j) + u(i+1,j) )/4 - ...
h^2*f(i,j)/4;
end
end
u = u_next;

end

Basically do some averaging at each step.

### Parallel version (5 point stencil)

 Boundary values: white Data on P0: green Ghost cell data: blue

### Parallel version (9 point stencil)

 Boundary values: white Data on P0: green Ghost cell data: blue

### Parallel version (5 point stencil)

Communicate ghost cells before each step.

### Parallel version (9 point stencil)

Communicate in two phases (EW, NS) to get corners.

### Parallel version (9 point stencil)

Communicate in two phases (EW, NS) to get corners.

### Gauss-Seidel on 2D Poisson

for step = 1:nsteps

for i = 2:n-1
for j = 2:n-1
u(i,j) = ...
( u(i,j+1) + u(i,j-1) + ...
u(i-1,j) + u(i+1,j) )/4 - ...
h^2*f(i,j)/4;
end
end

end

Bottom values depend on top; how to parallelize?

### Red-Black Gauss-Seidel

Red depends only on black, and vice-versa.
Generalization: multi-color orderings

### Red black Gauss-Seidel step

  for i = 2:n-1
for j = 2:n-1
if mod(i+j,2) == 0
u(i,j) = ...
end
end
end

for i = 2:n-1
for j = 2:n-1
if mod(i+j,2) == 1
u(i,j) = ...
end
end
end

### Parallel red-black Gauss-Seidel sketch

At each step

• Send black ghost cells
• Update red cells
• Send red ghost cells
• Update black ghost cells

### More Sophistication

• Successive over-relaxation (SOR): extrapolate G-S
• Block Jacobi: $$M$$ a block diagonal matrix from $$A$$
• Other block variants similar
• Alternating Direction Implicit (ADI): alternately solve on vertical lines and horizontal lines
• Multigrid

Mostly the opening act for Krylov methods.