HowTo: Use Block Operators and Matrices

Some problem formulations are defined in block form and need to be implemented in terms of block operators. Examples include saddle point problems (ex5.cpp), DPG discretization (ex8.cpp), and problems with multiple variables (ex19.cpp). The resulting discretized system is expressed in terms of block operators and vectors, which may be distributed in parallel. This article gives an overview of working with block operators and their matrix representations.

It should be noted in general that operators and matrices are appropriate in different situations, regardless of whether they are in block form. Generally, it is preferable to have an operator and not its matrix representation when only its action is needed and can be computed faster than matrix assembly, or when matrix storage requires too much memory. For example, this is the case for high-order FEM, when partial assembly (PA) is used for fast operator multiplication on GPUs without storing matrices. Also, matrix storage becomes increasingly expensive (more nonzeros per row) as FEM order increases, which is another reason to avoid matrix assembly and matrix-based preconditioners for very high order. On the other hand, for low-order FEM, matrices are necessary for example in order to use AMG preconditioning (e.g. with hypre). Thus there are cases where operators or matrices are preferable, in general and in block form.

First, it is important to understand how a single, monolithic operator or matrix is distributed in parallel in MFEM. Vectors, matrices, and operators are distributed consistently with hypre, which decomposes the rows of a parallel matrix (HypreParMatrix, see mfem/hypre.hpp) but stores all columns of the locally owned rows on each MPI rank. On each process, a Vector or HypreParVector is of size equal to the number of locally owned rows, and a HypreParMatrix stores the local rows. The parallel communication necessary for matrix-vector multiplication is performed in hypre. Similarly, an Operator should act on a Vector of local entries, perform any necessary communication, and compute a Vector of local entries.

In the case of block operators and vectors, a Vector stores the local entries for each block contiguously in its data. Offsets define where each block begins and ends. For example, in ex5.cpp, there are two blocks for spaces R_space and W_space, and block_offsets is of size three, storing offsets 0, R_space->GetVSize(), and R_space->GetVSize() + W_space->GetVSize(). The class BlockOperator (see mfem/linalg/blockoperator.hpp) can be used to form one operator from operators defining the blocks. It operates on vectors of local entries, stored block-wise.

Similarly, a monolithic HypreParMatrix can be constructed, using the function HypreParMatrixFromBlocks (see hypre.hpp), from blocks defined as HypreParMatrix pointers or null pointers for empty blocks. The blocks may be rectangular, but their sizes must be consistent. Scalar coefficients can optionally be used. The monolithic matrix will have copies of the entries from the blocks, so it can be modified or destroyed independently of the blocks.

The unit test mfem/tests/unit/linalg/test_matrix_rectangular.cpp provides an example that compares a BlockOperator and a monolithic HypreParMatrix. As noted above, it is not practical to have both an operator and a matrix, but this test illustrates the equivalence of the two approaches.

The capability to form a monolithic matrix is available only for HypreParMatrix, not for the serial class SparseMatrix.