ezp is a lightweight C++ wrapper for selected ScaLAPACK 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 |
|---|
| ✔️ | general (partial pivoting) | simple | PxGESV |
| ✔️ | general (partial pivoting) | expert | PxGESVX |
| ✔️ | symmetric/Hermitian positive definite | simple | PxPOSV |
| ✔️ | symmetric/Hermitian positive definite | expert | PxPOSVX |
| ✔️ | general band (partial pivoting) | simple | PxGBSV |
| ✔️ | general band (no pivoting) | simple | PxDBSV |
| ✔️ | symmetric/Hermitian positive definite band | simple | PxPBSV |
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 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 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;}