[petsc-users] block matrix in serial

[Lawrence Mitchell]

-----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.croucher at auckland.ac.nz <mailto:a.croucher at 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-----