🔬 materials science

Solvers for the Hermitian and the pseudo-Hermitian Bethe-Salpeter equation in the Yambo code: Implementation and Performance

This paper presents and benchmarks the implementation of both direct (ScaLAPACK/ELPA) and iterative (SLEPc) solvers for Hermitian and pseudo-Hermitian Bethe-Salpeter equations within the Yambo code, demonstrating the feasibility of handling dense matrices of up to $10^5$ elements on CPU and GPU clusters.

Original authors: Petru Milev, Blanca Mellado-Pinto, Muralidhar Nalabothula, Ali Esquembre Kucukalic, Fernando Alvarruiz, Enrique Ramos, Francesco Filippone, Alejandro Molina-Sanchez, Ludger Wirtz, Jose E. Roman, David

Published 2026-07-07

📖 5 min read🧠 Deep dive

CC BY 4.0

Original authors: Petru Milev, Blanca Mellado-Pinto, Muralidhar Nalabothula, Ali Esquembre Kucukalic, Fernando Alvarruiz, Enrique Ramos, Francesco Filippone, Alejandro Molina-Sanchez, Ludger Wirtz, Jose E. Roman, Davide Sangalli

Original paper licensed under CC BY 4.0 (http://creativecommons.org/licenses/by/4.0/). ✨ This is an AI-generated explanation of the paper below. It is not written or endorsed by the authors. For technical accuracy, refer to the original paper. Read full disclaimer

Imagine you are trying to predict how a new material will interact with light—whether it will glow, absorb color, or act as a mirror. To do this, scientists use a complex mathematical recipe called the Bethe-Salpeter Equation (BSE). Think of this equation as a massive, intricate puzzle where every piece represents an electron and a "hole" (a missing electron) interacting with each other.

The paper you provided is essentially a report card on how fast and efficiently a specific software tool, called Yambo, can solve these puzzles. The authors tested two different ways to solve the puzzle and compared them on different types of computer hardware (standard processors and powerful graphics cards).

Here is a breakdown of their findings using simple analogies:

1. The Two Types of Puzzles

The BSE puzzle comes in two flavors:

The "Simple" Version (Hermitian/TDA): This is like a puzzle where the pieces only interact in one direction. It's easier to solve.
The "Complex" Version (Pseudo-Hermitian): This is the full, realistic version where pieces interact in both directions, creating a much more tangled web. Historically, this was much harder and slower to solve.

2. The Two Solvers (The "How-To" Methods)

The researchers tested two main strategies to solve these puzzles:

The "Brute Force" Method (Direct Diagonalization):
Imagine trying to solve a jigsaw puzzle by laying out every single piece on a giant table, sorting them all at once, and finding the picture. This is very thorough and accurate, but if the puzzle gets too big (millions of pieces), the table runs out of space, and the process takes forever.
- The Tools: They used libraries called ScaLAPACK and ELPA. Think of these as different brands of "puzzle sorters." The paper found that ELPA is a much faster, more efficient sorter than ScaLAPACK, especially for the complex version of the puzzle.
The "Smart Search" Method (Iterative Solver):
Instead of looking at every single piece, imagine you only need to find the 100 most important pieces to see the main picture. You use a smart search algorithm to find just those pieces without sorting the whole table. This is much faster and uses less space.
- The Tool: They used a library called SLEPc. This is like a detective that only looks for the specific clues you need.

3. The Big Breakthrough: Fixing the "Complex" Puzzle

The most exciting part of the paper is how they handled the "Complex" (Pseudo-Hermitian) version.

Before: Solving the complex version was like trying to run a marathon in heavy boots. It was slow and inefficient.
Now: The authors implemented a new trick that recognizes the hidden symmetry in the complex puzzle. It's like realizing the heavy boots are actually made of a special material that makes you run faster than before.
The Result: By using this new trick, solving the complex puzzle is now almost as fast as solving the simple one. They managed to solve puzzles with 100,000 pieces (matrices of size $10^5$ ), which was previously very difficult.

4. The Hardware: CPUs vs. GPUs

They tested these methods on two types of computer "muscle":

CPUs (The General Workers): These are like a team of many generalists working together. The "Brute Force" method (ELPA) worked very well here, scaling up efficiently as they added more workers.
GPUs (The Super-Specialists): These are like a massive army of speedsters (graphics cards).
- For small puzzles, the speedsters were so fast that adding more of them didn't help much (they were waiting for instructions).
- For huge puzzles, the speedsters shined. The larger the puzzle, the better the GPUs performed, allowing them to handle massive calculations that would choke a standard computer.

5. The Memory Problem (The "Table Space")

A major hurdle in these calculations is memory (RAM).

The Issue: Storing the full puzzle table requires a massive amount of space. Sometimes, just storing the data takes up terabytes (thousands of gigabytes).
The Finding: The "Smart Search" method (SLEPc) is much better at saving space. It doesn't need to lay out the whole table; it just grabs what it needs. However, the "Brute Force" method (ELPA) is so efficient at organizing the data that it can still handle huge puzzles if you have enough computer nodes (workers) to share the load.

Summary

The paper claims that by combining the Yambo software with powerful external libraries (ELPA and SLEPc) and using a new mathematical trick to handle complex interactions, scientists can now solve massive, realistic models of how materials interact with light.

They have successfully removed the "speed barrier" that used to stop them from studying very large or complex materials. They can now handle puzzles with 100,000 pieces on both standard supercomputers and next-generation GPU clusters, making it feasible to predict optical properties for complex 2D materials, defects, and nanostructures with high accuracy.

Technical Summary: Solvers for the Hermitian and Pseudo-Hermitian Bethe-Salpeter Equation in the Yambo Code

Problem Statement
The Bethe-Salpeter Equation (BSE) is the state-of-the-art first-principles approach for computing neutral excitations and optical properties in condensed matter physics. However, solving the BSE eigenvalue problem presents significant computational challenges, particularly as system sizes increase. The BSE matrix, often reaching dimensions of $N \approx 10^5$ , is dense and structured. While the Tamm-Dancoff approximation (TDA) reduces the problem to a Hermitian eigenvalue problem, many physical cases (e.g., magnetic materials, systems without time-reversal symmetry) require the full pseudo-Hermitian formulation. This full form involves a non-Hermitian structured matrix that is computationally more demanding to diagonalize. Furthermore, while matrix generation has been optimized for High-Performance Computing (HPC) environments, the solver step has historically been a bottleneck, especially when eigenvectors (excitonic wavefunctions) are required. Standard iterative solvers often struggle with the specific structure of the full BSE matrix, and direct diagonalization scales as $O(N^3)$ , becoming prohibitive for large systems.

Methodology
This work analyzes and implements two distinct strategies for solving the BSE eigenvalue problem within the Yambo code, interfacing with established linear algebra libraries:

Direct Diagonalization: The authors interface Yambo with ScaLAPACK and ELPA libraries. For the full pseudo-Hermitian case (where coupling terms are non-zero), they implement a specialized algorithm that transforms the problem into a real skew-symmetric eigenvalue problem. This involves constructing a real symmetric matrix $M_r$ and utilizing a unitary transformation to map the original Hamiltonian $H_{2p}$ to a form solvable by standard Hermitian or skew-symmetric solvers. This approach avoids the computational cost of general non-Hermitian solvers.
Iterative Solver: The authors implement an iterative approach using the SLEPc library. They utilize structure-preserving Lanczos iterations (specifically the Shao method) designed for pseudo-Hermitian matrices. Two distribution schemes are tested:
- Shell approach: Reuses Yambo's existing parallel structure without explicit matrix redistribution, minimizing memory but lacking GPU support.
- Explicit approach: Redistributes the BSE matrix from Yambo's block-cyclic layout to SLEPc's standard 1D row-wise partitioning, enabling GPU acceleration via CUDA/HIP.

The performance of these solvers is evaluated on both CPU (ISM-CNR cluster) and GPU (Leonardo Booster) architectures. The study considers three matrix sizes ( $N \approx 10^4$ , $3 \cdot 10^4$ , and $10^5$ ) and compares Hermitian (TDA), non-Hermitian, and pseudo-Hermitian algorithmic implementations.

Key Contributions

Implementation of Pseudo-Hermitian Solvers: The paper details the implementation of specialized solvers for the full BSE matrix that exploit its pseudo-Hermitian structure. This allows the use of efficient real skew-symmetric or Hermitian-definite solvers rather than general non-Hermitian ones.
Library Integration: The work establishes robust interfaces between Yambo and ScaLAPACK, ELPA, and SLEPc, including a new library (Lumen fork) that handles the specific diagonalization routines required for the coupling case.
GPU Porting: The authors successfully ported the iterative solver (SLEPc explicit) and the diagonalization solver (ELPA) to GPU architectures, demonstrating feasibility on NVIDIA A100 GPUs.
Performance Benchmarking: Comprehensive benchmarks are provided for time complexity, strong scaling, and memory utilization across different matrix sizes and architectures.

Results

Algorithmic Efficiency: The new pseudo-Hermitian (PH) algorithms significantly outperform standard non-Hermitian (NH) approaches. For iterative solvers (SLEPc), the PH implementation yields a 5- to 6-fold speedup compared to NH methods, primarily due to improved scaling exponents (closer to the ideal $O(N^2)$ ). For direct diagonalization, the PH approach also offers a roughly 5x speedup over NH methods.
Scaling Behavior:
- CPU: ELPA consistently demonstrates superior scaling and efficiency compared to ScaLAPACK, particularly for the coupling case. The ELPA parallelization factor ( $p$ ) approaches the ideal $-1$ (values around $-0.8$ to $-0.9$), whereas ScaLAPACK and SLEPc show lower efficiency ( $p \approx -0.6$ to $-0.7$).
- GPU: SLEPc is orders of magnitude faster than ELPA for small to medium matrices but is limited by memory overhead and data transfer. ELPA shows excellent scaling on GPUs for the coupling case, with efficiency plateauing only at very large core counts.
Memory Usage: Iterative solvers (SLEPc) generally require less memory than direct diagonalization. However, the explicit SLEPc approach incurs a significant parallel overhead memory cost when transitioning from serial to parallel runs. The pseudo-Hermitian implementation in SLEPc improves memory efficiency compared to the non-Hermitian version, while the diagonalization case shows a slight increase in memory usage for PH compared to NH.
Feasibility: The results demonstrate that dense BSE matrices of the order of $10^5$ can now be handled effectively, provided sufficient MPI tasks and memory resources are available.

Significance and Claims
The authors claim that their work removes a critical barrier in BSE calculations: the "solver barrier" for large matrices. By leveraging the pseudo-Hermitian structure, they have made the solution of the full coupling BSE matrix nearly as efficient as the resonant (TDA) approximation. This enables the study of excitonic effects in complex systems (e.g., magnetic materials, 2D heterostructures) where the TDA is insufficient.

The paper emphasizes that while matrix generation has seen significant optimization, the solver step is now the limiting factor for extracting eigenvectors. The authors argue that using external, highly optimized libraries (ELPA, ScaLAPACK, SLEPc) is a viable and effective strategy for BSE solvers, contrasting with DFT where custom iterative solvers are often preferred. The successful GPU porting and the demonstration of handling $N \approx 10^5$ matrices suggest that the code is now capable of addressing large-scale optical property simulations that were previously computationally infeasible. The work does not claim to solve the memory storage issue for the largest systems (which can exceed 1 TB) but demonstrates that the diagonalization step can be performed efficiently on such scales.