Skip to content

Motivation

As of this writing, it is necessary to directly call ScaLAPACK subroutines via the C interface if one wants to solve linear systems using ScaLAPACK. It is, however, not a trivial task, compared to solving similar problems using LAPACK. One needs to figure out how to distribute the matrix over the process grid, and allocate the necessary buffers (with appropriate sizes) for the ScaLAPACK subroutines. It requires quite a bit of boilerplate code, which is not only tedious to write, but also error-prone.

Since both MPI and LAPACK/BLAS have fabulous C++ wrappers (see mpl and Armadillo), it would be nice if ScaLAPACK can be used in a similar OO manner. Such a wrapper is advantageous and beneficial in the context of C++ applications due to the following reasons.

  1. It avoids direct interactions with the C API, which requires manual memory management. The caller has to follow pretty much the same procedural style.
  2. A well-designed OO approach can significantly reduce the code size.

ezp is designed to provide an intuitive interface that resembles the non-MPI version of linear systems solvers. The majority of ScaLAPACK, BLACS and MPI communication details shall be well hidden from the caller so that solving a system can be as simple as solver.solve(A, B). Of course, a complete, working example would be longer than a single line of code, due to the distributed-memory nature of MPI. However, as shown in the front page, a minimum example spans merely a few lines of code.

C++
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
#include <ezp/pgesv.hpp>
#include <iostream>

int main() {
    constexpr auto N = 6, NRHS = 2;

    std::vector<double> A, B;

    if(0 == ezp::get_env<int>().rank()) {
        A.resize(N * N, 0.);
        B.resize(N * NRHS, 1.);

        const auto IDX = ezp::par_dgesv<int>::indexer{N};
        for(auto I = 0; I < N; ++I) A[IDX(I, I)] = I;
    }

    ezp::par_dgesv<int>().solve({N, N, A.data()}, {N, NRHS, B.data()});

    if(0 == ezp::get_env<int>().rank()) for(const auto i : B) std::cout << i << '\n';

    return 0;
}

Available Solvers

ezp implements all precisions (DSZC) with both 32-bit and 64-bit indexing.

Dense Matrix

ezp provides all solvers listed in ScaLAPCK page for full and band matrices. Solvers for tridiagonal matrices are not implemented as they can be stored in band formats and solved by the corresponding solver.

Sparse Matrix

ezp provides interface to PARDISO direct solver bundled in Intel® oneAPI Math Kernel Library (oneMKL).

ezp provides interface to MUMPS direct solver.

ezp provides interface to Lis iterative solver.

Why not SuperLU?

The codebase is messy and hard to maintain due to, for example, abuse of macros. And it is not significantly faster than other direct solvers. I do not like its code aesthetics.

References

  1. Configuration of a linear solver for linearly implicit time integration and efficient data transfer in parallel thermo-hydraulic computations
  2. A Parallel Geometric Multifrontal Solver Using Hierarchically Semiseparable Structure
  3. Caveats of three direct linear solvers for finite element analyses
  4. Evaluating Accuracy and Efficiency of HPC Solvers for Sparse Linear Systems with Applications to PDEs