How to efficiently fill in, in parallel, a PETSc matrix from a COO sparse matrix?

Considering the following COO sparse matrix format, with repeated indices:

std::vector<size_t> rows{0, 0, 1, 2, 3, 4};
std::vector<size_t> cols{0, 0, 1, 2, 3, 4};
std::vector<double> values{2, -1, 2, 3, 4, 5};


that represents the matrix

$$$$\begin{bmatrix} 1 & 0 & 0 & 0 & 0 \\ 0 & 2 & 0 & 0 & 0 \\ 0 & 0 & 3 & 0 & 0 \\ 0 & 0 & 0 & 4 & 0 \\ 0 & 0 & 0 & 0 & 5 \end{bmatrix}$$$$

So far, the code that I have is:

// fill_in_matrix.cc

static char help[] = "Fill in a parallel COO format sparse matrix.";

#include <petsc.h>
#include <vector>

int main(int argc, char **args)
{
Mat A;
PetscInt m = 5, i, Istart, Iend;

PetscCall(PetscInitialize(&argc, &args, NULL, help));

PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, m, m));
PetscCall(MatSetFromOptions(A));
PetscCall(MatSetUp(A));
PetscCall(MatGetOwnershipRange(A, &Istart, &Iend));

std::vector<PetscInt> II{0, 0, 1, 2, 3, 4};
std::vector<PetscInt> JJ{0, 0, 1, 2, 3, 4};
std::vector<PetscScalar> XX{2, -1, 2, 3, 4, 5};

for (i = Istart; i < Iend; i++)
PetscCall(MatSetValues(A, 1, &II.at(i), 1, &JJ.at(i), &XX.at(i), ADD_VALUES));

PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
PetscCall(MatView(A, PETSC_VIEWER_STDERR_WORLD));

PetscCall(MatDestroy(&A));
PetscCall(PetscFinalize());
return 0;
}


When running it with

petscmpiexec -n 4 ./fill_in_matrix


I get

Mat Object: 4 MPI processes
type: mpiaij
row 0: (0, 1.)
row 1: (1, 2.)
row 2: (2, 3.)
row 3: (3, 4.)
row 4:


Which is missing the entry of the last row.

What am I missing? Even better, which would be the best way to fill in this matrix?

• Perhaps better asked on the Petsc mailing list?
– NNN
Commented Jun 20, 2023 at 4:55

The problem is your loop iterates through all possible rows, but your COO data has more than that amount of data (only 5 rows, but you have 6 COO entries because you sum into position 0,0 twice).

There are a few ways to resolve this:

1. split your COO data so each rank adds a (roughly) equal amount.
2. split your COO data so each rank adds the data for all the rows they own.

If you're just given a COO data vector with repeats (even if it's not sorted), method 1 is easier to implement. Usually when generating matrix data though it's easier to use method 2.

This looks something like this:

int rank, nprocs;
MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
MPI_Comm_size(PETSC_COMM_WORLD, &nprocs);
for(i = rank * II.size() / nprocs; i < (rank+1) * II.size() / nprocs; ++i)
{
PetscCall(MatSetValues(A, 1, &II.at(i), 1, &JJ.at(i), &XX.at(i), ADD_VALUES));
}


edit:

PETSc's API for inserting/adding values to matrices and vectors using their API allows you to write to any entry from any rank; it doesn't have to be part of the ownership range of that rank.

So all we need to do is divide up the COO data such that every entry is added exactly once, no more and no less. It doesn't matter which rank adds an entry; you could add them all from rank 0, or split them up between different ranks to add.

My example code above splits the COO data up roughly evenly (as evenly as possible) between the ranks. When running with 4 ranks, you have:

rank 0 inserts entry 0 rank 1 inserts entries 1 and 2 rank 2 inserts entry 3 rank 3 inserts entries 4 and 5

It finds these non-overlapping bounds by taking advantage of integer division automatically rounding/truncating in a consistent fashion. You can try working out by hand how this partitioning strategy works with a few example lengths/number of ranks.

PETSc will do all the heavy lifting to make sure that the correct rank gets the data regardless of what original rank made the API call.

• Please, could you explain why your code works? I mean, how PETSc "knows" which rank is running and distribute the data to it? Can I use this approach to fill in any matrix/vector? Thank you for answering. Commented Jun 20, 2023 at 19:16
• What do I need to change if I want to use only one rank to fill in the matrix? Commented Jun 21, 2023 at 14:26