ezp

Distributed linear solvers, distilled into one line of C++.

Header-only · C++20 · MPI + ScaLAPACK, without the boilerplate.

Documentation · Quick start · Solvers · Standalone solvers

ezp is a lightweight, header-only C++ wrapper around a curated set of distributed linear-system solvers. It hides the BLACS, ScaLAPACK, MPI and low-level C/Fortran minutiae behind an interface that feels just like solving a system sequentially.

solver.solve(A, B); // that's the whole idea

pardiso

Solver for general sparse matrices.

Why <tt>ezp</tt>?

Features

	Feature
✨	Ergonomic interface — prepare $A$ and $B$ on the root process; everything else is automatic
📦	Drop-in, header-only library — no build step, no ABI to manage; just `#include`
🧮	Full precision coverage — all four `DSZC` precisions, with 32-bit and 64-bit indexing
🖥️	Standalone solver binaries — invoke the solver from any `MPI`-capable caller, including Python via `mpi4py`
🧠	Build once, solve many — factorization is retained and reused across right-hand sides
✅	Validated with randomly generated inputs — an extensive test suite backs every solver

The library implements all precisions with both 32-bit and 64-bit integer indexing, and supports both row- and column-major process grids.

Supported Solvers

Dense (ScaLAPACK)

availability	type of matrix	operation	solver	package
✅	general (partial pivoting)	simple	PxGESV	ScaLAPACK
✅	general (partial pivoting)	expert	PxGESVX	ScaLAPACK
✅	symmetric/Hermitian positive definite	simple	PxPOSV	ScaLAPACK
✅	symmetric/Hermitian positive definite	expert	PxPOSVX	ScaLAPACK
✅	general band (partial pivoting)	simple	PxGBSV	ScaLAPACK
✅	general band (no pivoting)	simple	PxDBSV	ScaLAPACK
✅	symmetric/Hermitian positive definite band	simple	PxPBSV	ScaLAPACK

Sparse

availability	type of matrix	operation	solver	package
✅	sparse (one- or zero-indexing CSR format)	direct	PARDISO	MKL
✅	sparse (one-indexing COO format)	direct	MUMPS	MUMPS
✅	sparse (zero-indexing CSR format)	iterative	Lis	Lis

Expert drivers additionally provide equilibration, iterative refinement, condition estimation and determinant computation where applicable.

Quick Start

ezp requires a C++20-compatible compiler and the following drivers:

an implementation of LAPACK / BLAS, such as OpenBLAS, MKL, AOCL, ...
an implementation of ScaLAPACK;
an implementation of MPI, such as OpenMPI, MPICH, Intel® MPI, ...

The easiest setup is the Intel® oneAPI toolkit, which supplies MPI, oneMKL (with a complete ScaLAPACK) in one shot:

git clone --recurse-submodules --depth 1 https://github.com/TLCFEM/ezp.git
mkdir ezp/build && cd ezp/build
. /opt/intel/oneapi/setvars.sh && cmake -DEZP_TEST=ON .. && cmake --build . --config Release

Using system libraries is also possible — see the compilation guide for openSUSE, Fedora and other distros.

Minimal Example

It is assumed that the root node (rank 0) prepares the left-hand side $A$ and right-hand side $B$. The solver then distributes the matrices across the available processes, solves the system, and returns the solution back to the root node.

All BLACS and ScaLAPACK details are hidden: prepare the matrices (on the root node), wrap them, and call solve. The result highly resembles how one would solve a linear system sequentially.

#include <ezp/pgesv.hpp>
#include <iomanip>
#include <iostream>
 
using namespace ezp;
 
int main() {
    // get the current blacs environment
    const auto env = get_env<int>();
 
    constexpr auto N = 6, NRHS = 2;
 
    // storage for the matrices A and B
    std::vector<double> A, B;
 
    // helper functor to convert 2D indices to 1D indices
    // it's likely the matrices are provided by some other subsystem
    const auto IDX = par_dgesv<int>::indexer{N};
 
    if(0 == env.rank()) {
        // the matrices are only initialized on the root process
        A.resize(N * N, 0.);
        B.resize(N * NRHS, 1.);
 
        for(auto I = 0; I < N; ++I) A[IDX(I, I)] = static_cast<double>(I + 1);
    }
 
    // create a parallel solver
    // it takes the number of rows and columns of the process grid as arguments,
    // or let the library automatically determine them
    auto solver = par_dgesv<int>();
 
    // wrap the data in full_mat objects:
    // number of rows, number of columns, and a pointer to the data
    // on non-root processes the pointer is nullptr, as the vectors are empty
    const auto info = solver.solve({N, N, A.data()}, {N, NRHS, B.data()});
 
    // reuse the stored factorization to compute the determinant
    const auto det = solver.det({N, N, A.data()});
 
    if(0 == env.rank() && 0 == info) {
        std::cout << std::setprecision(10) << "Determinant: " << det << '\n';
        std::cout << "Solution:\n";
        for(const auto i : B) std::cout << i << '\n';
    }
 
    return info;
}

More runnable examples live in examples/, covering every solver and input format.

Standalone Solvers

MPI can dynamically spawn worker processes via MPI_Comm_spawn. ezp exploits this to ship its solvers as standalone executables, cleanly separating the main application from the solver.

Standalone solver architecture

This is particularly handy when the main application runs a different parallelism model, has latency constraints, or simply wants a clear boundary around the solver. Because the boundary is a plain MPI interface, callers can be written in any MPI-compatible language — the bundled runner.py demonstrates doing it from Python with mpi4py.

Enable them at configure time with -DEZP_STANDALONE=ON; see the standalone guide.

Documentation

Full guides, per-solver references and the API docs are published at **https://tlcfem.github.io/ezp/**:

Motivation and Compilation
Usage, Custom matrices, MPI management, Standalone solvers
Per-solver reference for dense and sparse solvers

‍[!IMPORTANT]

Use branch master for development and testing.

Use branch lite-lib for integration into other projects:
git clone --recursive --branch lite-lib https://github.com/TLCFEM/ezp.git

License

ezp is distributed under the terms of the GNU General Public License v3.0. See [LICENSE](LICENSE) for details.