Standalone Solver

MPI provides a mechanism to dynamically spawn worker processes through MPI_Comm_spawn. By exploiting this feature, it is possible to separate the main application from the solver. The architecture is illustrated below.

architecture design

The ezp solvers are also available as standalone executables that follow this design. Such a form is particularly beneficial when one or more of the following holds.

The main application is not suitable for a distributed-memory model, for example due to latency constraints.
The main application adopts a different parallelism pattern.
One wants a clear boundary between the main application and the solver.

The clear boundary can also be a shortcoming. The root process must store the whole matrices in memory, which is subject to physical limitations. Because the spawned processes are not persistent, solving the same system for different right-hand sides may incur unnecessary data communication. If the matrices are read from file in a distributed manner, one may also need fine control over how the data is prepared, possibly by calling BLACS/ScaLAPACK/MPI functions directly.

Backbone

The pgesv solver is used below to illustrate the critical steps of implementing this architecture. ezp does not use the raw MPI functions directly; instead, the mpl library is used. It is necessary to obtain the default communicator and the parent inter-communicator.

solver.full.hpp:24:25
inline const auto& comm_world{mpl::environment::comm_world()};
inline const auto& parent = mpl::inter_communicator::parent();

Because mpl handles MPI environment initialisation and finalisation, ezp must be told not to finalise MPI. The function ezp::blacs_env<>::do_not_manage_mpi() must be called when MPI is managed externally.

solver.pgesv.cpp:54:55
    ezp::blacs_env<>::do_not_manage_mpi();

The parameters of the linear system are broadcast from the main application first.

solver.pgesv.cpp:61:65
    const auto all = mpl::communicator(parent, mpl::communicator::order_high);

    int config[3]{};

    all.bcast(0, config);

Knowing the problem sizes N and NRHS, the root process initialises the containers and receives the contents of the two matrices.

solver.full.hpp:28:40
    std::vector<DT> A, B;

    if(0 == comm_world.rank()) {
        A.resize(N * N);
        B.resize(N * NRHS);

        mpl::irequest_pool requests;

        requests.push(parent.irecv(A, 0, mpl::tag_t{0}));
        requests.push(parent.irecv(B, 0, mpl::tag_t{1}));

        requests.waitall();
    }

Solving the system follows the conventional approach: create a solver object, wrap the data, and call solve. The solution is then sent back to the caller.

solver.full.hpp:42:47
    const auto error = solver_t().solve({N, N, A.data()}, {N, NRHS, B.data()});

    if(0 == comm_world.rank()) {
        parent.send(error, 0);
        if(0 == error) parent.send(B, 0);
    }

The Caller

From the code snippets above, the caller must

broadcast config;
send A and B;
receive the error code and, if no error occurred, the solution X.

To this end, one declares the corresponding containers with std::vector.

runner.cpp:41:42
    std::vector<double> A, B(N * NRHS, 1.);
    std::vector<int> config;

The following creates a diagonal matrix for illustration.

runner.cpp:48:53
        solver = "solver.pgesv";

        config = {N, NRHS, 1};

        A.resize(N * N, 0.);
        for(auto I = 0; I < N; ++I) A[I * N + I] = I + 1.;

Once all data is ready to be communicated, the communication itself is concise and straightforward.

runner.cpp:114:129
    const auto& comm_world{mpl::environment::comm_world()};
    const auto worker = comm_world.spawn(0, std::abs(std::stoi(argv[2])), {solver});
    const auto all = mpl::communicator(worker, mpl::communicator::order_low);

    all.bcast(0, config.data(), mpl::contiguous_layout<int>(config.size()));

    mpl::irequest_pool requests;

    requests.push(worker.isend(A, 0, mpl::tag_t{0}));
    requests.push(worker.isend(B, 0, mpl::tag_t{1}));

    requests.waitall();

    int error = -1;
    worker.recv(error, 0);
    if(0 == error) worker.recv(B, 0);

Python Caller

Because there is a clear boundary between the main application and the solver, any MPI-compatible caller can exchange the problem data. For example, one can use mpi4py to call the solvers.

Note that the column-major memory layout must be preserved. NumPy arrays are row-major by default, so they must be created with order="F".

runner.py:41:76
def run(nprocs: int, N: int, NRHS: int):
    comm = MPI.COMM_WORLD

    # spawn the solver
    worker = comm.Spawn("./solver.pgesv", maxprocs=nprocs)
    all = worker.Merge()

    # broadcast the problem configuration
    all.Bcast(array("i", [N, NRHS, 1]), root=0)

    # send the matrices
    # ! numpy arrays are row-major by default
    # ! column-major order is required by (Sca)LAPACK
    # ! thus need to specify the order as "F"
    A = numpy.zeros((N, N), dtype=numpy.float64, order="F")
    for i in range(N):
        A[i, i] = i + 1
    print("A:\n", A)
    worker.Send(A, dest=0, tag=0)

    B = numpy.ones((N, NRHS), dtype=numpy.float64, order="F")
    for i in range(NRHS):
        B[:, i] = i + 1
    print("B:\n", B)
    worker.Send(B, dest=0, tag=1)

    # receive the error code and the solution if no error
    error = array("i", [-1])
    worker.Recv(error)
    if 0 == error[0]:
        X = numpy.ones((N, NRHS), dtype=numpy.float64, order="F")
        worker.Recv(X)
        print("Solution:\n", X)

    worker.Disconnect()
    all.Disconnect()

Calling the script with

./runner.py 2 6 2

produces the following correct solution.

A:
 [[1. 0. 0. 0. 0. 0.]
 [0. 2. 0. 0. 0. 0.]
 [0. 0. 3. 0. 0. 0.]
 [0. 0. 0. 4. 0. 0.]
 [0. 0. 0. 0. 5. 0.]
 [0. 0. 0. 0. 0. 6.]]
B:
 [[1. 2.]
 [1. 2.]
 [1. 2.]
 [1. 2.]
 [1. 2.]
 [1. 2.]]
Solution:
 [[1.         2.        ]
 [0.5        1.        ]
 [0.33333333 0.66666667]
 [0.25       0.5       ]
 [0.2        0.4       ]
 [0.16666667 0.33333333]]

Full Reference Implementation

The following is a full reference implementation of a standalone solver and its caller.

runner.cpp

runner.cpp
#include <iostream>
#include <mpl/mpl.hpp>

int main(int argc, char** argv) {
    if(argc < 4) {
        std::cout << "Usage: runner ge|gex|po|pox|gb|db|pb [n] [m]\n\n";
        std::cout << "    Generate a matrix of size [m] and solve it with [n] processes using ge|gex|po|pox|gb|db|pb solver.\n";
        std::cout << "Example: runner ge 3 100\n";
        return 0;
    }
    constexpr auto NRHS = 1;

    const auto N = std::stoi(argv[3]);

    std::vector<double> A, B(N * NRHS, 1.);
    std::vector<int> config;

    const auto type = std::string(argv[1]);
    std::string solver;

    if("ge" == type) {
        solver = "solver.pgesv";

        config = {N, NRHS, 1};

        A.resize(N * N, 0.);
        for(auto I = 0; I < N; ++I) A[I * N + I] = I + 1.;
    }
    else if("gex" == type) {
        solver = "solver.pgesvx";

        config = {N, NRHS, 1};

        A.resize(N * N, 0.);
        for(auto I = 0; I < N; ++I) A[I * N + I] = I + 1.;
    }
    else if("po" == type) {
        solver = "solver.pposv";

        config = {N, NRHS, 1};

        A.resize(N * N, 0.);
        for(auto I = 0; I < N; ++I) A[I * N + I] = I + 1.;
    }
    else if("pox" == type) {
        solver = "solver.pposvx";

        config = {N, NRHS, 1};

        A.resize(N * N, 0.);
        for(auto I = 0; I < N; ++I) A[I * N + I] = I + 1.;
    }
    else if("gb" == type) {
        solver = "solver.pgbsv";

        constexpr auto KL = 1, KU = 1;

        config = {N, KL, KU, NRHS, 1};

        A.resize(N * (KL + KU + 1), 0.);
        for(auto I = 0; I < N; ++I) A[KU + I * (KL + KU + 1)] = I + 1.;
    }
    else if("db" == type) {
        solver = "solver.pdbsv";

        constexpr auto KL = 1, KU = 1;

        config = {N, KL, KU, NRHS, 1};

        A.resize(N * (KL + KU + 1), 0.);
        for(auto I = 0; I < N; ++I) A[KU + I * (KL + KU + 1)] = I + 1.;
    }
    else if("pb" == type) {
        solver = "solver.ppbsv";

        constexpr auto KLU = 1;

        config = {N, KLU, NRHS, 1};

        A.resize(N * (KLU + 1), 0.);
        for(auto I = 0; I < N; ++I) A[I + I * KLU] = I + 1.;
    }
    else {
        std::cout << "Usage: runner ge|gex|po|pox|gb|db|pb [n] [m]\n";
        return 0;
    }

    const auto& comm_world{mpl::environment::comm_world()};
    const auto worker = comm_world.spawn(0, std::abs(std::stoi(argv[2])), {solver});
    const auto all = mpl::communicator(worker, mpl::communicator::order_low);

    all.bcast(0, config.data(), mpl::contiguous_layout<int>(config.size()));

    mpl::irequest_pool requests;

    requests.push(worker.isend(A, 0, mpl::tag_t{0}));
    requests.push(worker.isend(B, 0, mpl::tag_t{1}));

    requests.waitall();

    int error = -1;
    worker.recv(error, 0);
    if(0 == error) worker.recv(B, 0);

    for(auto I = 0; I < N; ++I) printf("x[%d] = %+.6f\n", I, B[I]);

    return 0;
}

solver.pgesv.cpp

solver.pgesv.cpp
#include "solver.full.hpp"

int main(int, char**) {
    ezp::blacs_env<>::do_not_manage_mpi();

    if(!parent.is_valid()) {
        printf("This program must be invoked by the host application.\n");
        return 0;
    }

    const auto all = mpl::communicator(parent, mpl::communicator::order_high);

    int config[3]{};

    all.bcast(0, config);

    const auto N = config[0];
    const auto NRHS = config[1];
    const auto FLOAT = config[2];

    if(FLOAT >= 10) return run<complex16, int_t, ezp::pgesv<complex16, int_t>>(N, NRHS);
    if(FLOAT >= 0) return run<double, int_t, ezp::pgesv<double, int_t>>(N, NRHS);
    if(FLOAT > -10) return run<float, int_t, ezp::pgesv<float, int_t>>(N, NRHS);

    return run<complex8, int_t, ezp::pgesv<complex8, int_t>>(N, NRHS);
}