# What heuristics can be used to minimize the asymptotic matrix bandwidth of a 5-point Laplacian discretization?

I can see that there are multiple heuristics to achieve a matrix with minimum bandwidth. As heuristics, they can't guarantee an optimal solution in polynomial time (after all, the problem is NP-complete).

I wonder, then, for the finite difference discretization of the Poisson problem, how does the bandwidth expand as the problem size increases?

Does the bandwidth approach the number of gridpoints as the spatial resolution decreases? Or does it remain constant? Which approaches produce the smallest bandwidth (asymptotically)? If the cost is too great to produce a minimum bandwidth, then which algorithms balance the computational cost with the width cost?

• Your question title is very general, but the body seems to indicate that you are interested specifically in 5-point discretization of the Laplacian in 2D. Could you change the title or the body so that they agree? – David Ketcheson Jan 19 '12 at 18:02
• How does this title sound? better? – Paul Jan 19 '12 at 18:06
• Yes, now they seem to agree. – David Ketcheson Jan 19 '12 at 18:12

It depends on the shape of your domain and presence of local refinement, but for isotropic domains in $d$ dimensions, the minimum achievable bandwidth scales as $N^{\frac{d-1}{d}}$ where $N$ is the total number of elements (or vertices). You can use a method like Reverse Cuthill-McKee to produce reasonable low-bandwidth orderings.

However, minimizing bandwidth is not the objective of ordering techniques for sparse direct solvers. Indeed, this will lead to $N^{\frac{2d-1}{d}}$ nonzeros in the factors, but optimal orderings (usually found by orderings such as nested dissection which is used by most packages) have $N \log N$ nonzeros in 2D and $N^{4/3}$ nonzeros in 3D. See George, Liu, and Ng's book for details on this result.

## Visual comparison of orderings

The following graphically compares orderings for the 5-point Laplacian on a $12\times 12$ regular grid. For each ordering, I report the number of nonzeros in the $L$ and $U$ factors for a $100\times 100$ grid. Note that the difference is much more pronounced in 3D. You can plot the fill matrices and diagnostics with PETSc:

$cd$PETSC_DIR/src/ksp/ksp/examples/tutorials && make ex2


$./ex2 -m 12 -n 12 -ksp_view -pc_type lu -pc_factor_mat_ordering_type nd -mat_view_draw -draw_pause 2 ### Matrix$A$and factor$L$in natural ordering (1990198 nonzeros) ### Matrix$A$and factor$L$in Reverse Cuthill-McKee ordering (1353100 nonzeros) ### Matrix$A$and factor$L\$ in Nested Dissection ordering (382934 nonzeros)

On the other hand, low-bandwidth orderings are often reasonable orderings (algorithmically) for incomplete factorization and relaxation like SOR and tend to reuse cache well (so are good for throughput). For maximum performance on structured grids, you can skip the column indices (because they have regular structure) and only store the stencil coefficients, or even recompute the coefficients on the fly. It is well worth careful profiling and analysis using a performance model to determine whether optimizations like this are worth the implementation time and reduced flexibility.

A comparison of reordering algorithms as well as their effect on iterative solvers is also here (a few pages down):

http://www.dealii.org/developer/doxygen/deal.II/namespaceDoFRenumbering.html

• This is a nice comparison! But I don't get one thing, why is number of iterations not an integer number? E.g. "we needed an average of 92.2 iterations", so you averaged it over runs, but why then does it change between them? – Alexander Mar 30 '12 at 7:15
• Like it says right above the figures, "The quality of the preconditioner can then be measured by the number of CG iterations required to solve a linear system with [the velocity-velocity block of the Stokes matrix]. For some of the reordering strategies below we record this number for adaptive refinement cycle 3, with 93176 degrees of freedom; because we solve several linear systems with [this matrix block] in the Schur complement [solver for the overall Stokes problem], the average number of iterations is reported." – Wolfgang Bangerth Mar 31 '12 at 8:46