What are d_nz and o_nz in PETSc's MatMPIAIJSetPreallocation function?

Question

I'm working from ex2.c of PETSc's example codes. On line 65, the code specifies:

MatMPIAIJSetPreallocation(A,5,PETSC_NULL,5,PETSC_NULL);  

I looked up the documentation for MatMPIAIJSetPreallocation(), and understand that its purpose is to allocate memory for the matrix in a compressed sparse row format in parallel. That is, to allocate memory only for the maximum number of non-zeros expected in each row. According to this documentation, the second argument is

d_nz - number of nonzeros per row in DIAGONAL portion of local submatrix

I'm not quite sure what is meant by "diagonal portion of local submatrix". Furthermore, why are 5 values needed for both the diagonal (d_nz) and off-diagonal (o_nz) portions of each row? As I understand, the 5-point stencil laplacian problem being solved has a linear system of equations which requires a maximum of 5 non-zero entries per row. The d_nz and o_nz values specified here seem to indicate to me that memory has been allocated for 10 values per row. I'm sure I'm missing something in the interpretation of d_nz and o_nz.

Answer 1

Distributed vector spaces and matrix storage formats

PETSc vectors are distributed across the communicator with each process owning a contiguous block of rows. The most common case is for all the parts to be the same size, but that is not necessary. Any partition of the $N$ entries into $P$ contiguous parts (some possibly empty) with the parts labeled sequentially is valid. The Mat class represents a linear transformation between two finite dimensional vector spaces, $A: X \to Y$. The spaces $X$ and $Y$ may be different (even different sizes), though Krylov methods and most linear solvers require that $X = Y$. There are routines for looking at the ownership range of the current process:

• VecGetOwnershipRange(x,&start,&end) returns range of owned entries for any vector $x$.
• MatGetOwnershipRangeColumn(A,&cstart,&cend) returns the owned part of the domain $X$.
• MatGetOwnershipRange(A,&rstart,&rend) returns the owned part of the range $Y$.

The storage format used by PETSc MPI*AIJ matrices distributes the entries according to ownership of the range space $Y$. In this model, there are generally some entries of the domain $X$ that affect the owned part of $Y$. To be precise, let $X_i \subset X$ be the part of the domain space $X$ that is owned by process (MPI rank) $i$ and similarly for $Y_i \subset Y$. The "diagonal block" of the matrix owned by rank $i$ is exactly the effect of $X_i$ on $Y_i$. When performing multiplication, the part of a vector in $X_i$ is available locally without any communication. The "off-diagonal block" represents the effect of $X \setminus X_i$ on $Y_i$, which cannot be applied without communicating entries. The operation $y \gets A x$ is computed as

1. Start communicating off-process entries of $x$ to $x_i^{\text{off}}$.
2. Apply the diagonal block of the matrix $y_i \gets A_i^{\text{diag}} x_i$.
3. Finish communicating off-process entries.
4. Add the action of the off-diagonal block $y_i \gets y_i + A_i^{\text{off}} x_i^{\text{off}}$

Geometrically, the off-diagonal part corresponds to the influence of the "ghost points" on the subdomain.

Preallocation

Since the diagonal part $A_i^{\text{diag}}$ is stored separately from $A_i^{\text{off}}$, we need separate preallocation information to avoid dynamic data structures. This can be done by providing one parameter indicating the maximum number of nonzeros in any row (*_nz) or by providing an array indicating the number of nonzeros in each row (*_nnz). Note that rows corresponding to "interior" points will generally have no nonzero entries in the off-diagonal part.

As for ex2.c, there are indeed never more than 5 nonzeros per row of the matrix, of which at least one (the diagonal entry) must fall in the diagonal block. Therefore it would have been safe to pass o_nz=4. Note that this extreme case o_nz=4 is only realized by this example for a single-point subdomain.

Answer 2

If you split a matrix so that it is stored across processes in an MPI program (say, the matrix is $N\times N$, and we have $M$ processors), my understanding is that PETSc stores it in such a way that every processor stores a matrix that represents the diagonal block (size $N/M \times N/M$) and one that represents the other entries of the $N/M$ rows of the matrix that are stored on the local processor. It's done this way since every processor stores $N/M$ entries of a vector, and doing a matvec operation is simplest if you can just do a matvec on the diagonal block of the matrix with the locally stored part of the matrix, and then fix up the rest with the remaining matrix entries stored locally and the vector entries imported from the other processors.

Given this storage scheme, I believe that the diagonal and non-diagonal matrix entries you preallocate correspond to the diagonal $N/M \times N/M$ matrix and the $N/M \times N-N/M$ remainder block.