Week 6: Benchmarking the Native PETSc HEOM Solver

Following last week’s successful integration of the PETSc backend directly into QuTiP’s HEOMSolver, this week was entirely focused on benchmarking. The goal was to compare the newly integrated distributed PETSc backend against the traditional scipy.sparse-based solver using a standard Ising model spin chain.

The Benchmark Setup

We set up a 1D Ising model with N=4 spins and a hierarchy depth of max_depth = 6. We then tested the performance by coupling the spins to Drude-Lorentz baths.

The tests were run locally using 4 MPI ranks (mpiexec -n 4) for the PETSc backend, and compared against the single-threaded traditional solver.

The Results

Case 1: Single Bath (N=4, 1 Bath)

In this scenario, only the first spin is coupled to a bath. The hierarchy generates 462 ADOs. - Traditional Solver: ~97.4 seconds - PETSc Backend (4 ranks): ~145.8 seconds

Case 2: Two Baths (N=4, 2 Baths)

Here, the first two spins are coupled to baths, drastically increasing the state space to 8008 ADOs. - Traditional Solver: ~1166.5 seconds - PETSc Backend (4 ranks): ~1901.6 seconds

Note: During subsequent local runs with multiple baths, the system and terminal experienced freezing, likely due to heavy context switching or memory exhaustion on the local machine.

Analyzing the Performance

At first glance, the PETSc backend running on 4 cores is slower than the traditional solver running on a single core. However, this is an expected hurdle when moving to distributed computing for relatively “small” problems:

1. MPI Communication Overhead: Distributing a matrix of 462 (or even 8008) rows across 4 processes means the MPI ranks are spending more time communicating boundary data than actually crunching numbers. The network/IPC latency dominates the compute time.
2. Solver Stiffness and Overhead: PETSc’s implicit solvers (like bdf with bcgs/bjacobi preconditioners) have significant initialization and per-step algebraic overhead. For small matrices, traditional explicit ODE solvers (like QuTiP’s default adams or lsoda) fly through the integration steps because the entire matrix fits neatly into L2/L3 CPU cache.
3. Local Machine Constraints: Running 4 MPI ranks on a local laptop or desktop often leads to resource contention.

Conclusion and Next Steps

These results confirm that PETSc is functioning natively inside QuTiP, but they also remind us that MPI is not a magic bullet for all problem sizes. The true power of PETSc lies in weak scaling and solving problems that are simply impossible for the traditional solver due to memory limits (e.g., matrices with hundreds of thousands of ADOs).

For Week 7, the plan is to:

1. Move to an HPC Cluster: Benchmarking on a proper computing cluster to eliminate local resource contention.
2. Push the Limits: Test with 3 or 4 baths, or much larger N, to find the “crossover point” where the traditional solver exhausts RAM, but PETSc successfully distributes the load.
3. Tune PETSc Options: Experiment with different preconditioners, tolerances, and explicit TS methods in PETSc to minimize overhead for intermediate-sized problems.