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. > 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 >> >>
signature.asc
Description: PGP signature
