|
| 1 | +# Usage |
| 2 | + |
| 3 | +Distributed linear algebra frameworks are the backbone for efficient parallel |
| 4 | +codes in data analytics, scientific computing and machine learning. The central |
| 5 | +idea is that vectors and matrices can be partitioned into potentially |
| 6 | +overlapping chunks which are distributed across a set of workers on which we |
| 7 | +define the usual operations like products and norms. |
| 8 | + |
| 9 | +## Basic example |
| 10 | + |
| 11 | +In this section we take a look on solving the finite difference discretization |
| 12 | +of a Laplace problem in 1D over the domain [0,1]. As a reminder, the Laplace |
| 13 | +problem states to find function u(x) such that Δu(x) = 0 for all x ∈ [0,1]. |
| 14 | +Without boundary conditions the problem is not well-posed, hence we introduce |
| 15 | +the Dirichlet condition u(0) = 1. |
| 16 | + |
| 17 | +Applying the finite difference method with length 0.25 we discretize the problem |
| 18 | +into linear system with 5 unkowns (u₁,...,u₅), which we call degrees of freedom: |
| 19 | +```math |
| 20 | +\frac{1}{4} |
| 21 | +\begin{pmatrix} |
| 22 | +1 & 0 & 0 & 0 & 0 \\ |
| 23 | +0 & -2 & 1 & 0 & 0 \\ |
| 24 | +0 & 1 & -2 & 1 & 0 \\ |
| 25 | +0 & 0 & 1 & -2 & 1 \\ |
| 26 | +0 & 0 & 0 & 1 & -1 |
| 27 | +\end{pmatrix} |
| 28 | +\begin{pmatrix} |
| 29 | +u₁ \\ |
| 30 | +u₂ \\ |
| 31 | +u₃ \\ |
| 32 | +u₄ \\ |
| 33 | +u₅ |
| 34 | +\end{pmatrix} |
| 35 | += |
| 36 | +\begin{pmatrix} |
| 37 | + 1 \\ |
| 38 | +-1 \\ |
| 39 | + 0 \\ |
| 40 | + 0 \\ |
| 41 | + 0 |
| 42 | +\end{pmatrix} |
| 43 | +``` |
| 44 | + |
| 45 | +A detailed derivation can be found in standard numerical analysis lecture notes and books e.g. [these](https://people.sc.fsu.edu/~jburkardt/classes/math2071_2020/poisson_steady_1d/poisson_steady_1d.pdf). The linear system is then solved with |
| 46 | +conjugate gradients. |
| 47 | + |
| 48 | +### Commented Code |
| 49 | + |
| 50 | +To distribute the problem across two workers we have do choose a partitioning. |
| 51 | +Here we arbitrarily assign the first 3 columns and rows to worker 1 and the |
| 52 | +remaining 2 rows and columns to worker 2. |
| 53 | + |
| 54 | +First include the packages which are used. |
| 55 | +```julia |
| 56 | +using PartitionedArrays, SparseArrays, IterativeSolvers |
| 57 | +``` |
| 58 | + |
| 59 | +We want a partitioning into 2 pieces and chose the sequential backend to handle |
| 60 | +the task sequentially so that the code can be executed in a standard Julia REPL (e.g., to simplify debugging). |
| 61 | +```julia |
| 62 | +np = 2 |
| 63 | +backend = SequentialBackend() |
| 64 | +``` |
| 65 | + |
| 66 | +Most of the codes using `PartitionedArrays` start creating a distributed object that for each part contains its part id. We call it `parts`. |
| 67 | +```julia |
| 68 | +parts = get_part_id(backend,np) |
| 69 | +``` |
| 70 | + |
| 71 | +Now, we generate a partitioning of rows and columns. Note that the entry in row 3 |
| 72 | +column 4 is visible to the first worker |
| 73 | +```julia |
| 74 | +neighbors, row_partitioning, col_partitioning = map_parts(parts) do part |
| 75 | + if part == 1 |
| 76 | + ( |
| 77 | + Int32[2], |
| 78 | + IndexSet(part, [1,2,3], Int32[1,1,1]), |
| 79 | + IndexSet(part, [1,2,3,4], Int32[1,1,1,2]) |
| 80 | + ) |
| 81 | + else |
| 82 | + ( |
| 83 | + Int32[1], |
| 84 | + IndexSet(part, [3,4,5], Int32[1,2,2]), |
| 85 | + IndexSet(part, [3,4,5], Int32[1,2,2]) |
| 86 | + ) |
| 87 | + end |
| 88 | +end |
| 89 | +``` |
| 90 | + |
| 91 | +We create information exchangers to manage the synchronization of visible |
| 92 | +shared portions of the sparse matrix and the actual row/col |
| 93 | +```julia |
| 94 | +global_number_of_dofs = 5 |
| 95 | +row_exchanger = Exchanger(row_partitioning,neighbors) |
| 96 | +rows = PRange(global_number_of_dofs,row_partitioning,row_exchanger) |
| 97 | + |
| 98 | +col_exchanger = Exchanger(col_partitioning,neighbors) |
| 99 | +cols = PRange(global_number_of_dofs,col_partitioning,col_exchanger) |
| 100 | +``` |
| 101 | + |
| 102 | +Next we create the sparse matrix entries in COO format in their worker-local |
| 103 | +numbering. A note about the exact values of the sparse matrices can be found |
| 104 | +in the subsection below. |
| 105 | +```julia |
| 106 | +I, J, V = map_parts(parts) do part |
| 107 | + if part == 1 |
| 108 | + ( |
| 109 | + [ 1, 1, 2, 2, 2, 3, 3, 3], |
| 110 | + [ 1, 2, 1, 2, 3, 2, 3, 4], |
| 111 | + 0.25*Float64[1, 0, 0,-2, 1, 1,-1, 0] |
| 112 | + ) |
| 113 | + else |
| 114 | + ( |
| 115 | + [ 1, 1, 2, 2, 2, 3, 3], |
| 116 | + [ 1, 2, 1, 2, 3, 2, 3], |
| 117 | + 0.25*Float64[-1, 1, 1,-2, 1, 1,-1]) |
| 118 | + end |
| 119 | +end |
| 120 | +A = PSparseMatrix(I, J, V, rows, cols, ids=:local) |
| 121 | +``` |
| 122 | + |
| 123 | +Since the previous lines created the local prtions we have to trigger sync |
| 124 | +between the workers. |
| 125 | +```julia |
| 126 | +assemble!(A) |
| 127 | +``` |
| 128 | + |
| 129 | +Construct the right hand side. Note that the first entry of the rhs of worker 2 |
| 130 | +is shared with worker 1. |
| 131 | +```julia |
| 132 | +b = PVector{Float64}(undef, A.rows) |
| 133 | +map_parts(parts,local_view(b, b.rows)) do part, b_local |
| 134 | + if part == 1 |
| 135 | + b_local .= [1.0, -1.0, 0.0] |
| 136 | + else |
| 137 | + b_local .= [0.0, 0.0, 0.0] |
| 138 | + end |
| 139 | +end |
| 140 | +``` |
| 141 | + |
| 142 | +Now the sparse matrix and right hand side of the linear system are assembled |
| 143 | +globally and we can solve problem with cg. With the end in the last line we |
| 144 | +close the parallel environment. |
| 145 | +```julia |
| 146 | +u = IterativeSolvers.cg(A,b) |
| 147 | +``` |
| 148 | + |
| 149 | +### Parallel Code |
| 150 | + |
| 151 | +Now changing the backend to the MPI backend we can solve the problem in parallel. |
| 152 | +This just requires to change the line |
| 153 | +```julia |
| 154 | +backend = SequentialBackend() |
| 155 | +``` |
| 156 | +to |
| 157 | +```julia |
| 158 | +backend = MPIBackend() |
| 159 | +``` |
| 160 | +and including and initializing MPI. Now launching the script with MPI makes the run parallel. |
| 161 | + |
| 162 | +```sh |
| 163 | +$ mpirun -n 2 julia my-script.jl |
| 164 | +``` |
| 165 | + |
| 166 | +Hence the full MPI code is given in the next code box. Note that we have used the `prun` function that automatically includes and initializes MPI for us. |
| 167 | +```julia |
| 168 | +using PartitionedArrays, SparseArrays, IterativeSolvers |
| 169 | + |
| 170 | +np = 2 |
| 171 | +backend = MPIBackend() |
| 172 | + |
| 173 | +prun(backend,np) do parts |
| 174 | + # Construct the partitioning |
| 175 | + neighbors, row_partitioning, col_partitioning = map_parts(parts) do part |
| 176 | + if part == 1 |
| 177 | + ( |
| 178 | + Int32[2], |
| 179 | + IndexSet(part, [1,2,3], Int32[1,1,1]), |
| 180 | + IndexSet(part, [1,2,3,4], Int32[1,1,1,2]) |
| 181 | + ) |
| 182 | + else |
| 183 | + ( |
| 184 | + Int32[1], |
| 185 | + IndexSet(part, [3,4,5], Int32[1,2,2]), |
| 186 | + IndexSet(part, [3,4,5], Int32[1,2,2]) |
| 187 | + ) |
| 188 | + end |
| 189 | + end |
| 190 | + |
| 191 | + global_number_of_dofs = 5 |
| 192 | + |
| 193 | + row_exchanger = Exchanger(row_partitioning,neighbors) |
| 194 | + rows = PRange(global_number_of_dofs,row_partitioning,row_exchanger) |
| 195 | + |
| 196 | + col_exchanger = Exchanger(col_partitioning,neighbors) |
| 197 | + cols = PRange(global_number_of_dofs,col_partitioning,col_exchanger) |
| 198 | + |
| 199 | + # Construct the sparse matrix |
| 200 | + I, J, V = map_parts(parts) do part |
| 201 | + if part == 1 |
| 202 | + ( |
| 203 | + [ 1, 1, 2, 2, 2, 3, 3, 3], |
| 204 | + [ 1, 2, 1, 2, 3, 2, 3, 4], |
| 205 | + 0.25*Float64[1, 0, 0,-2, 1, 1,-1, 0] |
| 206 | + ) |
| 207 | + else |
| 208 | + ( |
| 209 | + [ 1, 1, 2, 2, 2, 3, 3], |
| 210 | + [ 1, 2, 1, 2, 3, 2, 3], |
| 211 | + 0.25*Float64[-1, 1, 1,-2, 1, 1,-1]) |
| 212 | + end |
| 213 | + end |
| 214 | + A = PSparseMatrix(I, J, V, rows, cols, ids=:local) |
| 215 | + assemble!(A) |
| 216 | + |
| 217 | + # Construct the dense right hand side |
| 218 | + b = PVector{Float64}(undef, A.rows) |
| 219 | + map_parts(parts,local_view(b, b.rows)) do part, b_local |
| 220 | + if part == 1 |
| 221 | + b_local .= [1.0, -1.0, 0.0] |
| 222 | + else |
| 223 | + b_local .= [0.0, 0.0, 0.0] |
| 224 | + end |
| 225 | + end |
| 226 | + |
| 227 | + # Solve the linear problem |
| 228 | + u = IterativeSolvers.cg(A,b) |
| 229 | +end |
| 230 | +``` |
| 231 | + |
| 232 | +### Note on Local Matrices |
| 233 | + |
| 234 | +It should be noted that the local matrices are constructed as if they were |
| 235 | +locally assembled on a process without knowledge of the remaining processes. |
| 236 | +Dropping the coefficient 0.25 the global and local matrices look as follows: |
| 237 | + |
| 238 | +``` |
| 239 | + Global Matrix |
| 240 | + P1 P1 P1 P2 P2 |
| 241 | +P1 1 0 0 0 0 |
| 242 | +P1 0 -2 1 0 0 |
| 243 | +P1 0 1 -2 1 0 |
| 244 | +P2 0 0 1 -2 1 |
| 245 | +P2 0 0 0 1 -1 |
| 246 | +
|
| 247 | + = |
| 248 | +
|
| 249 | + Process 1 Portion |
| 250 | + P1 P1 P1 P2 P2 |
| 251 | +P1 1 0 0 0 0 |
| 252 | +P1 0 -2 1 0 0 |
| 253 | +P1 0 1 -1 0 0 |
| 254 | +P2 x x x x x |
| 255 | +P2 x x x x x |
| 256 | +
|
| 257 | + + |
| 258 | +
|
| 259 | + Process 2 Portion |
| 260 | + P1 P1 P1 P2 P2 |
| 261 | +P1 x x x x x |
| 262 | +P1 x x x x x |
| 263 | +P1 0 0 -1 1 0 |
| 264 | +P2 0 0 1 -2 1 |
| 265 | +P2 0 0 0 1 -1 |
| 266 | +``` |
| 267 | + |
| 268 | +## Advanced example |
| 269 | + |
| 270 | +A more complex example can be found in the package [PartitionedPoisson.jl](https://github.com/fverdugo/PartitionedPoisson.jl), |
| 271 | +which describes the assembly of the finite element discretization of a |
| 272 | +Poisson problem in 3D. |
0 commit comments