Skip to content

ezp

codecov master

ezp is a lightweight C++ wrapper for selected distributed solvers for linear systems.

You May Be Interested if You...

  1. want to solve linear systems on distributed clusters
  2. want to avoid the hassles of learning how to call various C/Fortran solvers

Solving a linear system shall be as easy as calling a method solver.solve(A, B).

Features

  1. easy to use interface
  2. drop-in header-only library
  3. standalone solver binaries that can be invoked by various callers
  4. implementation validated with randomly generated inputs

The following solvers are implemented.

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 (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

[!IMPORTANT] * Use branch master for development and testing. * Use branch lite-lib for integration into other projects.

Dependency

The ezp library requires C++ 20 compatible compiler. The following drivers are needed.

  1. an implementation of LAPACK and BLAS, such as OpenBLAS, MKL, etc.
  2. an implementation of ScaLAPACK
  3. an implementation of MPI, such as OpenMPI, MPICH, etc.

Example

It is assumed that the root node (rank 0) prepares the left hand side \(\(A\)\) and right hand side \(\(B\)\). The solvers distribute the matrices to available processes and solve the system, return the solution back to the master node.

The solvers are designed in such a way that all BLACS and ScaLAPACK details are hidden. One shall prepare the matrices (on the root node) and call the solver. The following is a typical example. It highly resembles the sequential version of how one would typically solve a linear system.

The following is a working example.

C++
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
#include <ezp/pgesv.hpp>
#include <iomanip>
#include <iostream>

using namespace ezp;

int main() {
    // get the current blacs environment
    const auto rank = get_env<int>().rank();

    constexpr auto N = 6, NRHS = 2;

    // storage for the matrices A and B
    std::vector<double> A, B;

    if(0 == rank) {
        // the matrices are only initialized on the root process
        A.resize(N * N, 0.);
        B.resize(N * NRHS, 1.);

        // 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};

        for(auto I = 0; I < N; ++I) A[IDX(I, I)] = static_cast<double>(I);
    }

    // create a parallel solver
    // it takes the number of rows and columns of the process grid as arguments
    // or let the library automatically determine as follows
    // need to wrap the data in full_mat objects
    // it requires the number of rows and columns of the matrix, and a pointer to the data
    // on non-root processes, the data pointer is nullptr as the vector is empty
    // par_dgesv<int>().solve(full_mat{N, N, A.data()}, full_mat{N, NRHS, B.data()});
    par_dgesv<int>().solve({N, N, A.data()}, {N, NRHS, B.data()});

    if(0 == rank) {
        std::cout << std::setprecision(10) << "Solution:\n";
        for(auto i = 0; i < B.size(); ++i) std::cout << B[i] << '\n';
    }

    return 0;
}