[petsc-users] Newbie Question - Really slow - PetscMemCpy

[Jed Brown]

On Thu, 22 Apr 2010 16:32:06 +0100, "Parker, Andrew (UK Filton)" <Andrew.Parker2 at baesystems.com> wrote:
> Well that worked.  Now it errors like nobodies business.  To be clear
> though, the sparsity pattern is just an int, indicating the number of
> non-zero blocks in a row, including the diagonal?  In my case the
> non-zero entries are due to entries from neighbouring cells, and myself,
> so can't see how I messed that up.  I know the answer is that I have,
> but just can't see why?  

The sparsity is an array with one int per row, indicating the number of
nonzero blocks in that row (including the diagonal).  You can provide
column indices by calling MatSetValues(Blocked) on each row, but that is
not necessary since it is cheap for the matrix implementation to shuffle
entries around within a row as you set them, the expensive part is if a
row overflows the space that has been allocated in which case it needs
to reallocate the entire thing (usually hundreds of MB) and copy over.

> My concern is that I've not understood what the sparsity is doing.  In
> my case I have a fully unstructured graph, so the neighbour cell numbers
> are in no way anywhere near in the number range to myself, the cell of
> interest.  Therefore, telling it that I have 10 neighbours, does not let
> it know the actual indices of those ten cells, which are not
> contiguous.

It's sparse storage, not banded storage, so this doesn't matter.  In
parallel, you have to specify how many entries are in the "diagonal
block" (the column index corresponds to a node that is owned) and how
many are in the "off-diagonal block" (corresponding to a node that is
owned by a remote process).

> Unless there's some funky mapping going on down in the guts.  Anyway, if
> the only number I need to give it is the number of non-zero blocks in a
> row, then I've done that and I think I need to rule that out and look
> for something else.

Can you make the problem size trivially small so that you can look at
the matrix to see what is different between the matrix you preallocated
for and the thing that actually gets assembled?

Note that if you preallocate for M nonzeros in a row, then assemble the
matrix after inserting only N < M entries in that row, then the matrix
will only keep space for N entries.  So if your first assembly is not
going to put something in all the slots that later assemblies will use,
you should assemble once with all the nonzeros you might need (setting
all the column indices, use MatSeqBAIJSetColumnIndices() if you have
them in an array).

Jed