How to manually create sparsity pattern for PETSc sparsity matrix in parallel

[Pai Liu]

Hi all,

I need to manually create a PETSc sparsity matrix (not the system matrix) in parallel computing for other purposes. I only know the number of nonzero values per row before I initialize the sparsity matrix. Thus I use the following codes to initialize the PETSc sparsity matrix and add entries to this matrix:

SparsityPattern sparsity_pattern;
sparsity_pattern.reinit(dof_handler_DG.n_dofs(), dof_handler_DG.n_dofs(), n_nonzero_per_row);

sparsity_pattern.compress();
filter_matrix.reinit(    locally_owned_dofs_DG,
                                 locally_owned_dofs_DG,
                                 sparsity_pattern,
                                 mpi_communicator);

...

filter_matrix.add(current_dof_index, col_indices, weights);

Then I got the following errors:

[0]PETSC ERROR: Argument out of range
[0]PETSC ERROR: Inserting a new nonzero at global row/column (0, 1) into matrix

I guess it may be caused by the improperly settings of the sparsity pattern. How can I correct these errors, and should I use dynamic sparsity pattern or just the SparsityPattern? Any help will be appreciated.

[Wolfgang Bangerth]

On 4/2/19 9:28 AM, Pai Liu wrote:
>
> I need to manually create a PETSc sparsity matrix (not the system
> matrix) in parallel computing for other purposes.

I would start by asking where the nonzero entries are going to lie. What
kind of operator do you plan on implementing for which you don't know
where the nonzeros are?

Best
W.

--
------------------------------------------------------------------------
Wolfgang Bangerth email: bang...@colostate.edu
www: http://www.math.colostate.edu/~bangerth/

[Pai Liu]

Hi Wolfgang,

 

I would start by asking where the nonzero entries are going to lie. What
kind of operator do you plan on implementing for which you don't know
where the nonzeros are?

I want to multiply this matrix with a PETSc vector (which is also stored in parallel).

When I initialize the PETSc sparsity matrix, I don't know where the nonzeros are.

Then I loop over the locally owned cells on each processor, and guadually know for each row of the matrix the locations and values of the nonzeros. Do you mean that I can't directly use the PETSc_sparsity_matrix.add() function to insert nonzero values into the matrix without telling the sparsity pattern the nonzeros locations?

If I have to tell the sparsity pattern where the nonzeros are manually, what are the necessary procedures and functions I should follow and use (since the matrix is not a system matrix, I can't use the DoFTools::make_sparsity_pattern function with an dofhandler as input)?

For example, for row 0, how can I tell the sparsity pattern the elements at (0,3) and (0,5) arenonzeros? I search the documentation, and find maybe I should use dynamic_spasity_pattern? I am so confused.

[Juan Carlos Araujo Cabarcas]

Hey Pai Liu,

At some point in my research I had to do something similar to what you want to do. In the application of the so-called Dirichlet-to-Neumann maps for Helmholtz problems, one has to assemble a dense block corresponding to the boundary degrees of freedom, so I had to include extra entries in my FE matrices. This is what I did (omitting some non-related details):

PETScWrappers::SparseMatrix                       Q;

std::vector<unsigned int> map_dofs_boundary;
DoFTools::map_dof_to_boundary_indices (    dof_handler,
                                                                                map_dofs_boundary    );

std::vector<unsigned int> dofs_boundary_lenghts ( map_dofs_boundary.size() );

Q.reinit (par.dofs, par.dofs, dofs_boundary_lenghts);
MatSetOption(Q, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);

FE assembly routine ...

in a loop i, j ...
    Q.set ( boundary_idx[i], boundary_idx[j], my_value (i,j) );

Q.compress (VectorOperation::add);

and that works for me.

I hope this may be of help to you.

Juan Carlos Araújo Cabarcas

Doctor candidate in Computational Science with emphasis in mathematics,

Umeå Universitet

[Wolfgang Bangerth]

In addition to the excellent idea by Juan Carlos:

> When I initialize the PETSc sparsity matrix, I don't know where the nonzeros are.
> Then I loop over the locally owned cells on each processor, and guadually know
> for each row of the matrix the locations and values of the nonzeros. Do you
> mean that I can't directly use the PETSc_sparsity_matrix.add() function to
> insert nonzero values into the matrix without telling the sparsity pattern the
> nonzeros locations?

The way we set up these data structures is that we want to know up front where
the nonzeros are. It's far more efficient to do it that way, but it is
possible to add nonzeros on the fly as shown in the post by Juan Carlos.

How do we know where the nonzeros are? We "simulate" what matrix assembly
would look like, i.e., in a first stage, we record those locations that we
want to write into later on, and this gives us the sparsity pattern (e.g.,
using a DynamicSparsityPattern). Then, we can in a second stage do the actual
assembly into a matrix where we already know where the nonzeros are going to be.

[Pai Liu]

Hi Juan Carlos Araújo Cabarcas,

Thank you for your code, and it is really an excellent way. I also tried to use the dynamic sparsity pattern in dealii, and use the add_entries() function to indicate the nonzeros' locations, then I successfully insert the values in the sparse matrix as you did in your code with the set() function.

[Pai Liu]

Hi Wolfgang,

Thank you so much for your kind help. I tried the dynamic sparsity pattern, and with the sparsity_pattern.add_entries() and sparse_matrix.set() function, I sucesefully insert values in the sparse matrix.