Fast Sparse Matrix Permutation for Mesh-Based Direct Solvers

TL;DR

This work targets the dominant preprocessing bottleneck in sparse Cholesky factorization for mesh-based systems: finding a fill-reducing permutation. It introduces a patch-based, quotient-graph nested-dissection approach that constructs the elimination tree concurrently with the permutation, reuses coarse structure across recursive levels, and relaxes stringent balance and separator-size requirements to speed up computation. The method preserves the ND hierarchy needed for symbolic analysis and integrates with vendor solvers on CPU and GPU, delivering substantial end-to-end speedups (notably up to 6.62x on GPUs and 2.55x on CPUs in various workloads). By trading some permutation quality for speed, the approach significantly reduces preprocessing overhead while maintaining solver robustness, with strong practical impact for triangle-mesh pipelines in graphics applications. The results demonstrate the method’s scalability, and the authors outline extensions to volumetric meshes and GPU-based separator computation as future work.

Abstract

We present a fast sparse matrix permutation algorithm tailored to linear systems arising from triangle meshes. Our approach produces nested-dissection-style permutations while significantly reducing permutation runtime overhead. Rather than enforcing strict balance and separator optimality, the algorithm deliberately relaxes these design decisions to favor fast partitioning and efficient elimination-tree construction. Our method decomposes permutation into patch-level local orderings and a compact quotient-graph ordering of separators, preserving the essential structure required by sparse Cholesky factorization while avoiding its most expensive components. We integrate our algorithm into vendor-maintained sparse Cholesky solvers on both CPUs and GPUs. Across a range of graphics applications, including single factorizations, repeated factorizations, our method reduces permutation time and improves the sparse Cholesky solve performance by up to 6.27x.

Figures (21)

• Figure 1: Scalable fill-reducing permutation computation on a triangular mesh. Our patch-based, nested-dissection-style method compresses the triangular mesh domain into a small quotient graph and reuses this compression across recursive calls, improving scalability over state-of-the-art methods. We achieve a 6.62$\times$ speedup over the highly optimized NVIDIA cuDSS by accelerating its permutation computation.
• Figure 2: Fraction of time spent on matrix permutation to end-to-end time for solving mean curvature flow Desbrun:1999:IFO for different inputs.
• Figure 3: Nested dissection example. Illustration of a single nested-dissection step. From left to right: the input matrix $A$, its graph $G$, the corresponding elimination tree $T$, the permuted graph $G'$, and the sparsity pattern of the Cholesky factor $L$. The separator vertices $S$ are ordered after the vertices in $g_1$ and $g_2$. The elimination tree encodes the dependencies induced by this order: computations for $g_1$ and $g_2$ are independent, while computations for $S$ can proceed only after both have completed.
• Figure 4: Elimination tree construction. We compute the root separator $\{7,2,5\}$. For the resulting subgraphs $\textit{g}_\ell=\{1,4,6\}$ and $\textit{g}_r=\{3,8,9\}$, we recurse to compute the next-level separators (e.g., $\{4\}$ and $\{8\}$), continuing until reaching the desired depth. We store the full binary tree in a 1D array.
• Figure 5: Separator computation. For the current subgraph $\textit{g}$, we form a quotient graph $\textit{q}$, bipartition $\textit{q}$, lift the partition back to $V_\textit{g}$ to obtain a separator superset, refine it to reduce its size while maintaining balance, and then extract the left and right subgraphs for the next recursion level.