How should I build a 2D 5-point stencil Laplacian matrix in parallel?

Question

I'm making a simple eigenvalue solver with SLEPc, using a 5-point stencil and the finite difference method. I want to be able to assemble the matrix in parallel.

My first thought was just to use MatGetOwnershipRange() to get range_begin and range_end, then use tests inside the loop to determine if that row represents a point on the edge. All those conditionals make the loop extremely slow, and elsewhere in the PETSc documentation it recommends building separate loops for the corners and edges.

How can this be done in a parallel program?

Is there a better way that works for parallel code not mentioned in the laplacian example code?

Answer 1

Those conditionals are the least of your worries. If matrix insertion/assembly is slow, it almost always means that the matrix was not correctly preallocated. After that, check that the parallel distribution is good (e.g. use a 2D decomposition instead of 1D) and that your assembly loop mostly sets values in rows owned by the given process. The 2D decomposition is somewhat tricky to program yourself, so you might want to look at an example like this which uses a DMDA to manage the parallel decomposition. Finally, prefer using one call to MatSetValues() per row or per element, instead of inserting scalar entries.

Answer 2

Consider defining a boundary and setting its values to zero. In 2d:

0 0 0 0 0
0 . . . 0
0 . . . 0
0 . . . 0
0 0 0 0 0

This way the inner loop can be free of conditionals, and you can carve up your matrix as desired for parallelisation.