-----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 > <a.crouc...@auckland.ac.nz <mailto:a.crouc...@auckland.ac.nz>> > 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 2. (Always), the lgmap has block size of 1, irrespective of the discovered block size from the section. 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-----