[petsc-users] Red-black block Jacobi in PETSc?

[Ben Yee]

I am interested in implementing it in the SOR framework.  My ultimate goal
is to find something that performs well for my problem, and I think I have
enough evidence from my simple 1D codes that it should work.  If I'm going
to assess performance, I think it would be good to implement it reasonably
optimally.

The manual page for PCSOR says that it is "Not a true parallel SOR, in
parallel this implementation corresponds to block Jacobi with SOR on each
block".  If I'm reading this correctly, then that description sounds like
what I want to do except that it is missing the coloring scheme -- is that
correct?  If that is correct, then is there a way to define a block size on
a MPIAIJ matrix, or would I have to use a BAIJ matrix?  (With BJACOBI,
there's the set total number of blocks function for defining the block
sizes.  Is there something equivalent for PCSOR?)

Also, could you provide a brief overview of how I should go about altering
SOR to do the coloring?  Thanks!



On Tue, Aug 29, 2017 at 11:35 PM, Barry Smith <bsmith at mcs.anl.gov> wrote:

>
> > On Aug 29, 2017, at 10:31 PM, Jed Brown <jed at jedbrown.org> wrote:
> >
> > Barry Smith <bsmith at mcs.anl.gov> writes:
> >
> >>  Ok, you should be using the BAIJ matrix it supports point-block Jacobi
> and point-block Gauss-Seidel.
> >>
> >>   We do not have a red-black Jacobi/Gauss-Seidel but you could easily
> add it. You will add it, not by using a any shell objects but by adding
> functionality to the PCPBJACOBI code in PETSc which is in
> src/ksp/pc/impls/pbjacobi/pbjacobi.c
> >>
> >>   First you will need to add a routine to supply the colors (make it
> general for any number of colors since the code is basically the same as
> for only two colors) call it say
> >>
> >>   PCPBJacobiSetColors(PC pc,PetscInt number of colors, PetscInt
> *sizeofeachcolor, PetscInt **idsforeachcolor);
> >>
> >>  you will have to add additional fields in PC_PBJacobi to contain this
> information.
> >>
> >>  Then copy the static PetscErrorCode PCApply_PBJacobi_N(PC pc,Vec x,Vec
> y)  for your block size N (for example 3) and modify it to do what you
> want, so for example
> >>
> >> static PetscErrorCode PCApply_PBJacobi_2_Color(PC pc,Vec x,Vec y)
> >> {
> >>  PC_PBJacobi     *jac = (PC_PBJacobi*)pc->data;
> >>  PetscErrorCode  ierr;
> >>  PetscInt        i,m = jac->mbs;
> >>  const MatScalar *diag = jac->diag;
> >>  PetscScalar     x0,x1,*yy;
> >>  const PetscScalar *xx;
> >>
> >>  PetscFunctionBegin;
> >>  if (!jac->b) {
> >>     ierr = VecDuplicate(x,&jac->b);
> >>     ierr = VecDuplicate(x,&jac->work);
> >>  }
> >>
> >>  ierr = VecCopy(x,b);CHKERRQ(ierr);
> >>  for (j=0; j<jac->numberofcolors; j++) {
> >>
> >>  ierr = VecGetArrayRead(b,&xx);CHKERRQ(ierr);
> >>  ierr = VecGetArray(y,&yy);CHKERRQ(ierr);
> >>
> >>  for (i=0; i<jac->sizeofeachcolor[j]; i++) {
> >>    ii = jac->idsforeachcolor[j][i];
> >>    diag = jac->diag + 4*ii;
> >>    x0        = xx[2*ii]; x1 = xx[2*ii+1];
> >>    yy[2*ii]   = diag[0]*x0 + diag[2]*x1;
> >>    yy[2*ii+1] = diag[1]*x0 + diag[3]*x1;
> >>  }
> >>  ierr = VecRestoreArrayRead(b,&xx);CHKERRQ(ierr);
> >>  ierr = VecRestoreArray(y,&yy);CHKERRQ(ierr);
> >>
> >>  /* update residual */
> >>  if (i < jac->sizeofeachcolor[j]-1) {
> >>     ierr = MatMult(pc->matrix,y,work2);
> >>     ierr = VecAYPX(b,-1,work1);
> >>  }
> >>  }
> >
> > Why this way instead of adding the functionality to SOR?  This global
> > residual update is horribly inefficient compared to merely computing the
> > residual on each color.
>
>    True, but if this is for theoretical studies of convergence and not
> optimal performance it is easier. For example you do this and see the
> convergence is much better with coloring than without and then conclude it
> is worth spending the time to do it properly within an SOR framework.
>
>    Barry
>
> >
> >>  PetscFunctionReturn(0);
> >> }
> >>
> >>  Finally in PCSetUp_PBJacobi() you would set the apply function pointer
> to the "color" version if the user has provided the coloring information.
> >>
> >>  Pretty simple.
> >>
> >>  Barry
> >>
> >>
> >>> On Aug 29, 2017, at 6:47 PM, Ben Yee <bcyee at umich.edu> wrote:
> >>>
> >>> I'm solving a coupled set of equations, so each "block" corresponds to
> a set of unknowns for a particular spatial cell.  The matrix is structured
> such that all of the unknowns for a given spatial cell have adjacent global
> matrix indices (i.e., they're next to each other in the global solution
> vector).  Effectively, I want to do red-black Gauss Seidel, but with
> blocks.  Alternatively, it's the same as applying block Jacobi for all the
> red cells and then applying block Jacobi for all the black cells.
> >>>
> >>> The color of the block is determined from the geometry of the problem
> which is stored in various structures in the code I'm working on,
> independent of petsc.  (Physically, I generally have a nice 3d cartesian
> spatial grid and the coloring is just a checkerboard in that case.)
> >>>
> >>> The reason I want to do this is for research purposes.  I've
> implemented my own interpolation matrices for PCMG, and, in my simpler 1d
> codes and convergence analyses, I've found that doing a red-black smoothing
> significantly improved convergence for my particular problem (though I'm
> aware that this generally leads to poor cache efficiency).
> >>>
> >>> On Aug 29, 2017 7:33 PM, "Barry Smith" <bsmith at mcs.anl.gov> wrote:
> >>>
> >>>  Ben,
> >>>
> >>>   Please explain more what you mean by "a red-black block Jacobi
> smoother". What is your matrix structure? What are the blocks? How do you
> decide which ones are which color?  Why do you wish to use some a smoother.
> >>>
> >>>  Barry
> >>>
> >>>> On Aug 29, 2017, at 6:19 PM, Ben Yee <bcyee at umich.edu> wrote:
> >>>>
> >>>> Hi,
> >>>>
> >>>> For the smoother in PCMG, I want to use a red-black block Jacobi
> smoother.  Is this available with the existing PETSc options?  If so, how
> do I designate which blocks are red and which are black?
> >>>>
> >>>> If it is not possible, what would be the best way to implement this?
> Would I use KSPRICHARDSON+PCSHELL?
> >>>>
> >>>> Thanks!
> >>>>
> >>>> --
> >>>> Ben Yee
> >>>>
> >>>> NERS PhD Candidate, University of Michigan
> >>>> B.S. Mech. & Nuclear Eng., U.C. Berkeley
> >>>
> >>>
>
>


-- 
Ben Yee

NERS PhD Candidate, University of Michigan
B.S. Mech. & Nuclear Eng., U.C. Berkeley
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20170830/e057d4f8/attachment.html>