[petsc-users] block matrix in serial

[Adrian Croucher]

hi

Matt, did you (or anyone else) ever get a chance to look into these two bugs?

AFAICT the blocksize discovery code hasn't been changed since then.

I'm currently running a problem where the linear solver is struggling (running out of iterations),
and I'm wondering if it might perform better if the matrix blocksize were being set correctly.
May be nothing to do with it of course, but would at least be useful to eliminate that possibility.

Cheers, Adrian
//>//>/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 
<https://lists.mcs.anl.gov/mailman/listinfo/petsc-users> 
<mailto:a.croucher at auckland.ac.nz 
<https://lists.mcs.anl.gov/mailman/listinfo/petsc-users>>> />/ > 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

-- 
Dr Adrian Croucher
Senior Research Fellow
Department of Engineering Science
University of Auckland, New Zealand
email: a.croucher at auckland.ac.nz
tel: +64 (0)9 923 84611

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20160418/05eb8588/attachment.html>