ezp is a lightweight C++ wrapper for selected distributed solvers for linear systems.
- easy to use interface
- drop-in header-only library
- standalone solver binaries that can be invoked by various callers
- random tested implementation
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 |
Theezp library requires C++ 20 compatible compiler.The following drivers are needed.
- an implementation of
LAPACK andBLAS, such asOpenBLAS,MKL, etc. - an implementation of
ScaLAPACK - an implementation of
MPI, such asOpenMPI,MPICH, etc.
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 masternode.
The solvers are designed in such a way that allBLACS andScaLAPACK 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>usingnamespaceezp;intmain() {// get the current blacs environmentconstauto rank = get_env<int>().rank();constexprauto 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 subsystemconstauto 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'; }return0;}