[petsc-dev] Dealing with off-processor aij-matrices not having preallocation information

[Karl Rupp]

Hi guys,

if a distributed aij-matrix is created and no preallocation hints are 
given, the current behavior is to silently allow preallocation even if 
the user requested to error out if not preallocated.
 From MatMPIAIJSetPreallocation_MPIAIJ() in mpiaij.c:3285:
(slightly modified using style violations for better visual appearance)


   if (d_nz == PETSC_DEFAULT || d_nz == PETSC_DECIDE) d_nz = 5;
   if (o_nz == PETSC_DEFAULT || o_nz == PETSC_DECIDE) o_nz = 2;

   /* Quite some more code here */

   ierr = MatSeqAIJSetPreallocation(b->A,d_nz,d_nnz);CHKERRQ(ierr);
   ierr = MatSeqAIJSetPreallocation(b->B,o_nz,o_nnz);CHKERRQ(ierr);

   /* Do not error if the user did not give real preallocation 
information. Ugly because this would overwrite a previous user call to 
MatSetOption(). */
   if (!d_realalloc)
     MatSetOption(b->A,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);
   if (!o_realalloc)
     MatSetOption(b->B,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);

As the comment noted, this is quite a hack and actually breaks the 
current implementation for GPU matrices. I think that the problem can be 
completely resolved by not checking for PETSC_DEFAULT and PETSC_DECIDE 
in the MPI preallocation above, but just forward d_nz and o_nz to the 
sequential preallocation routine. This way, the logic remains in 
MatSeqAIJSetPreallocation(), where MAT_NEW_NONZERO_ALLOCATION_ERR is 
only set if the user provided an explicit preallocation scheme.

Do I miss something?

Best regards,
Karli