[petsc-users] Beginner questions : MatCreateMPIAIJWithSeqAIJ, MatCreateMPIAIJWithSplitArrays

[Mark Olesen]

I'm still at the beginning phase of looking at PETSc and accordingly 
have some beginner questions. My apologies if these are FAQs, but I 
didn't find much to address these specific questions.

My simulation matrices are sparse, and will generally (but not always) 
be generated in parallel. There is currently no conventional internal 
storage format (something like a COO variant), but lets just assume that 
I have CSR format for the moment.

I would like the usual combination of convenience and high efficiency, 
but efficiency (speed, memory) is the main criterion.

For serial, MatCreateSeqAIJWithArrays() looks like the thing to be 
using. It would provide a very thin wrapper around my CSR matrix without 
much additional allocation. The only extra allocation appears to be a 
precompute column range size per row (ilen) instead of doing it 
on-the-fly. If my matrix is actually to be considered symmetric, then 
use MatCreateSeqSBAIJWithArrays() instead.
This all seems pretty clear.


For parallel, MatCreateMPIAIJWithSplitArrays() appears to be the 
equivalent for efficiency, but I also read the note discouraging its 
use, which I fully appreciate. It also leads neatly into my question. I 
obviously will have fairly ready access to my on-processor portions of 
the matrix, but collecting the information for the off-processor 
portions is required. What would a normal or recommended approach look like?

For example,
====
Mat A = MatCreateSeqAIJWithArrays() to wrap the local CSR.

Mat B = MatCreateSeqAIJ(). Do some preallocation for num non-zeroes, use 
  MatSetValues() to fill in. Need extra garray[] as linear lookup for 
the global column numbers of B.

Or as an alternative, calculate the off-diagonal as a CSR by hand and 
use Mat B = MatCreateSeqAIJWithArrays() to wrap it.

Finally,
Use MatCreateMPIAIJWithSeqAIJ() to produce the full matrix.

Assuming that I used MatCreateSeqAIJWithArrays() to create both the A 
and B matrices, then they both hold a shallow copy of my own storage.
In MatCreateSeqAIJWithArrays(), I can't really tell what happens to the 
A matrix. For the B matrix, it appears that its column entries are 
changed to be those of the global columns and its data values are handed 
off to another MatCreateSeqAIJ() as the off-diagonal bits. The original 
B matrix is tagged up to avoid any deletion, and the shallow copied part 
is tagged to be deleted. If I follow this properly, this implies that if 
I was managing the storage of the original B matrix myself, I now have 
double deletion?

I would have expected something like this instead (around line 3431 of 
mpiaij.c in master):

   /* Retain original memory management */
   bnew->singlemalloc = b->singlemalloc;
   bnew->free_a       = b->free_a;
   bnew->free_ij      = b->free_ij;

   /* B arrays are shared by Bnew */
   b->singlemalloc = PETSC_FALSE;
   b->free_a       = PETSC_FALSE;
   b->free_ij      = PETSC_FALSE;
   ierr = MatDestroy(&B);CHKERRQ(ierr);


Have I gone off in completely the wrong direction here?
Is there a better method of approaching this?

Cheers,
/mark