PaperlayerView exact revision on arXiv

Source-linked AI summary

An efficient multi-core implementation of a novel HSS-structured multifrontal solver using randomized sampling

Pieter Ghysels, Xiaoye S. Li, Francois-Henry Rouet, Samuel Williams, Artem Napov

arXiv:1502.07405v1cs.MS

TL;DR

Dense frontal operations can bottleneck sparse-system solvers, despite low-rank structure in many applications. The paper uses randomized HSS compression and fast ULV factorization in a multifrontal solver; it reports speedups up to 7× over the pure multifrontal algorithm.

  • Problem

    Dense operations on multifrontal frontal matrices can become a bottleneck, although many applications have low-rank structure in those matrices.

    Source: §1 · p. 1

  • Method

    The solver approximates larger frontal matrices with HSS using randomized sampling and interpolative decompositions, then applies ULV factorization.

    Source: §1 · pp. 1–2

  • Results

    Speedups up to 7× over the pure multifrontal algorithm, with favorable comparison to commercial PARDISO.

    Source: §7 · p. 22

  • Takeaways & Limitations

    The solver has lower computational complexity than the pure multifrontal method, but the constants involved are much larger.

    Source: §7 · p. 22

    For many consecutive, highly accurate solves, pure multifrontal factorization is probably optimal; the HSS solver is faster in the solve phase when only a few digits of accuracy are required.

    Source: §6.1 · p. 16

Abstract

from arXiv · show

We present a sparse linear system solver that is based on a multifrontal variant of Gaussian elimination, and exploits low-rank approximation of the resulting dense frontal matrices. We use hierarchically semiseparable (HSS) matrices, which have low-rank off-diagonal blocks, to approximate the frontal matrices. For HSS matrix construction, a randomized sampling algorithm is used together with interpolative decompositions. The combination of the randomized compression with a fast ULV HSS factorization leads to a solver with lower computational complexity than the standard multifrontal method for many applications, resulting in speedups up to 7 fold for problems in our test suite. The implementation targets many-core systems by using task parallelism with dynamic runtime scheduling. Numerical experiments show performance improvements over state-of-the-art sparse direct solvers. The implementation achieves high performance and good scalability on a range of modern shared memory parallel systems, including the Intel Xeon Phi (MIC). The code is part of a software package called STRUMPACK -- STRUctured Matrices PACKage, which also has a distributed memory component for dense rank-structured matrices.

1 Introduction

The paper develops a fully structured, parallel multifrontal solver for general sparse linear systems, using randomized HSS compression to address the cost of dense frontal operations. It combines this approach with an adaptive rank construction and task parallelism for modern multicore systems.

1 Introduction

  • The solver targets general, including non-symmetric, sparse systems by approximating large frontal matrices with HSS and using randomized sampling for compression.The earlier algorithm in [55] handles only symmetric positive definite systems; this implementation uses the randomized method from [39].Source: §1 · p. 1
  • HSS is a hierarchical rank-structured format that approximates off-diagonal blocks as low-rank and exploits hierarchically nested generator bases.Unlike H, H2, and BLR, HSS and HODLR approximate only off-diagonal blocks; nested bases can enable asymptotically faster factorization for some matrices.Source: §1 · p. 2
  • Randomized HSS construction avoids explicitly forming larger frontal matrices as dense, while the randomized method with fast matrix-vector products can have linear or almost linear complexity, depending on rank structure.The method combines interpolative decompositions with the ULV decomposition from.Source: §1 · p. 2
  • Adaptive HSS construction computes the maximum submatrix rank during compression instead of assuming it in advance, and task parallelism supports performance and scalability on multicore and many-core systems.The code uses OpenMP task scheduling; recursively generated tasks mean the DAG is never explicitly constructed.Source: §1 · p. 2
  • The implementation is fully structured and uses randomized sampling to mitigate the HSS extend-add obstacle faced by partially structured approaches.The earlier Hsolver was limited to regular-mesh discretizations and was only partially structured.Source: §1 · pp. 2–3
  • The code is parallel, unlike the algebraic HSS solver of Napov and Li, and reports detailed results using randomized-sampling HSS compression.The earlier experiments did not include randomized-sampling HSS compression.Source: §1 · p. 3

2 HSS: Hierarchically Semi-Separable matrices

The section introduces HSS matrices as a data-sparse representation and describes their randomized construction, adaptive rank selection, and ULV-like factorization. The methods use sampled matrix actions and interpolative decompositions to build and solve with HSS representations.

2.1 Overview of the HSS matrix format

  • HSS stores dense diagonal leaf blocks and approximates sibling off-diagonal blocks with low-rank factors whose bases are nested through the tree.Its storage and matrix-vector multiplication costs are O(rN), compared with O(N^2) for dense multiplication.Source: §2.1 · pp. 3–4

2.2 Fast HSS construction through randomized sampling

  • Randomized HSS construction compresses off-diagonal blocks from samples AR and A*R, using interpolative decompositions to form the representation.The sampling approach needs a fast matrix-vector routine and selected matrix elements rather than explicit access to every entry.Source: §2.2 · pp. 4–5

2.2.1 Adaptive scheme to determine the HSS-rank

  • When the HSS rank is unknown, the algorithm increases the sample width and reuses successful descendant compressions while updating their new sample columns.The procedure repeats when a node fails to meet the requested accuracy.Source: §2.2.1 · p. 5

2.2.2 Implementation issues

  • The implementation uses rank-revealing QR in the interpolative decomposition; replacing it with strong rank-revealing QR is left as future work for possible accuracy and rank gains.The passage also names adaptive cross approximation and a matrix-free approach as alternatives.Source: §2.2.2 · p. 6

2.3 ULV-like factorization and solve

  • The paper uses a ULV-like decomposition that exploits the special structure of the HSS generators rather than relying on unitary transformations.The factorization is adapted from the ULV methods for symmetric and non-symmetric matrices described in prior work.Source: §2.3 · p. 6

3 Multifrontal sparse LU factorization

The section describes sparse LU as a sequence of partial factorizations and Schur-complement updates organized by an elimination tree. It covers preprocessing, frontal-matrix assembly and elimination, and the traversals used to solve the factored system.

3 Multifrontal sparse LU factorization

  • Multifrontal LU reduces sparse factorization to partial factorizations of smaller dense matrices and Schur-complement updates, organized by an elimination tree.Preprocessing scales and permutes the matrix for stability, then applies a fill-reducing permutation based on nested dissection.Source: §3 · p. 7, §3.1 · p. 7

3.2 Numerical factorization

  • A node's frontal matrix is assembled from the corresponding matrix entries and its children's extended update matrices, then its fully summed variables are eliminated and a Schur-complement update is formed.Child updates with different index sets require alignment by padding with zeros during extend-add.Source: §3.2 · pp. 7–8

3.3 Solution

  • After factorization, forward substitution traverses the elimination tree bottom-up and backward substitution traverses it top-down.The forward and backward steps solve triangular systems with L and U, respectively.Source: §3.3 · p. 8

4 Multifrontal solver with HSS frontal matrices

The section embeds HSS compression and factorization in a multifrontal solver to reduce complexity and storage. It describes frontal-matrix selection, separator ordering, randomized extend-merge, HSS factorization and solve-related procedures, and complexity considerations.

4.1 Selection of HSS frontal matrices

  • A switch level determines which frontal matrices remain dense and which are HSS-compressed; performance depends crucially on its choice.The proposed analysis-based rule balances factorization costs above and below the switch, but the passage says that rule is not very practical.Source: §4.1 · p. 8

4.2 Separator reordering

  • Alternative reorderings to nested dissection are possible, and the influence of reordering on off-diagonal ranks is studied in cited work.The supplied text does not state which reordering produces the best ranks.Source: §4.2 · p. 9

4.3 Skinny extend-add

  • For HSS fronts, child random matrices are merged at the parent and missing entries are generated; row-wise seeding by global index keeps sibling random entries consistent.This randomized construction lets the fronts above the switch level remain fully HSS-compressed without being explicitly formed as dense matrices.Source: §4.3 · p. 9

4.4 ULV factorization and low-rank Schur complement update

  • After HSS compression, the solver uses ULV for an HSS frontal matrix and computes the low-rank Schur-complement update through HSS matrix-vector multiplication and dense rectangular products.Dense fronts instead use pivoted LU and explicitly form the update.Source: §4.4 · p. 10

4.5 Extracting elements from an HSS matrix

  • Selected columns of an HSS matrix can be extracted by multiplying it by unit vectors, with parts of the tree traversal pruned when possible.At leaf nodes, the required columns can be selected directly from V*.Source: §4.5 · p. 10

4.6 Preconditioning versus iterative refinement

  • The HSS multifrontal method is used as a GMRES preconditioner rather than with iterative refinement, and increasing compression tolerance can increase outer iterations.For the same number of multifrontal solve steps, the passage reports that a Krylov solver typically gives smaller residuals than iterative refinement.Source: §4.6 · p. 10

4.7 Solver complexity

  • Standard multifrontal complexity is dominated by the largest dense fronts, while HSS-solver complexity depends on compression costs and the frontal-matrix rank pattern.The section summarizes complexity for two PDE types and two sparse solvers in a table not included among the supplied passages.Source: §4.7 · p. 10

5 Shared memory parallel implementation

The shared-memory implementation uses OpenMP tasks to coordinate parallel work across the elimination tree, HSS hierarchy, and dense operations within tree nodes. It shifts from tree-level to in-node parallelism as concurrency in the trees decreases.

5 Shared memory parallel implementation

  • The shared-memory code uses C++ and OpenMP tasks, with BLAS, LAPACK, and graph-partitioning libraries as dependencies.OpenMP was chosen for its ease of use, performance, and support; the authors note that similar task-parallel runtimes could also be used.Source: §5 · pp. 10–11

5.1 Task based tree parallelism

  • Bottom-up tree traversals expose parallelism across independent subtrees, with task creation capped at depth log2(number of threads) + 3 to limit overhead.At that setting, the maximum number of live tasks is about 8 times the thread count.Source: §5.1 · p. 11

5.2 Hybrid node and tree parallelism

  • As tree-level parallelism declines near the root, task parallelism within tree nodes lets the runtime shift work toward dense algebra while balancing cores.The approach avoids calling multithreaded BLAS or LAPACK from multiple tasks, which can cause oversubscription.Source: §5.2 · pp. 11–12, Figure 3 · p. 13

5.3 Parallel BLAS and LAPACK

  • Recursive task-based BLAS splits matrix multiplications into smaller operations, stopping at a tuning threshold to call vendor-optimized sequential GEMM.For rank-revealing QR, the modified routine stops once the ε-rank is detected rather than computing the full decomposition.Source: §5.3 · pp. 12, 14
  • Recursive LU scales better than reference netlib getrf with parallel BLAS, but worse than MKL optimized multithreaded getrf, while retaining subtree concurrency.Calling MKL multithreaded getrf from multiple threads simultaneously risks oversubscription and a performance penalty.Source: Figure 4 · p. 14, §5.3 · p. 14

5.4 Scaling bottlenecks

  • Memory allocation can bottleneck parallel symbolic factorization, and not all compilers support the latest OpenMP 4.0 standard.The authors report improved performance with scalable allocators such as TCMalloc or TBBMalloc.Source: §5.4 · pp. 14–15
  • On a 60-core Intel Xeon Phi, symbolic factorization runs up to 56× faster with TBBMalloc than with the default allocator.The parallel symbolic factorization can make memory allocation a scaling bottleneck.Source: §5.4 · pp. 14–15

6 Numerical experiments

The experiments benchmark HSS-augmented multifrontal solving on regular-grid PDEs and application matrices, and assess its parallel scaling.

6.1 PDEs on a regular grid

  • For 2D Poisson, HSS gives about 2× factorization speedup at 7 levels, while 3D permits more aggressive compression and benefits more overall.For 50002 unknowns, 2D HSS reduces flops by about 3.3×, but increased solve work and lower flop rates limit time savings; for 125^3 3D Poisson at ε = 0.9, flops fall to 4.4% of MF. Parameters were tuned for one solve, and accurate repeated solves may favor exact MF.Source: §6.1 · p. 16, Figure 5 · p. 17, Figure 6 · p. 18
  • On the tested 2D and 3D PDE grids, HSS compression reduces the predicted flop-scaling exponent, with larger constants than pure multifrontal solving.The fitted exponents are sensitive to the data and not very reliable.Source: §6.1 · p. 16, Figure 6 · p. 18

6.2 Matrices from various applications

  • Application-matrix benchmarks compare MF+HSS, pure MF, and PARDISO on a 12-core Ivy Bridge system; the tested matrices come from 2D or 3D PDEs.The regular-grid benchmarks use geometric nested dissection; the application matrices come from the Florida Sparse Matrix Collection and SciDAC projects.Source: §6.1 · p. 15, §6 · p. 15, §6.2 · p. 20
  • MF+HSS consistently outperforms pure MF on the selected application matrices and outperforms PARDISO; for dt0, pure MF ran out of memory in double precision.For A22, separator reordering takes substantial time; omitting link-two edges instead makes HSS ranks too large to yield a net performance benefit.Source: §6.2 · p. 20, Figure 9 · p. 21, §6.3 · p. 21

6.3 Many-core parallel performance

  • MF+HSS scales well for numerical factorization and reasonably for solving on Ivy Bridge and Xeon Phi, while sequential reordering becomes a thread-scaling bottleneck.The Xeon Phi's 8 GB memory excludes the larger Table 3 problems.Source: §6.3 · pp. 20, 22, Table 3 · p. 21, Figure 10 · p. 22

7 Conclusions & Outlook

The work presents an initial HSS-structured multifrontal solver, reporting speedups up to 7× and favorable comparisons with PARDISO. Its lower computational complexity comes with larger constants, which the authors identify as a focus for improvement.

7 Conclusions & Outlook

  • The solver reports speedups up to 7× over pure multifrontal solving and compares favorably with PARDISO, but has larger constants and lower floating-point throughput.The authors identify reducing HSS ranks as one possible route to lower constants.Source: §7 · p. 22
  • The authors plan to reduce HSS ranks and address larger problems with a distributed-memory implementation.Possible rank-reduction strategies include power iteration, improved separator reordering, and stronger rank-revealing QR.Source: §7 · p. 22
  • STRUMPACK currently combines a sparse shared-memory solver with a dense distributed-memory solver; the longer-term goal is one scalable hybrid-parallel code for both.The authors describe this work and the distributed HSS code developed for as steps toward that goal.Source: §7 · p. 23
  • Nested HSS approximation and matrix-free direct-solver preconditioners are identified as possible future directions.The authors note that a matrix-free direct-solver-based preconditioner may require knowledge of the sparsity pattern.Source: §7 · p. 23
Loading 1502.07405v1…