[petsc-users] Algorithms for counting nonzeros!

[Satish Balay]

On Wed, 13 Nov 2013, Eric Chamberland wrote:

> Hi all,
> 
> we are in the process of adding new functionalities to our parallel
> finite-element code and are facing new challenges in computing the number of
> non-zeros of a sparse matrix.  Since we use Petsc, we want to be able to give
> this information before the first assembly (to pass to
> Mat{MPI,SEQ}AIJSetPreallocation).
> 
> Willing to rewrite the nonzeros counting algorithm from scratch if necessary,
> we asked ourselves a number of questions:
> 
> 1) Do algorithms in general count the *exact* number of non-zeros?  or do they
> just count a maximum and let Petsc release the memory after the first
> assembly?
> 
> 2) How bad is it to reserve too-much non-zeros before the first assembly?

If you are counting extra nonzeros - it won't affect the performance
[as no additional mallocs will be required].

However the extra memory is not freed at the assembly time.

Satish

> 
> 3) Where could we find some hints/articles examples or information about that
> kind of algorithm (working in parallel of course) ?
> 
> Here is the context that we have up to now: we can manage a mixture of
> interpolations (P1, P2, discontinuous interpolations, etc) and of dimensions
> (scalar, 3D, etc) for all unknowns in a global matrix (example: a P2-P1
> [Taylor-Hood] mixed formulation for Stokes problem).
> 
> Since these parameters are determined at run-time (in fact they are user
> parameters), we have to think of a general algorithm to compute the "mostly"
> exact number of non-zeros for each matrix lines.
> 
> Right now, we have a working, but memory consuming algorithm which computes
> the *exact* number on non-zeros per line (we exploded the memory for a problem
> with a little more than 2^31 unknowns with it -- using 64 bits indices of
> course).  In fact, we manage to store the number of neighbors for each type of
> mesh components (ex: for a vertice, we have the following number of neighbors:
> vertices, edges, faces, and each type of elements).  It is costly and moreover
> will not handle the new functionalities we want in the code (XFEM for
> example).
> 
> Thanks a lot!
> 
> Eric
>