Hi Barry, Am 17.07.2017 um 10:42 schrieb Barry Smith: > The intention with the AO is you map everything to the PETSc ordering and > then you just work in the PETSc ordering completely, this way there is no > "sending around with MPI preallocation information". That is your matrix > assembly and preallocation (as well as vectors) always works in the PETSc > numbering. In the PETSc ordering there is no need to send information between > processes about row allocations since the information is always on the > correct process.
The problem I have, is that the information needed to preallocate AO row 1 => PO row 8 is located at rank 0 of my application. So I either need to send the data around to be able to generate the prealloc information of PO row 8 on rank 1 OR send the preallocation data to prealloc PO row 8 on rank 1, albeit generated on rank 0. Is that correct? Best, Florian > > Barry > >> On Jul 16, 2017, at 9:26 PM, Florian Lindner <[email protected]> wrote: >> >> 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 <[email protected]> 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 >>>>> >>> >> >
