[petsc-users] Understanding preallocation for MPI

[Florian Lindner]

Hello,

Am 11.07.2017 um 02:45 schrieb Barry Smith:
> 
>   You might consider using http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Mat/MatPreallocateInitialize.html and http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Mat/MatPreallocateSetLocal.html#MatPreallocateSetLocal and friends
> 
>    These take out some of the busywork of setting the preallocation arrays. They are macros in petscmat.h so even if you don't use them you can see the code they use.

Thanks for your input, I took a look at it and learned from it.

I have this preallocation routine which works perfectly for a index set on the columns and non (resp. an identity index
set) on the columns:

void doPreallocation(Mat &A, ISLocalToGlobalMapping &ISmapping, double nz_ratio, int nz_number)
{
  PetscInt d_nnz[local_rows], o_nnz[local_rows];

  const PetscInt *mapIndizes;
  ISLocalToGlobalMappingGetIndices(ISmapping, &mapIndizes);

  for (int row = row_range_start; row < row_range_end; row++) {
    // PetscInt relative_row = mapIndizes[row- row_range_start];
    PetscInt relative_row = row - row_range_start;
    d_nnz[relative_row] = 0;
    o_nnz[relative_row] = 0;
    int colStart, colEnd;
    std::tie(colStart, colEnd) = colBorders(row);
    for (int col = 0; col < global_cols; col++) {
      if (col >= colStart and col < colEnd) {
        int petsc_col = mapIndizes[col];
        if (petsc_col >= col_range_start and petsc_col < col_range_end) {
          d_nnz[relative_row]++;
        } else {
          o_nnz[relative_row]++;
        }
      }
    }
  }
  ierr = ISLocalToGlobalMappingRestoreIndices(ISmapping, &mapIndizes); CHKERRV(ierr);
  ierr = MatMPIAIJSetPreallocation(A, 0, d_nnz, 0, o_nnz); CHKERRV(ierr); // works only on minimum 2 ranks
}

colStart and colEnd give a defined number of diagonals for the matrix, just for testing. *_range_* are from the PETSc
routines.

But I'm puzzled how to do it when I have also a (the same) index set for the rows.

Let's say I have 2 procs, a global size of 10 and a mapping [0, 8, 9, 1, 3] on proc 0 and [2, 5, 4, 7, 6] on proc 1.

row_range is [0, 5] on the first proc.

On the first row-iteration I generate data for row 0 (application ordering AO) which is mapped to row 0 (petsc ordering
PO), I add these nnz to o_/d_nnz[0]

On the second row-iteration I generate data for row 1 (AO) which is mapped to row 8 (PO). As with the columns I expect
the nnz numbering to be in PO, which means either

1) nnz[8] (which is not existing), what I would be do in the commented out line.
2) nnz[4] on the second proc (which is unreachable from the first proc)
3) nnz[1] on the first proc, thus in AO, which produces new mallocs.

2. still seems to most probable for me, but before implementing sending the stuff around using MPI, I wanted to ask you
if it is correct that way? If it needs to be done that way, is there some PETSc magic which helps exchanging that array?

If the answer lies within MatPreallocate* routines, I failed to see :-/

Best Thanks,
Florian


>> On Jul 10, 2017, at 3:22 AM, Florian Lindner <mailinglists at xgm.de> wrote:
>>
>> Hey,
>>
>> one more question about preallocation:
>>
>> I can determine if a column index is diagonal or off-diagonal using that code
>>
>> if (col >= col_range_start and col < col_range_end)
>>    d_nnz[relative_row]++;
>> else
>>    o_nnz[relative_row]++;
>>
>>
>> My code, however uses index sets from which a ISlocalToGlobalMapping created:
>>
>>  // Creates a mapping from permutation and use that for the cols
>>  ierr = ISCreateGeneral(PETSC_COMM_WORLD, local_cols, permutation.data(), PETSC_COPY_VALUES, &ISlocal); CHKERRV(ierr);
>>  ierr = ISSetPermutation(ISlocal); CHKERRV(ierr);
>>  ierr = ISAllGather(ISlocal, &ISglobal); CHKERRV(ierr); // Gather the IS from all processors
>>  ierr = ISLocalToGlobalMappingCreateIS(ISglobal, &ISmapping); CHKERRV(ierr); // Make it a mapping
>>
>>  // Create an identity mapping and use that for the rows of A.
>>  ierr = ISCreateStride(PETSC_COMM_WORLD, local_rows, row_range_start, 1, &ISidentity); CHKERRV(ierr);
>>  ierr = ISSetIdentity(ISidentity); CHKERRV(ierr);
>>  ierr = ISAllGather(ISidentity, &ISidentityGlobal); CHKERRV(ierr);
>>  ierr = ISLocalToGlobalMappingCreateIS(ISidentityGlobal, &ISidentityMapping); CHKERRV(ierr);
>>
>>  ierr = MatSetLocalToGlobalMapping(A, ISidentityMapping, ISmapping); CHKERRV(ierr);
>>
>> since SetPreallocation routines define the diagonal / off-diagonal blocks from the petsc ordering, I have to translate
>> the col to a petsc_col.
>>
>> What is the best / fastest way to do that?
>>
>> Is that the way to go?
>>
>>  PetscInt mapSize;
>>  ISLocalToGlobalMappingGetSize(ISmapping, &mapSize);
>>  const PetscInt *mapIndizes;
>>  ISLocalToGlobalMappingGetIndices(ISmapping, &mapIndizes);
>>
>> Thanks,
>> Florian
>>
>>
>>
>> Am 07.07.2017 um 17:31 schrieb Florian Lindner:
>>> Hello,
>>>
>>> I'm having some struggle understanding the preallocation for MPIAIJ matrices, especially when a value is in off-diagonal
>>> vs. diagonal block.
>>>
>>> The small example program is at https://pastebin.com/67dXnGm3
>>>
>>> In general it should be parallel, but right now I just run it in serial.
>>>
>>> According to my understanding of
>>>
>>> http://www.mcs.anl.gov/petsc/petsc-3.7/docs/manualpages/Mat/MatMPIAIJSetPreallocation.html
>>>
>>> a entry is in the diagonal submatrix, if its row is in the OwnershipRange and its column is in OwnershipRangeColumn.
>>> That also means that in a serial run, there is only a diagonal submatrix.
>>>
>>> However, having MAT_NEW_NONZERO_ALLOCATION_ERR set, I get an error when
>>>
>>> Inserting 6 elements in row 2, though I have exactly
>>>
>>> 2 o_nnz = 0, d_nnz = 6 (means 6 elements allocated in the diagonal submatrix of row 2)
>>>
>>> Error is:
>>>
>>> [0]PETSC ERROR: Argument out of range
>>> [0]PETSC ERROR: New nonzero at (2,5) caused a malloc
>>> Use MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE) to turn off this check
>>>
>>>
>>> What is wrong with my understanding?
>>>
>>> Thanks,
>>> Florian
>>>
>