Function to restrict matrices

[peterrum]

This PR introduces functions that allow to extract parts of a distributed sparse matrix. In the context of Schwarz and Domain Decomposition methods this is known as restriction: A_i = R_i A R_i^T, where R_i is Boolean matrix. The granularity of the restriction is either on a cell level (restrict_to_cells()) or can be specified by an index set (restrict_to_serial_sparse_matrix()).

I have create a new namespace with the name SparseMatrixTools. Is there already a more suitable place for such functionality?

[peterrum]

@kronbichler Do you thinks such a function would be helpful? We have used prefix sums already at so many places! If yes, I would generalize the function so that it also works for vectors (as we need to process data for all multigrid levels simultaneously).

[kronbichler]

I think that function would make sense, yes.

[drwells]

Could we put these in MatrixTools instead? We don't use that namespace much but I would prefer not to introduce a new one that is similar.

[peterrum]

Could we put these in MatrixTools instead? We don't use that namespace much but I would prefer not to introduce a new one that is similar.

I know the MatrixTools and MatrixCreator namespace. But both are in the numerics folder (which does not fit for the current use case) and the purpose of the functions are somewhat different (modify the matrix according to some boundary condition and create a matrix for a special type of operator (mass matrix/Laplace matrix)).

[kronbichler]

Regarding the namespace, I have no strong opinion. SparseMatrixTools sounds OK in my opinion. One question I'm wondering, without having looked at it: Did you compare this to the sparse_vanka.h file? I would assume something similar to happen there, but the last time I looked in that file was in 2009.

[peterrum]

Regarding the namespace, I have no strong opinion. SparseMatrixTools sounds OK in my opinion. One question I'm wondering, without having looked at it: Did you compare this to the sparse_vanka.h file? I would assume something similar to happen there, but the last time I looked in that file was in 2009.

@kronbichler I went through the SparseVanka code. It is serial and I fail how it would help. Indeed, the code extracts submatrices from a (serial) sparse matrix:

dealii/include/deal.II/lac/sparse_vanka.templates.h

Line 199 in 05b61d7

inverses[row]->extract_submatrix_from(*matrix, local_indices, local_indices);

which is the following way:

dealii/include/deal.II/lac/full_matrix.h

Lines 1207 to 1225 in 900dc61

	template <typename number>
	template <typename MatrixType, typename index_type>
	inline void
	FullMatrix<number>::extract_submatrix_from(
	const MatrixType & matrix,
	const std::vector<index_type> &row_index_set,
	const std::vector<index_type> &column_index_set)
	{
	AssertDimension(row_index_set.size(), this->n_rows());
	AssertDimension(column_index_set.size(), this->n_cols());
	const size_type n_rows_submatrix = row_index_set.size();
	const size_type n_cols_submatrix = column_index_set.size();
	for (size_type sub_row = 0; sub_row < n_rows_submatrix; ++sub_row)
	for (size_type sub_col = 0; sub_col < n_cols_submatrix; ++sub_col)
	(*this)(sub_row, sub_col) =
	matrix.el(row_index_set[sub_row], column_index_set[sub_col]);
	}

That code, however, is not able to handle the two major issues tackled in the new functions in this PR:

dealii/include/deal.II/lac/sparse_matrix_tools.h

Lines 130 to 134 in 0579b4c

	* @note In a certain sense, this is the reversion of the cell loop during
	* matrix assembly. However, doing this on a distributed matrix is not
	* trivial, since 1) rows might be owned by different processes and 2) degrees
	* of freedoms might be constrained, resulting in "missing" entries in the
	* matrix.

[peterrum]

I have updated PR and removed the dependency to #14124. Ready from my side.

[kronbichler]

Can you explain which communicator you refer to? The one of the system_matrix?

[kronbichler]

I think this needs to be adjusted, because we also have ChunkSparseMatrix and some similar classes, that are similar to dealii::SparseMatrix.

[peterrum]

Done.

[kronbichler]

I think this needs to be documented - we support the case of empty indices, in which case the result is also empty (clearly).

[peterrum]

@kronbichler I have fixed the typos and updated the documentation.