#include <mpi.h>
#include <hpc/mpi/matrix.h>
#include <hpc/matvec/gematrix.h>
#include <hpc/matvec/print.h>
#include <hpc/matvec/apply.h>
int main(int argc, char** argv) {
MPI_Init(&argc, &argv);
int nof_processes; MPI_Comm_size(MPI_COMM_WORLD, &nof_processes);
int rank; MPI_Comm_rank(MPI_COMM_WORLD, &rank);
using namespace hpc::matvec;
using namespace hpc::mpi;
using Matrix = GeMatrix<double>;
using Index = typename Matrix::Index;
int share = 3;
int num_rows = nof_processes * share;
int num_cols = 5;
Matrix B(share, num_cols, StorageOrder::RowMajor); /* individual share */
MPI_Datatype rowtype_B = get_row_type(B);
if (rank == 0) {
Matrix A(num_rows, num_cols); /* entire matrix */
apply(A, [](double& val, Index i, Index j) -> void {
val = i * 100 + j;
});
MPI_Datatype rowtype_A = get_row_type(A);
MPI_Scatter(&A(0, 0), share, rowtype_A,
&B(0, 0), share, rowtype_B,
0, MPI_COMM_WORLD);
apply(B, [=](double& val, Index i, Index j) -> void {
val += 10000 * (rank + 1);
});
MPI_Gather(&B(0, 0), share, rowtype_B,
&A(0, 0), share, rowtype_A,
0, MPI_COMM_WORLD);
print(A, "A");
} else {
MPI_Scatter(nullptr, 0, nullptr, /* ignored parameters */
&B(0, 0), share, rowtype_B,
0, MPI_COMM_WORLD);
apply(B, [=](double& val, Index i, Index j) -> void {
val += 10000 * (rank + 1);
});
MPI_Gather(&B(0, 0), share, rowtype_B,
nullptr, 0, nullptr, /* ignored parameters */
0, MPI_COMM_WORLD);
}
MPI_Finalize();
}