Hei,
I am currently struggling a bit with my understanding of local and
global rows/columns in distributed dense matrices in PETSc.
The attached short program is creating a local matrix in every thread
with [matrix_size_rows, matrix_size_cols] as global size. Afterwards, I
flatten the matrix into a vector, and split this vector using the
boost-function split() when running using several MPI threads. In
addition I create a dense MPI-distributed matrix based on PETSc, using
the same input sizes, and set all values to zero (after I'm not
interested in the content itself, just the shape). Finally, I retrieve
the global and local values for the number of rows and columns using
MatGetLocalSize, and the number of the first and last row in each matrix
using MatGetOwnershipRange() and print them. According to the
documentation I would expect that the PETSc-based matrix is split up
Now, for an input size of [8, 8] and a single thread I get
/Rank 0 has a total size of 8, 8
Rank 0 has an initial size of 64
Rank 0 has a local matrix size of 8, 8
Rank 0 has a global matrix size of 8, 8
Rank 0 has a matrix spanning from row 0 to 8
Rank 0 has a vector size of 64///
which is what I expect. For two threads, though, I get
/Rank 0 has a total size of 8, 8//
//Rank 0 has an initial size of 64//
//Rank 0 has a local matrix size of 4, 4//
//Rank 0 has a global matrix size of 8, 8//
//Rank 0 has a matrix spanning from row 0 to 4//
//Rank 0 has a vector size of 32//
//Rank 1 has a total size of 8, 8//
//Rank 1 has an initial size of 64//
//Rank 1 has a local matrix size of 4, 4//
//Rank 1 has a global matrix size of 8, 8//
//Rank 1 has a matrix spanning from row 4 to 8//
Rank 1 has a vector size of 32/
Here, most entries make sense, except the size of the local matrices.
Why do I get a size of [4, 4], and not a size of [4, 8]? Each row should
be contiguous in the local process, and therefore each row should
contain all columns, not only a part of it.
Is there a misunderstanding about how to use MatGetLocalSize(), or
something else?
Thanks!
Regards,
Roland Richter
/
/
#include <iostream>
#include <armadillo>
#include <petscts.h>
#include <petscmat.h>
#include <petscvec.h>
#include <fftw3.h>
#include <fftw3-mpi.h>
#include <metis.h>
#include <chrono>
#include <boost/numeric/odeint.hpp>
void test_matrix_distribution(const size_t matrix_size_rows, const size_t matrix_size_cols) {
using namespace boost::numeric::odeint;
typedef std::complex<double> data_type;
typedef std::vector<data_type> value_type;
typedef value_type state_type;
typedef mpi_state<value_type> mpi_state_type;
PetscMPIInt rank, size;
MPI_Comm_size(PETSC_COMM_WORLD, &size);
MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
arma::cx_mat initial_mat = arma::ones<arma::cx_mat>(matrix_size_rows, matrix_size_cols);
arma::cx_colvec initial_condition_vec = arma::vectorise(initial_mat);
boost::mpi::communicator boost_comm(PETSC_COMM_WORLD, boost::mpi::comm_attach);
mpi_state_type x_split(boost_comm);
value_type x_initial = arma::conv_to<value_type>::from(initial_condition_vec);
Mat input_mat;
MatCreateDense(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, matrix_size_rows, matrix_size_cols, NULL, &input_mat);
MatZeroEntries(input_mat);
MatAssemblyBegin(input_mat, MAT_FINAL_ASSEMBLY);
MatAssemblyEnd(input_mat, MAT_FINAL_ASSEMBLY);
PetscInt first_row, last_row;
MatGetOwnershipRange(input_mat, &first_row, &last_row);
PetscInt local_cols, global_cols, local_rows, global_rows;
MatGetLocalSize(input_mat, &local_rows, &local_cols);
MatGetSize(input_mat, &global_rows, &global_cols);
split(x_initial, x_split);
std::cout << "Rank " << rank << " has a total size of " << matrix_size_rows << ", " << matrix_size_cols << '\n';
std::cout << "Rank " << rank << " has an initial size of " << initial_condition_vec.n_elem << '\n';
std::cout << "Rank " << rank << " has a local matrix size of " << local_rows << ", " << local_cols << '\n';
std::cout << "Rank " << rank << " has a global matrix size of " << global_rows << ", " << global_cols << '\n';
std::cout << "Rank " << rank << " has a matrix spanning from row " << first_row << " to " << last_row << '\n';
std::cout << "Rank " << rank << " has a vector size of " << x_split().size() << '\n';
unsplit(x_split, x_initial);
MPI_Bcast(x_initial.data(), 2 * x_initial.size(), MPI_DOUBLE, 0, PETSC_COMM_WORLD);
MatDestroy(&input_mat);
}
int main(int argc, char **args) {
PetscMPIInt rank, size;
PetscInitialize(&argc, &args, (char*) 0, NULL);
MPI_Comm_size(PETSC_COMM_WORLD, &size);
MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
test_matrix_distribution (8, 8);
PetscFinalize();
return 0;
}
