[petsc-users] block matrix in serial

[Matthew Knepley]

On Wed, Sep 16, 2015 at 7:18 AM, Lawrence Mitchell <
lawrence.mitchell at imperial.ac.uk> 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
> > <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
>

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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20150916/0099906c/attachment.html>