Jed has it pretty much right. See my SC2002 paper (my web page) or Prometheus source: src/gauss_seid.C
The local equation are partitioned into ~20 sets and GS is applied to these sets one at a time, with communication between them. It is completely asynchronous and it is a true G-S. The reverse stuff is need for symmetry so in the back sweep these ISs need to be processed in reverse order and they need to process the equation in reverse order. It somewhat complicated code but the control code is all in gauss_seid.C, I just have to move it to C (including a linked list). As far as PETSc design, the current SOR is not really SOR, it is block Jacobi with SOR local solver. So I think you are taking my name space! THe current BJ/SOR is very practical so I don't like making everyone change their command line args ... its up to you. Mark On Jun 9, 2012, at 9:50 PM, Jed Brown wrote: > I would be inclined to put the complexity in MatAIJ, but hopefully have it > reuse the special-purpose comm routines. > > On Jun 9, 2012 8:47 PM, "Barry Smith" <bsmith at mcs.anl.gov> wrote: > > On Jun 9, 2012, at 8:26 PM, Jed Brown wrote: > > > I think he wants to manage the messages at a higher level and call a matrix > > kernel for the partial sweeps. > > > > I would prefer to write the communication primitives so that GS for each > > matrix format is simple. > > I am not sure if this is the right model. But it is an alternative approach. > > Do you envision PCApply_SOR() managing all the stages of MPI sends and > receives? Or a new PCApply_PGS() managing that process. > > I admit I was thinking the Mat is responsible for its various SORs, but > other alternatives as you propose seem reasonable. > > From the user prospective I think it is better not bot have both PCSOR and > PCSPGS but to have options for PCSOR to make it it support PGS. How it is > handled internally with complicated PCApply_SOR or complicated > MatSOR_MPIAIJ() etc is more a detail I am less concerned about. > > > Barry > > > > > > On Jun 9, 2012 8:21 PM, "Barry Smith" <bsmith at mcs.anl.gov> wrote: > > > > On Jun 9, 2012, at 8:18 PM, Jed Brown wrote: > > > > > He is going to alternate between smoothing some points and sending > > > messages. > > > > Fine but that is all INSIDE a single SOR sweep? So I was wrong to say > > PCSORSetIS() maps to MatSORSetIS() it is PCSORSSetISs() maps to > > MatSORSetISs() he won't be calling MatSOR() for each piece but ONCE for the > > whole process (yes the MatSOR_SeqAIJ() is more complicated in this case. > > > > Barry > > > > > > > > On Jun 9, 2012 8:07 PM, "Barry Smith" <bsmith at mcs.anl.gov> wrote: > > > > > > On Jun 9, 2012, at 8:01 PM, Jed Brown wrote: > > > > > > > Parallel Gauss-Seidel. > > > > > > But if you know in advance the IS that you are providing (that > > > determines the order of the nodes smoothed) then why would you change it > > > the next iteration? That is, if you are providing the IS then it is in > > > no way asynchronous so that the fact that it is "parallel" Gauss-Seidel > > > doesn't affect the ordering. Hence I consider your response humorous but > > > non responsive :-) > > > > > > Barry > > > > > > > > > > > > > > > > > > > > On Jun 9, 2012 7:56 PM, "Barry Smith" <bsmith at mcs.anl.gov> wrote: > > > > > > > > On Jun 9, 2012, at 7:47 PM, Jed Brown wrote: > > > > > > > > > Fine, but I think Mark is going to change the IS every time MatSOR is > > > > > called. > > > > > > > > Surely not. What kind of weird-ass algorithm would that be? > > > > > > > > Barry > > > > > > > > > Either will work, but a separate call is awkward if it's not useful > > > > > to be persistent. > > > > > > > > > > On Jun 9, 2012 7:45 PM, "Barry Smith" <bsmith at mcs.anl.gov> wrote: > > > > > > > > > > On Jun 9, 2012, at 6:51 PM, Jed Brown wrote: > > > > > > > > > > > On Sat, Jun 9, 2012 at 6:43 PM, Mark F. Adams <mark.adams at > > > > > > columbia.edu> wrote: > > > > > > 1) I need a G-S kernel that takes an IS of indices to process and a > > > > > > flag to process them in forward or reverse order. How should I > > > > > > proceed to do this. Should I just clone sor? > > > > > > > > > > > > You are going to have several of these index sets? You could have a > > > > > > PCSORSetIS(). Probably need to add a MatOp for MatSORIS(). Barry > > > > > > might have other ideas. > > > > > > > > > > PCSORSetIS() would then go down to MatSORSetIS() and then the call > > > > > to MatSOR() would using the IS ordering if provided, otherwise use > > > > > the default natural ordering? > > > > > > > > > > I don't see a need to add a MatSORIS(). > > > > > > > > > > Barry > > > > > > > > > > > > > > > > > > > > > > 2) I don't want to use Richardson iterations for G-S. Should I > > > > > > make a G-S KPS method? I don't want to take a residual in the > > > > > > iterator (KSP) and if symmetric G-S is requested then it should > > > > > > drive this I think. > > > > > > > > > > > > Look at PCApplyRichardson_SOR(). > > > > > > > > > > > > SOR does two sweeps in each application; I'm not wild about that > > > > > > because a good way to run G-S in a V(1,1) cycle is to do a forward > > > > > > sweep in pre smoothing and a backward sweep in post smoothing. > > > > > > > > > > > > Well, MatSOR() has this flag MatSORType that can specify forward > > > > > > and reverse. You have one PC for the down-smoother and another for > > > > > > the up-smoother, then configure one to be a forward sweep and the > > > > > > other to be reverse. > > > > > > > > > > > > > > > -------------- next part -------------- An HTML attachment was scrubbed... URL: <http://lists.mcs.anl.gov/pipermail/petsc-dev/attachments/20120610/e47b7182/attachment.html>
