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

[Jed Brown]

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.

>   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
>> 
>> 
-------------- next part --------------
A non-text attachment was scrubbed...
Name: signature.asc
Type: application/pgp-signature
Size: 832 bytes
Desc: not available
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20170829/11d31046/attachment.sig>