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 <[email protected]> wrote: > > > On Aug 29, 2017, at 9:10 PM, Ben Yee <[email protected]> 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 <[email protected]> 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 <[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 > > -- Ben Yee NERS PhD Candidate, University of Michigan B.S. Mech. & Nuclear Eng., U.C. Berkeley
