On Wed, Sep 16, 2015 at 7:18 AM, Lawrence Mitchell < [email protected]> wrote:
> -----BEGIN PGP SIGNED MESSAGE----- > Hash: SHA1 > > On 16/09/15 12:40, Matthew Knepley wrote: > > On Tue, Sep 15, 2015 at 9:05 PM, Adrian Croucher > > <[email protected] <mailto:[email protected]>> > > wrote: > > > > hi > > > > I have a test code (attached) that sets up a finite volume mesh > > using DMPlex, with 2 degrees of freedom per cell. > > > > I then create a matrix using DMCreateMatrix(), having used > > DMSetMatType() to set the matrix type to MATBAIJ or MATMPIBAIJ, to > > take advantage of the block structure. > > > > This works OK and gives me the expected matrix structure when I > > run on > 1 processor, but neither MATBAIJ or MATMPIBAIJ works if I > > run it in serial: > > > > > > Plex should automatically discover the block size from the Section. > > If not, it uses block size 1. I have to look at the example to see > > why the discovery is not working. Do you have any constraints? > > It looks like block discovery in parallel effectively always > determines a block size of 1. Running with -mat_view ::ascii_info gives: > > Mat Object: 2 MPI processes > type: mpibaij > rows=20, cols=20 > total: nonzeros=112, allocated nonzeros=112 > total number of mallocs used during MatSetValues calls =0 > block size is 1 > ^^^ > > The block size discovery does this: > > for (p = pStart; p < pEnd; ++p) { > ierr = PetscSectionGetDof(sectionGlobal, p, &dof);CHKERRQ(ierr); > ierr = PetscSectionGetConstraintDof(sectionGlobal, p, > &cdof);CHKERRQ(ierr); > if (dof-cdof) { > if (bs < 0) { > bs = dof-cdof; > } else if (bs != dof-cdof) { > /* Layout does not admit a pointwise block size */ > bs = 1; > break; > } > } > } > > In parallel, some of the dofs are remote (encoded as negative). So we > run through seeing (dof - cdof) == 2, until we hit a "remote" point at > when we see (dof - cdof) != 2 and then we break out and set bs = 1. > > In serial, we correctly determine bs == 2. But then > DMGetLocalToGlobalMapping always does > > ierr = ISLocalToGlobalMappingCreate(PETSC_COMM_SELF, 1,size, > ltog, PETSC_OWN_POINTER, &dm->ltogmap);CHKERRQ(ierr); > > > i.e. is set with block size == 1. > > So there are two bugs here. > > 1. In parallel, block size discovery in Mat creation has gone wrong > Crap. Hoist on my own petard. Okay I will fix this. > 2. (Always), the lgmap has block size of 1, irrespective of the > discovered block size from the section. > Yep. This can also be fixed. It should work regardless, but would be better with blocking. Thanks Matt > Lawrence > > > -----BEGIN PGP SIGNATURE----- > Version: GnuPG v1 > > iQEcBAEBAgAGBQJV+V4tAAoJECOc1kQ8PEYvqO4IAN4+oIgtBevvDAughPgUVOzq > kESJkb0Bx4a7+y47IPsY/SOuOMjVgErz2SO3tGd7+K/U5fstojJbIC7zZqPIERn0 > S68lH3s+y1pVqmcIMFIorz6is+u46M2xIGS6MLe6d0DluslyThaaA7lMhuIvJkIX > FC8QmtqCsIvHv10VuNll/81UJ3pXSZ+E81+Rs6pBhoGMCnRSXdMgFEtfWBBGL2JR > byEOAueTmY+YZXJ5JINxsHG1C1lep5wyYfAERDiRaD9bhCsE4mf9z/yFdvhiv0e3 > JC5WB3Q2/t2dRXhKNOSxQJd5mBKGhHgKptTjzmnW8HF0uxCdy7XpKZ6vS4Kt8Gk= > =2JQe > -----END PGP SIGNATURE----- > -- What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead. -- Norbert Wiener
