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

[Ben Yee]

Oh, I could be mistaken, but I don't think the PCApply_PBJacobi_N()
function is in the master branch of petsc.  (At least, it's not in the
pbjacobi.c file that I downloaded for petsc-3.7.6)  I did a quick search
and found it in this old mat-taij branch: https://bitbucket.org/
petsc/petsc/branch/debo/mat-taij .  Is this the function you were referring
to?

Regarding the block density, I think the matrix is relatively diagonally
dominant, but the density is generally a bit higher than 30%.  The blocks
are of the same size.  I think you're right that the factorization is
faster, I will try the factorization first.  However, I am concerned about
the memory cost of the factorizing all of the blocks for problems with a
large number of equations, so I think I may have to resort to an iterative
solver or doing Gaussian elimination in those cases.

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

>
> > On Aug 29, 2017, at 9:10 PM, Ben Yee <bcyee at umich.edu> wrote:
> >
> > Thank you for your detailed response!
> >
> > In my case, the number of equations (the block size) depends on the user
> input (there is one equation per energy group, and the number of groups
> depends on how the user wants to discretize the energy variable).  It
> typically is around 50, but can be as high as a few hundred.  In
> pbjacobi.c, it seems that it is hard-coded to work for a fixed block size
> N, but I would like it to work for general N.  Moreover, I would like to
> solve the blocks iteratively (using SOR for example) since the block sizes
> can get rather large.
> >
> > I haven't tried this yet, but it seems that it wouldn't be too hard to
> modify what you have suggested to work as I have described above.  I could
> add an extra PetscBool to pbjacobi that indicates that I want to use my own
> PCApply_PBJacobi which would work for general block size N.   And, in that
> function, instead of applying a stored inverse block, I could implement SOR
> iterations.
>
>    There is a general PCApply_PBJacobi_N() you can start with no need for
> a bool.
>
>    Almost for sure if the blocks are more than 30% full you will do better
> to to do the factorization and inversion and NOT do SOR. How dense are your
> blocks?  I am assuming each block is of the same size?
>
>
>    Barry
>
> >
> > Does that sound like a good plan, or do you suggest an alternative
> approach?
> >
> > Thanks again for your help on this!
> >
> > On Tue, Aug 29, 2017 at 9:11 PM, Barry Smith <bsmith at mcs.anl.gov> wrote:
> >
> >   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);
> >   }
> >   }
> >
> >   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
>
>


-- 
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/20170829/e3fb430f/attachment-0001.html>