[petsc-users] Beginner questions : MatCreateMPIAIJWithSeqAIJ, MatCreateMPIAIJWithSplitArrays

[Smith, Barry F.]

   Mark,

     I think you are over-optimizing your matrix assembly leading to complicated, fragile code. Better just to create the matrix and use MatSetValues() to set values into the matrix and not to work directly with the various sparse matrix data structures. If you wish to work directly with the sparse matrix data structures then you are largely on your own and need to figure out yourself how to use them. Plus you will only get a small benefit time wise in going the much more complicated route.

    Barry


> On Jul 24, 2018, at 6:52 AM, Mark Olesen <Mark.Olesen at esi-group.com> wrote:
> 
> 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