Topics about the deal.II finite element library class "SparsityPattern"

Question

When learning the deal.II FE library, I am a bit confused about the mechanism of its "SparsityPattern" class. Through reading the documentation, I only got to know that it uses the Compressed Row Storage format to store indices of nonzero entries of a sparse matrix.

To put my confusion be explicit, suppose I have a quare 10x10 sparse matrix which stores the values corresponding to 10 degrees of freedom, namely dof_handler.n_dofs = 10:

\begin{pmatrix} 1 & 2 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\ 3 & 4 & 5 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 6 & 7 & 8 & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 9 & 10 & 11 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 12 & 13 & 14 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 15 & 16 & 17 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & 18 & 19 & 20 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & 0 & 21 & 22 & 23 & 0\\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 24 & 25 & 26\\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 27 & 28\\ \end{pmatrix}

But at first, suppose we don't know exactly where the nonzero entries are and can only get a maximal estimate of their number at each row, say max_per_row, and let max_per_row = 5. Then the sparsity pattern would be:

\begin{pmatrix} X & X & 0 & 0 & 0 & 0 & 0 & X & X & X\\ X & X & X & 0 & 0 & 0 & 0 & 0 & X & X\\ X & X & X & X & X & 0 & 0 & 0 & 0 & 0\\ 0 & X & X & X & X & X & 0 & 0 & 0 & 0\\ 0 & 0 & X & X & X & X & X & 0 & 0 & 0\\ 0 & 0 & 0 & X & X & X & X & X & 0 & 0\\ 0 & 0 & 0 & 0 & X & X & X & X & X & 0\\ 0 & 0 & 0 & 0 & 0 & X & X & X & X & X\\ X & 0 & 0 & 0 & 0 & 0 & X & X & X & X\\ X & X & 0 & 0 & 0 & 0 & 0 & X & X & X\\ \end{pmatrix}

One way to get the sparsitypattern is:

SparsityPattern  sparsity_pattern;
sparsity_pattern.reinit(10, 10, 5);
DoFTools::make_sparsity_pattern(dof_handler, sparsity_pattern);
sparsity_pattern.compress()

where dof_handler stores the information of degrees of freedom.

And it can also be implemented as below:

DynamicSparsityPattern dynamic_pattern (dof_handler.n_dofs());
DoFTools::make_sparsity_pattern (dof_handler, dynamic_pattern);
constraints.condense (dynamic_pattern);
SparsityPattern sp;
sp.copy_from (dynamic_pattern);

So, my question is:

*1) Here, when we call sparsity_pattern.reinit(10, 10, 5), does it create an empty 10x10 matrix (2D vector) or just a 1x100 vector storing the values row by row? And if it is 10x10 or 1x100, then what is the functionality of third parameter max_per_row=5?

*2) I think sparsity_pattern.reinit(10, 10, 5) creates a 10x10 matrix and then using DoFTools::make_sparsity_pattern(dof_handler, sparsity_pattern) to create a CRS format, i.e, two vectors storing the row columm index of nonzero entries at each row and the index starting a new row respectively, according to the dof information stored in dof_handler. Does this understanding make sense?

*3) I know in C++ STL, we can compress a vector to its actuall capacity using its member function vector.shrink_to_fit. But here if DoFTools::make_sparsity_pattern has already make a compressed format of 5 entries at each row (according to *2), how can it be compressed again using sparsity_pattern.compress()?

*4) Does DynamicSparsityPattern dynamic_pattern (dof_handler.n_dofs()) creates a (n_dofs) x (n_dofs) matrix, then what's the difference between DoFTools::make_sparsity_pattern(dof_handler, sparsity_pattern) and DoFTools::make_sparsity_pattern (dof_handler, dynamic_pattern).

*5) What about constraints.condense (dynamic_pattern)?

*6) In C++ STL, we know that we can copy a vector v1 of size 4, capacity 8 to a vector v2, then v2 will be of capacity 4. In this way, can can also compress a vector to its actuall capacity. Is it the same mechanism to use `sp.copy_from (dynamic_pattern)'?

I know these questions are so dummy to experienced users of deal.II. But to green learners as me, they are really great hurdles to jump. I sincerely hope someone could be so kind to give some help, and any comments would be greatly appreciated. Thanks in advance!

Answer 1

I will try answer based on my experience with deal.ii.

1. The max_per_row=5 means that at most there will be 5 non-zeros per row in the matrix. Since we now know this then we do not need to have a $1\times{}100$ matrix but rather a $1\times{}50$. In other words this parameter sets an upper-bound on the memory needed. In reality it is not stored as one vector but rather as two: a row pointer vector and a column indicies vector in accordance with CRS format. See 2).
2. Yes the sparsity pattern holds the row pointers and the column indices as two vectors.
3. We must potentially compress again because while we set maximum of 5 entries per row, there may in fact be less for any given row. One row might have 5 while a different row might only have 3. This results in extra entries in the column index vector that should be removed by compression. Note that these extra entries in the column index vector are often written as $-1$'s since a column index can never be negative and thus if it is $-1$ we know that no non-zero has been assigned for that column index.
4. As the documentation says the DynamicSparsityPattern is used to find the sparsity pattern while being compressed at all times. This reduces memory overhead at the expense of cpu time.

This class acts as an intermediate form of the SparsityPattern class. From the interface it mostly represents a SparsityPattern object that is kept compressed at all times. However, since the final sparsity pattern is not known while constructing it, keeping the pattern compressed at all times can only be achieved at the expense of either increased memory or run time consumption upon use. The main purpose of this class is to avoid some memory bottlenecks, so we chose to implement it memory conservative. The chosen data format is too unsuited to be used for actual matrices, though. It is therefore necessary to first copy the data of this object over to an object of type SparsityPattern before using it in actual matrices.

Another viewpoint is that this class does not need up front allocation of a certain amount of memory, but grows as necessary. An extensive description of sparsity patterns can be found in the documentation of the Sparsity patterns module.

5. The documentation has this to say:

Condense a sparsity pattern. The name of the function mimics the name of the function we use to condense linear systems, but it is a bit of a misnomer for the current context. This is because in the context of linear systems, we eliminate certain rows and columns of the linear system, i.e., we "reduce" or "condense" the linear system. On the other hand, in the current context, the functions does not remove nonzero entries from the sparsity pattern. Rather, it adds those nonzero entry locations to the sparsity pattern that will later be needed for the process of condensation of constrained degrees of freedom from a linear system.

Since this function adds new nonzero entries to the sparsity pattern, the given sparsity pattern must not be compressed. The constraint matrix (i.e., the current object) must be closed. The sparsity pattern is compressed at the end of the function.

6. We copy the DynamicSparsityPattern to a regular SparsityPattern using the copy call because apparently the dynamic one isn't in a nice enough format for later operations. How this is done internally may involve vector copies as you suggest but I am not sure.

Concrete Example

Lets run through your first SparsityPattern code snippet (keeping in mind that the exact details of how this is done in deal.ii might differ, but the basic idea I think is correct.)

SparsityPattern  sparsity_pattern;
sparsity_pattern.reinit(10, 10, 5);
DoFTools::make_sparsity_pattern(dof_handler, sparsity_pattern);
sparsity_pattern.compress()

1. Call SparsityPattern sparsity_pattern;. This just creates our sparsity pattern object with default constructor. The sparsity_pattern object contains member variables like row and col both of which are vectors (currently of size 0). Note: row and col are names I made up. They are called something else in deal.ii.
2. Call sparsity_pattern.reinit(10, 10, 5);. Our sparsity pattern initializes the row pointer and column index vectors. This might look something like: row.resize(10+1,0); and col.resize(10*5,-1);
3. Call make_sparsity_pattern. This determines where the non-zeros in our sparse matrix will be, however some of the col entries will still be -1's - i.e. cases where there was less then 5 non-zeros in the row. Otherwise row is correctly filled and col contains both -1's and other positive column indices indicating the columns where non-zeros exist.
4. Call compress. This compresses the col vector by removing any -1's.