PaperlayerView exact revision on arXiv

Source-linked AI summary

Accelerating Atom Simulations with Variable-Block Sparse Matrix Library

Zhanghao Zhouyin, Hong Guo

arXiv:2609.04397v1cond-mat.mtrl-scics.DCquant-ph

TL;DR

Localized-orbital atomistic operators have variable-shaped blocks, while scalar CSR expands their entries and obscures that structure. VBCSR maps operator types to structure-matched sparse backends and batches equal-shaped blocks for dense kernels. It is competitive with or faster than tested references and is demonstrated on an InP nanoparticle exceeding one million atoms.

  • Problem

    Conventional scalar CSR expands atom-pair blocks into individual entries, obscuring their structure, while uniform-block storage pads smaller interactions to the largest block size.

    Source: §I · p. 1

  • Method

    VBCSR automatically maps scalar, uniform-basis, and species-dependent operators to CSR, BSR, or VBCSR, grouping equal-shaped blocks for optimized dense kernels.

    Source: §I · p. 1, B. Three Acceleration Mechanisms · p. 2

  • Results

    Across scalar, fixed-block, and variable-block operators, VBCSR was competitive with or faster than tested references; variable-block sparse matrix-matrix multiplication reached up to 8.4× with one thread and 3.8× with 16 threads.

    Source: §I · p. 1

  • Takeaways & Limitations

    The benchmarks demonstrate these backends on the tested CPU platform, while the InP application demonstrates capacity for large atomistic systems.

    Source: §VI · p. 7

    The reported speedups are specific to the tested operators, block-size range, thread counts, and references, and do not establish general superiority.

    Source: A. Kernel Efficiency · p. 5

Abstract

from arXiv · show

Modern atomistic simulations increasingly employ localized orbitals to represent quantum operators, yielding sparse block matrices whose block shapes vary with chemical species and basis choice. Conventional scalar sparse formats store the entries of each block individually, obscuring this local structure and limiting the use of efficient block algorithms. We present VBCSR, a distributed sparse matrix library that preserves variable-size atomic blocks and accelerates the core linear algebra of large-scale atomistic simulations. A unified interface automatically maps scalar, uniform-basis, and multispecies operators to compressed sparse row (CSR), block sparse row (BSR), or variable-block compressed sparse row (VBCSR). Our advanced acceleration method groups blocks of equal shape and dispatches them to optimized dense kernels. In the reported benchmarks, VBCSR outperforms the tested Python-accessible reference implementations for several block-sparse benchmarks. We further demonstrate VBCSR in an InP nanoparticle application containing more than \(10^6\) atoms.

I. INTRODUCTION

Localized-orbital atomistic simulations produce sparse operators whose block structure is obscured by scalar storage or padded by uniform-block storage. VBCSR preserves that structure through a unified distributed interface, with competitive benchmark results and an InP demonstration exceeding one million atoms.

I. INTRODUCTION

  • Scalar CSR expands dense atom-pair couplings into individual entries, while uniform-block storage pads smaller interactions to the largest block size.These representation costs grow with system size and can restrict the scale of simulations.Source: §I · p. 1
  • Atomistic developers need a lightweight distributed operator layer that preserves species-dependent blocks, scales to large simulations, and embeds in specialized workflows.Without such a layer, applications devote substantial effort to managing quantum-operator matrices instead of implementing their physical models.Source: §I · p. 1
  • VBCSR uses an atom-centred distributed graph and automatically selects CSR, BSR, or VBCSR for scalar, uniform-basis, or species-dependent operators.Equal-shaped blocks are grouped for optimized dense kernels.Source: §I · p. 1
  • VBCSR was competitive with or faster than tested references across scalar, fixed-block, and variable-block operators, with variable-block sparse matrix-matrix multiplication accelerated by up to 8.4× with one thread and 3.8× with 16 threads.The InP nanoparticle application computes density of states and size-dependent band gaps for systems containing up to 1,063,609 atoms.Source: §I · p. 1

II. VBCSR FOR ATOMISTIC SIMULATION

Localized atomic orbitals produce sparse operators with dense atom-pair blocks, whose dimensions vary with the atoms’ orbital counts. VBCSR matches storage and execution to that structure through backend selection, shape-based batching, and graph-guided parallel work.

A. From Localized Atomic Orbitals to Block-Sparse Operators

  • Localized-orbital interactions form dense blocks of size b_i × b_j; multispecies systems therefore produce variable block shapes that VBCSR retains without scalar expansion or maximum-size padding.Scalar operators suit CSR, uniform-sized blocks suit BSR, and heterogeneous blocks form the VBCSR structure.Source: A. From Localized Atomic Orbitals to… · p. 2

B. Three Acceleration Mechanisms

  • VBCSR automatically selects CSR, BSR, or VBCSR according to orbital counts, then groups equal-shaped blocks into contiguous pages for SIMD or BLAS kernels.The routing layer dispatches dense block operations to a specialized native SIMD kernel or a BLAS routine according to shape.Source: B. Three Acceleration Mechanisms · p. 2
  • Graph-derived partitions and exchange schedules are reused across operators, while OpenMP block batches and NUMA first-touch allocation keep threaded work close to data.These choices reduce avoidable representation changes, expose regular dense work, and keep computation close to data.Source: B. Three Acceleration Mechanisms · p. 2

C. Python Simplicity, Optimized Execution

VBCSR exposes distributed sparse operations through a single Python matrix object, while hiding backend selection, ghost synchronization, and threaded block execution. Conversion and solver interfaces support integration with existing Python scientific workflows.

C. Python Simplicity, Optimized Execution

  • A single Python matrix object supports construction, assembly, vector or multivector application, arithmetic, filtering, transposition, and matrix-function utilities.Matrix application and sparse multiplication use the @ operator.Source: C. Python Simplicity, Optimized… · p. 2
  • In serial mode, VBCSR interoperates with SciPy solvers through LinearOperator and converts to and from NumPy arrays and SciPy sparse matrices.These paths support integration with existing atomistic codes.Source: C. Python Simplicity, Optimized… · p. 2
  • A single Python expression can trigger ghost synchronization, backend selection, and threaded block kernels without requiring users to coordinate those steps.The Python layer keeps workflows operator-level while using specialized compiled execution.Source: C. Python Simplicity, Optimized… · p. 2

A. Data Model and Distributed Ownership

VBCSR’s distributed graph records block sizes, ownership, connectivity, and ghosts, making the matrix’s partition and communication boundary explicit. Reusable graph-derived schedules support forward synchronization and adjoint reduction.

A. Data Model and Distributed Ownership

  • Each graph node records its scalar block size, and each MPI rank stores owned block rows plus locally represented remote ghost nodes.Owned nodes precede ghosts, which are ordered by owner rank and then global identifier to make ownership explicit and communication buffers rank-contiguous.Source: A. Data Model and Distributed Ownership · p. 2
  • The graph builds reusable send and receive schedules for matrix operators, DistVector, and DistMultiVector, with matrix-vector multiplication synchronizing ghost entries.The schedule follows node connectivity in the distribution pattern.Source: A. Data Model and Distributed Ownership · p. 2
  • Adjoint operations return contributions accumulated for ghost columns to their owners for reduction, reusing the graph-derived communication boundary.The figure presents a shared distributed graph and unified interface across CSR, BSR, and VBCSR backends, with MPI communication, thread parallelism, and dense kernels.Source: A. Data Model and Distributed Ownership · p. 3, Figure 1 · p. 3

B. Structure-Matched Sparse Algorithms

VBCSR matches storage and kernels to block structure, routing compatible shapes to dense operations. Its distributed matrix products combine ghost communication with shape-aware computation and optional block-threshold filtering.

B. Structure-Matched Sparse Algorithms

  • VBCSR stores one active CSR, BSR, or VBCSR backend; its variable-block backend maps heterogeneous blocks into shape-specific buffers for direct batching without repacking.Construction registers each distinct block shape and maps each heterogeneous block to its shape-specific buffer.Source: B. Structure-Matched Sparse Algorithms · p. 3
  • Shape-based routing uses vectorized specialized kernels for small matrices and linked BLAS routines for larger blocks, with OpenMP distributing batches among threads.Output construction follows the same thread partitions, and threads touch pages in their corresponding NUMA domains.Source: B. Structure-Matched Sparse Algorithms · p. 3
  • Forward matrix products synchronize ghost inputs before traversing owned block rows; adjoint products reduce ghost-column contributions to their owning ranks.The stored block columns in row i are denoted J(i).Source: B. Structure-Matched Sparse Algorithms · p. 3
  • Thresholded multiplication exchanges remote row metadata and block norms, then rejects a candidate when its accumulated upper bound is below the nonnegative threshold ε.Filtering is applicable when physical locality or matrix decay makes sufficiently small blocks negligible within the calculation’s accuracy requirements.Source: B. Structure-Matched Sparse Algorithms · p. 4
  • The surviving symbolic pattern determines the result graph and needed remote blocks; numerical products are grouped by row, inner, and column dimensions before accumulation and final block-norm filtering.Grouping yields compatible dense products of shape (bi × bk)(bk × bj).Source: B. Structure-Matched Sparse Algorithms · p. 4

C. Atomistic Interfaces

AtomicData and ImageContainer connect species-dependent atom data and periodic translations to the distributed block-graph operator model. Reciprocal-space sampling combines translation-resolved images into a complex block-sparse operator.

C. Atomistic Interfaces

  • AtomicData records each atom’s identifier, type, position, and species-dependent orbital count, which becomes its node block size.Distributed neighbour construction assigns owned atoms to ranks and records adjacent local or remote atoms, with periodic edges carrying lattice shifts.Source: C. Atomistic Interfaces · p. 4
  • ImageContainer groups periodic real-space operator edges by lattice translation, with each translation represented by a block-sparse image following the atom partition and block-size model.The union of image graphs defines AtomicData connectivity.Source: C. Atomistic Interfaces · p. 4
  • Reciprocal-space sampling combines the images into a complex block-sparse operator; R is an integer lattice translation and k uses fractional reciprocal coordinates.Equation (4) is the default R convention; the alternative R+τ convention phases block (i, j) using R + τj − τi.Source: C. Atomistic Interfaces · p. 4

IV. BENCHMARK EXPERIMENTS

The experiments evaluate kernel efficiency against references and assess how operator families scale across MPI ranks. They test three sparse structural domains and operations under a fixed approximate storage budget.

IV. BENCHMARK EXPERIMENTS

  • The benchmark asks whether VBCSR retains or improves kernel efficiency against Python-accessible and optimized references, and how operator families scale across MPI ranks.The scaling question includes distribution and ghost communication.Source: §IV · p. 4
  • Tests cover CSR, BSR, and VBCSR on atom-like periodic geometric graphs with finite cutoffs and comparable mean connectivity, using SpMV, SpMM, and SpGEMM in double precision.Experiments used AMD EPYC 7532 nodes with 32 cores per node, Open MPI and OpenMP, and MKL with customized BLAS kernels.Source: §IV · p. 4
  • Kernel-efficiency and strong-scaling graph sizes targeted approximately 500 MiB of matrix storage, and each kernel time was the median of seven repetitions after three warm-ups.The graph size was chosen separately for each structural domain.Source: §IV · p. 4

A. Kernel Efficiency

VBCSR is competitive with or faster than the tested references across scalar, fixed-block, and variable-block operators. The largest reported gains are in variable-block matrix products, with results bounded to the tested configurations and references.

A. Kernel Efficiency

  • SpMV remained relatively memory-bound, while SpMM allowed greater matrix-value reuse and SpGEMM created more opportunities to batch products of equal shape.These qualitative trends were reported as consistent with the described acceleration mechanisms.Source: A. Kernel Efficiency · pp. 4–5
  • Variable-block speedups relative to the fastest reference reached approximately 8.4× for one-thread SpGEMM; with 16 threads, speedups were approximately 3.3 (SpMV), 3.2 (SpMM), and 3.8 (SpGEMM).Figure 2 reports median kernel times and speedup labels for the tested operations and structural domains.Source: Figure 2 · p. 5, A. Kernel Efficiency · p. 5
  • At 16 threads, scalar CSR speedups ranged from approximately 1.4 to 2.0 depending on the operation, while one-thread cases were near parity with references.The authors report these results as indicating no large backend-dispatch penalty for the tested scalar operators.Source: A. Kernel Efficiency · p. 5
  • For block size 8, one-thread BSR SpMV and SpMM were approximately 2.1–2.5 times faster than references, with 16-thread speedups of approximately 1.3–1.6.The measurements are described as consistent with fixed-block dense kernels, SIMD-aware batched operations, and NUMA-aware memory placement.Source: A. Kernel Efficiency · p. 5
  • The observed advantage beyond SpMV is consistent with preserving heterogeneous blocks, avoiding scalar CSR representation, and grouping compatible products by shape.SpMM reuses operator data across right-hand sides, while SpGEMM presents shape-batched block products to dense kernels.Source: A. Kernel Efficiency · p. 5
  • The reported speedups are specific to the tested operators, block-size range, thread counts, and references, and are not evidence of general superiority.The authors characterize them as support for VBCSR’s competitive performance.Source: A. Kernel Efficiency · p. 5

B. Distributed Scaling and Reproducibility

The distributed tests measure strong scaling with a fixed global matrix and weak scaling with a fixed local matrix footprint. Strong-scaling performance depends on operation arithmetic intensity, while weak-scaling efficiency remains close to ideal over the tested range.

B. Distributed Scaling and Reproducibility

  • Strong scaling increased MPI ranks from 1 to 64 at fixed global problem size, whereas weak scaling held the local matrix footprint approximately fixed as the global problem grew with rank count.Figure 3 defines speedup as S(p) = T(1)/T(p) and weak-scaling efficiency as Ew(p) = T(32)/T(p).Source: B. Distributed Scaling and… · p. 5, Figure 3 · p. 6
  • Strong scaling was good across all three structural domains, with SpGEMM generally scaling further than SpMV and SpMM as communication, synchronization, and scheduling overheads grew more prominent.The fixed global matrix storage target was approximately 500 MiB per structural domain.Source: B. Distributed Scaling and… · pp. 5–6, Figure 3 · p. 6
  • Weak-scaling efficiencies ranged from roughly 0.92 to 1.06 as ranks increased from 32 to 128, while runtime changed by less than 10% as the global problem quadrupled.The local matrix footprint remained approximately fixed.Source: B. Distributed Scaling and… · p. 6

V. LARGE-SCALE INP NANOPARTICLE DEMONSTRATION

The InP demonstration applies VBCSR to hydrogen-passivated nanoparticles spanning three sizes, including a completed particle with more than one million atoms. The calculated gap shifts agree with a previously reported size-dependent fit.

V. LARGE-SCALE INP NANOPARTICLE DEMONSTRATION

  • The largest completed hydrogen-passivated particle contains 1,063,609 atoms and has an effective diameter of 36.44 nm.The three systems contain 12,083, 120,250, and 1,063,609 atoms, at diameters of 7.71, 17.28, and 36.44 nm.Source: §V · p. 6, Figure 4 · p. 7
  • The DOS calculation uses a second-nearest-neighbour In(sp3)–P(sp3d5)–H(s) tight-binding Hamiltonian and 147,456 KPM moments, giving 0.987 meV energy resolution.The gap is estimated from the low-DOS interval using a threshold of 10% of the resolution.Source: §V · pp. 6–7, Figure 4 · p. 7
  • The estimated InP nanoparticle gap shifts are consistent with the reported Sapra et al. fit.The gap shift is estimated from the KPM density of states relative to the bulk band gap.Source: §V · p. 7, Figure 4 · p. 7

VI. CONCLUSION

VBCSR provides a single interface for distributed atom-centred sparse operators, selecting storage and kernels suited to scalar, uniform-basis, and multispecies systems. Benchmarks on the tested CPU platform and the InP application demonstrate performance and scale, while finite-element PDE use remains future work.

VI. CONCLUSION

  • One interface supports CSR, BSR, and variable-block CSR backends for single-orbital, uniform-basis, and multispecies systems without changing application-level workflow.The atom-centred graph maps atoms to localized degrees of freedom and physical interactions to dense matrix blocks.Source: §VI · p. 7
  • The benchmarks demonstrate backend performance on the tested CPU platform, and the InP application demonstrates capacity for large atomistic systems.Source: §VI · p. 7
  • The graph-and-block abstraction may also support variable-space finite-element PDE discretizations, which the paper leaves for future exploration.Source: §VI · p. 7
Loading 2609.04397v1…