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 <[email protected]> wrote: > > > On Aug 29, 2017, at 10:31 PM, Jed Brown <[email protected]> wrote: > > > > Barry Smith <[email protected]> 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 <[email protected]> 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" <[email protected]> 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 <[email protected]> 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
