<div dir="ltr">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.<div><br></div><div>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?)<br><br>Also, could you provide a brief overview of how I should go about altering SOR to do the coloring? Thanks!</div><div><br></div><div><br></div></div><div class="gmail_extra"><br><div class="gmail_quote">On Tue, Aug 29, 2017 at 11:35 PM, Barry Smith <span dir="ltr"><<a href="mailto:bsmith@mcs.anl.gov" target="_blank">bsmith@mcs.anl.gov</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div class="HOEnZb"><div class="h5"><br>
> On Aug 29, 2017, at 10:31 PM, Jed Brown <<a href="mailto:jed@jedbrown.org">jed@jedbrown.org</a>> wrote:<br>
><br>
> Barry Smith <<a href="mailto:bsmith@mcs.anl.gov">bsmith@mcs.anl.gov</a>> writes:<br>
><br>
>> Ok, you should be using the BAIJ matrix it supports point-block Jacobi and point-block Gauss-Seidel.<br>
>><br>
>> 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/<wbr>pbjacobi.c<br>
>><br>
>> 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<br>
>><br>
>> PCPBJacobiSetColors(PC pc,PetscInt number of colors, PetscInt *sizeofeachcolor, PetscInt **idsforeachcolor);<br>
>><br>
>> you will have to add additional fields in PC_PBJacobi to contain this information.<br>
>><br>
>> 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<br>
>><br>
>> static PetscErrorCode PCApply_PBJacobi_2_Color(PC pc,Vec x,Vec y)<br>
>> {<br>
>> PC_PBJacobi *jac = (PC_PBJacobi*)pc->data;<br>
>> PetscErrorCode ierr;<br>
>> PetscInt i,m = jac->mbs;<br>
>> const MatScalar *diag = jac->diag;<br>
>> PetscScalar x0,x1,*yy;<br>
>> const PetscScalar *xx;<br>
>><br>
>> PetscFunctionBegin;<br>
>> if (!jac->b) {<br>
>> ierr = VecDuplicate(x,&jac->b);<br>
>> ierr = VecDuplicate(x,&jac->work);<br>
>> }<br>
>><br>
>> ierr = VecCopy(x,b);CHKERRQ(ierr);<br>
>> for (j=0; j<jac->numberofcolors; j++) {<br>
>><br>
>> ierr = VecGetArrayRead(b,&xx);<wbr>CHKERRQ(ierr);<br>
>> ierr = VecGetArray(y,&yy);CHKERRQ(<wbr>ierr);<br>
>><br>
>> for (i=0; i<jac->sizeofeachcolor[j]; i++) {<br>
>> ii = jac->idsforeachcolor[j][i];<br>
>> diag = jac->diag + 4*ii;<br>
>> x0 = xx[2*ii]; x1 = xx[2*ii+1];<br>
>> yy[2*ii] = diag[0]*x0 + diag[2]*x1;<br>
>> yy[2*ii+1] = diag[1]*x0 + diag[3]*x1;<br>
>> }<br>
>> ierr = VecRestoreArrayRead(b,&xx);<wbr>CHKERRQ(ierr);<br>
>> ierr = VecRestoreArray(y,&yy);<wbr>CHKERRQ(ierr);<br>
>><br>
>> /* update residual */<br>
>> if (i < jac->sizeofeachcolor[j]-1) {<br>
>> ierr = MatMult(pc->matrix,y,work2);<br>
>> ierr = VecAYPX(b,-1,work1);<br>
>> }<br>
>> }<br>
><br>
> Why this way instead of adding the functionality to SOR? This global<br>
> residual update is horribly inefficient compared to merely computing the<br>
> residual on each color.<br>
<br>
</div></div> 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.<br>
<span class="HOEnZb"><font color="#888888"><br>
Barry<br>
</font></span><div class="HOEnZb"><div class="h5"><br>
><br>
>> PetscFunctionReturn(0);<br>
>> }<br>
>><br>
>> Finally in PCSetUp_PBJacobi() you would set the apply function pointer to the "color" version if the user has provided the coloring information.<br>
>><br>
>> Pretty simple.<br>
>><br>
>> Barry<br>
>><br>
>><br>
>>> On Aug 29, 2017, at 6:47 PM, Ben Yee <<a href="mailto:bcyee@umich.edu">bcyee@umich.edu</a>> wrote:<br>
>>><br>
>>> 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.<br>
>>><br>
>>> 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.)<br>
>>><br>
>>> 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).<br>
>>><br>
>>> On Aug 29, 2017 7:33 PM, "Barry Smith" <<a href="mailto:bsmith@mcs.anl.gov">bsmith@mcs.anl.gov</a>> wrote:<br>
>>><br>
>>> Ben,<br>
>>><br>
>>> 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.<br>
>>><br>
>>> Barry<br>
>>><br>
>>>> On Aug 29, 2017, at 6:19 PM, Ben Yee <<a href="mailto:bcyee@umich.edu">bcyee@umich.edu</a>> wrote:<br>
>>>><br>
>>>> Hi,<br>
>>>><br>
>>>> 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?<br>
>>>><br>
>>>> If it is not possible, what would be the best way to implement this? Would I use KSPRICHARDSON+PCSHELL?<br>
>>>><br>
>>>> Thanks!<br>
>>>><br>
>>>> --<br>
>>>> Ben Yee<br>
>>>><br>
>>>> NERS PhD Candidate, University of Michigan<br>
>>>> B.S. Mech. & Nuclear Eng., U.C. Berkeley<br>
>>><br>
>>><br>
<br>
</div></div></blockquote></div><br><br clear="all"><div><br></div>-- <br><div class="gmail_signature" data-smartmail="gmail_signature"><div dir="ltr"><div><div dir="ltr">Ben Yee<div><font size="1"><br></font><div><font size="1">NERS PhD Candidate, University of Michigan</font></div><div><font size="1">B.S. Mech. & Nuclear Eng., U.C. Berkeley</font></div></div></div></div></div></div>
</div>