ezp

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

Features

1. easy to use interface
2. drop-in header-only library
3. standalone solver binaries that can be invoked by various callers

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 distrbute 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.

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