PETScMatrix allocation assumes sparsity pattern divisible by block size

[Jan Blechta]

Issue #721 revealed that PETScMatrix::init assumes here that sparsity pattern is regular, divisible by block size, i.e. that the block contributes to the pattern as a whole, full block, without inner structure. On the other hand sparsity pattern implementation does not use block size information in any way and allows allocations with any granularity within the block.

The assumption is not met at least in two cases:

1. diagonal option of SparsityPatternBuilder is used, then 1-diagonal (rather than bs-diagonal) is allocated; this caused #721 and was worked around by ceiling division in 3365a1ccf50;

2. user input to SparsityPattern.insert_foo() is not divisible by blocks; then PETScMatrix::init will assume a block matching the first row and column of the block which might not be the case;

The question is whether we want to fix the flaw 2. even though SparsityPattern interface is probably very scarcely used by users.

Possible fixes:

1. either make SparsityPattern::insert_foo work block-wise, only allowing nodes matching to full blocks (involves interface change in SparsityPattern: indexing modulo blocks); advantage is the this will shrink sparsity pattern data by multiple of 1/bs,

2. or fix PETScMatrix::init to take into account all the rows and columns from every block into account.

Also don't forget to check the situation in other backends.