<div dir="ltr"><div class="gmail_extra"><div class="gmail_quote">On Sun, Apr 17, 2016 at 11:31 PM, Adrian Croucher <span dir="ltr"><<a href="mailto:a.croucher@auckland.ac.nz" target="_blank">a.croucher@auckland.ac.nz</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
  
    
  
  <div bgcolor="#FFFFFF" text="#000000">
    hi<br>
    <pre>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.</pre></div></blockquote><div><br></div><div>I have finally pushed what I think are the right fixes into master. Can you take a look and let me know</div><div>if it is fixed for you?</div><div><br></div><div>  Thanks,</div><div><br></div><div>     Matt</div><div> </div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div bgcolor="#FFFFFF" text="#000000"><pre>Cheers, Adrian
<span class=""><i>
</i>><i>
</i>><i> On 16/09/15 12:40, Matthew Knepley wrote:
</i></span><span class="">><i> > On Tue, Sep 15, 2015 at 9:05 PM, Adrian Croucher
</i></span>><i> > <<a href="https://lists.mcs.anl.gov/mailman/listinfo/petsc-users" target="_blank">a.croucher at auckland.ac.nz</a> <mailto:<a href="https://lists.mcs.anl.gov/mailman/listinfo/petsc-users" target="_blank">a.croucher at auckland.ac.nz</a>>>
</i><span class="">><i> > wrote:
</i>><i> >
</i>><i> > hi
</i>><i> >
</i>><i> > I have a test code (attached) that sets up a finite volume mesh
</i>><i> > using DMPlex, with 2 degrees of freedom per cell.
</i>><i> >
</i>><i> > I then create a matrix using DMCreateMatrix(), having used
</i>><i> > DMSetMatType() to set the matrix type to MATBAIJ or MATMPIBAIJ, to
</i>><i> > take advantage of the block structure.
</i>><i> >
</i>><i> > This works OK and gives me the expected matrix structure when I
</i>><i> > run on > 1 processor, but neither MATBAIJ or MATMPIBAIJ works if I
</i>><i> > run it in serial:
</i>><i> >
</i>><i> >
</i>><i> > Plex should automatically discover the block size from the Section.
</i>><i> > If not, it uses block size 1. I have to look at the example to see
</i>><i> > why the discovery is not working. Do you have any constraints?
</i>><i>
</i></span><div><div class="h5">><i> It looks like block discovery in parallel effectively always
</i>><i> determines a block size of 1.  Running with -mat_view ::ascii_info gives:
</i>><i>
</i>><i> Mat Object: 2 MPI processes
</i>><i>   type: mpibaij
</i>><i>   rows=20, cols=20
</i>><i>   total: nonzeros=112, allocated nonzeros=112
</i>><i>   total number of mallocs used during MatSetValues calls =0
</i>><i>       block size is 1
</i>><i>                     ^^^
</i>><i>
</i>><i> The block size discovery does this:
</i>><i>
</i>><i>         for (p = pStart; p < pEnd; ++p) {
</i>><i>           ierr = PetscSectionGetDof(sectionGlobal, p, &dof);CHKERRQ(ierr);
</i>><i>           ierr = PetscSectionGetConstraintDof(sectionGlobal, p,
</i>><i> &cdof);CHKERRQ(ierr);
</i>><i>           if (dof-cdof) {
</i>><i>             if (bs < 0) {
</i>><i>               bs = dof-cdof;
</i>><i>             } else if (bs != dof-cdof) {
</i>><i>               /* Layout does not admit a pointwise block size */
</i>><i>               bs = 1;
</i>><i>               break;
</i>><i>             }
</i>><i>           }
</i>><i>         }
</i>><i>
</i>><i> In parallel, some of the dofs are remote (encoded as negative).  So we
</i>><i> run through seeing (dof - cdof) == 2, until we hit a "remote" point at
</i>><i> when we see (dof - cdof) != 2 and then we break out and set bs = 1.
</i>><i>
</i>><i> In serial, we correctly determine bs == 2.  But then
</i>><i> DMGetLocalToGlobalMapping always does
</i>><i>
</i>><i>       ierr = ISLocalToGlobalMappingCreate(PETSC_COMM_SELF, 1,size,
</i>><i> ltog, PETSC_OWN_POINTER, &dm->ltogmap);CHKERRQ(ierr);
</i>><i>
</i>><i>
</i>><i> i.e. is set with block size == 1.
</i>><i>
</i>><i> So there are two bugs here.
</i>><i>
</i>><i> 1. In parallel, block size discovery in Mat creation has gone wrong
</i>><i>
</i>
Crap. Hoist on my own petard. Okay I will fix this.


><i> 2. (Always), the lgmap has block size of 1, irrespective of the
</i>><i> discovered block size from the section.
</i>><i>
</i>
Yep. This can also be fixed. It should work regardless, but would be better
with blocking.

  Thanks

     Matt</div></div></pre><span class="HOEnZb"><font color="#888888">
    </font></span><pre cols="72"><span class="HOEnZb"><font color="#888888">-- 
Dr Adrian Croucher
Senior Research Fellow
Department of Engineering Science
University of Auckland, New Zealand
email: <a href="mailto:a.croucher@auckland.ac.nz" target="_blank">a.croucher@auckland.ac.nz</a></font></span><span class="">
tel: +64 (0)9 923 84611
</span></pre>
  </div>

</blockquote></div><br><br clear="all"><div><br></div>-- <br><div class="gmail_signature">What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.<br>-- Norbert Wiener</div>
</div></div>