[petsc-users] Algorithms for counting nonzeros!

[Eric Chamberland]

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?

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