June 2012 - petsc-dev - mailman3.cels.anl.gov

Re: [petsc-dev] Parallel Gauss - Seidel
by Mark F. Adams 10 Jun '12

10 Jun '12

On Jun 9, 2012, at 7:51 PM, Jed Brown wrote: > On Sat, Jun 9, 2012 at 6:43 PM, Mark F. Adams <mark.adams(a)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? O(20) > You could have a PCSORSetIS(). Probably need to add a MatOp for MatSORIS(). Barry might have other ideas. > > 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. I can have two sets of ISs and avoid the reverse flag. I was avoiding that option but now that I think about it it simplifies the control code.

1 0

Re: [petsc-dev] How to dump matrix in petsc
by Jed Brown 10 Jun '12

10 Jun '12

On Sun, Jun 10, 2012 at 10:07 AM, Dave Nystrom <dnystrom1(a)comcast.net>wrote: > Is there an easy way to dump a matrix and rhs in petsc? I have a > particular > system that has been challenging to solve with petsc other than using > direct > solves. Also, if I could dump the linear system, is it acceptable to post > it > on this email list so that someone who knows more about petsc and linear > solvers could investigate it if they had the time and inclination? > Run with -ksp_view_binary (it produces a file called 'binaryoutput') Don't email large files to the list. You can send up to a few MB to petsc-maint. If you can't host larger files yourself, anonymous ftp to ftp.mcs.anl.gov, cd incoming, and upload the file (preferably with a unique name), then tell us the file name.

1 0

How to dump matrix in petsc
by Dave Nystrom 10 Jun '12

10 Jun '12

Is there an easy way to dump a matrix and rhs in petsc? I have a particular system that has been challenging to solve with petsc other than using direct solves. Also, if I could dump the linear system, is it acceptable to post it on this email list so that someone who knows more about petsc and linear solvers could investigate it if they had the time and inclination? Thanks, Dave

1 0

Re: [petsc-dev] Problem with VecLoad and .info files
by Blaise Bourdin 10 Jun '12

10 Jun '12

On Jun 10, 2012, at 9:06 AM, Jed Brown wrote: > Can you set a different prefix for each vector? (I'm surprised that more people don't run into this problem of incompatibility, but don't know an easy generic way to fix the problem.) > > $ cat vector.dat.info > -vecload_block_size 1 > -vecload_block_size 3 Very good point. That fixes the problem. and for existing files, I can always delete the .info file. I need to switch to the built-in binary vtk viewer anyway... Blaise > > On Sun, Jun 10, 2012 at 8:37 AM, Blaise Bourdin <bourdin(a)lsu.edu> wrote: > Hi, > > I think I found another bug related to Vec blocksize... > > It looks like loading vectors with different block sizes from a binary fine is broken in petsc-3.3: > try the attached example: > > iMac:Misc blaise$ ./TestVecLoad-dev -ndofU 1 -ndofV 3 -nx 2 -ny 2 -mode 0 > Vector Object:UVec 1 MPI processes > type: mpi > Process [0] > 0 > 10 > 1 > 11 > Vector Object:VVec 1 MPI processes > type: mpi > Process [0] > 0 > -100 > -200 > -10 > -110 > -210 > -1 > -101 > -201 > -11 > -111 > -211 > Writing vectors to file > iMac:Misc blaise$ ./TestVecLoad-dev -ndofU 1 -ndofV 3 -nx 2 -ny 2 -mode 1 > Reading vectors from file > Processor [0] M 2 N 2 m 1 n 1 w 1 s 1 > X range of indices: 0 2, Y range of indices: 0 2 > Processor [0] M 2 N 2 m 1 n 1 w 3 s 1 > X range of indices: 0 2, Y range of indices: 0 2 > [0]PETSC ERROR: --------------------- Error Message ------------------------------------ > [0]PETSC ERROR: Arguments are incompatible! > [0]PETSC ERROR: Local size 4 not compatible with block size 3! > [0]PETSC ERROR: ------------------------------------------------------------------------ > [0]PETSC ERROR: Petsc Release Version 3.3.0, Patch 0, unknown > [0]PETSC ERROR: See docs/changes/index.html for recent updates. > [0]PETSC ERROR: See docs/faq.html for hints about trouble shooting. > [0]PETSC ERROR: See docs/index.html for manual pages. > [0]PETSC ERROR: ------------------------------------------------------------------------ > [0]PETSC ERROR: ./TestVecLoad-dev on a Darwin-gc named iMac.local by blaise Sun Jun 10 08:35:07 2012 > [0]PETSC ERROR: Libraries linked from /opt/HPC/petsc-3.3/Darwin-gcc4.2-g/lib > [0]PETSC ERROR: Configure run at Fri Jun 8 22:01:00 2012 > [0]PETSC ERROR: Configure options CFLAGS=-Wno-unused --download-hdf5=1 --download-metis=1 --download-parmetis=1 --download-sowing=1 --download-triangle=1 --with-cmake=cmake --with-debugging=1 --with-gnu-compilers=1 --with-mpi-dir=/opt/HPC/mpich2-1.4.1p1-gcc4.2 --with-pic --with-shared-libraries=1 --with-x11=1 > [0]PETSC ERROR: ------------------------------------------------------------------------ > [0]PETSC ERROR: PetscLayoutSetBlockSize() line 459 in /opt/HPC/petsc-3.3/src/vec/vec/impls/mpi/pmap.c > [0]PETSC ERROR: VecSetBlockSize() line 1526 in /opt/HPC/petsc-3.3/src/vec/vec/interface/vector.c > [0]PETSC ERROR: VecLoad_Binary() line 104 in /opt/HPC/petsc-3.3/src/vec/vec/utils/vecio.c > [0]PETSC ERROR: VecLoad_Default() line 348 in /opt/HPC/petsc-3.3/src/vec/vec/utils/vecio.c > [0]PETSC ERROR: VecLoad() line 1111 in /opt/HPC/petsc-3.3/src/vec/vec/interface/vector.c > [0]PETSC ERROR: main() line 86 in TestVecLoad-dev.c > application called MPI_Abort(MPI_COMM_WORLD, 75) - process 0 > [unset]: aborting job: > application called MPI_Abort(MPI_COMM_WORLD, 75) - process 0 > > but if one destroys the vector.dat.info file, everything seems fine: > > iMac:Misc blaise$ ./TestVecLoad-dev -ndofU 1 -ndofV 3 -nx 2 -ny 2 -mode 0 > Vector Object:UVec 1 MPI processes > type: mpi > Process [0] > 0 > 10 > 1 > 11 > Vector Object:VVec 1 MPI processes > type: mpi > Process [0] > 0 > -100 > -200 > -10 > -110 > -210 > -1 > -101 > -201 > -11 > -111 > -211 > Writing vectors to file > iMac:Misc blaise$ rm vector.dat.info > iMac:Misc blaise$ ./TestVecLoad-dev -ndofU 1 -ndofV 3 -nx 2 -ny 2 -mode 1 > Reading vectors from file > Processor [0] M 2 N 2 m 1 n 1 w 1 s 1 > X range of indices: 0 2, Y range of indices: 0 2 > Processor [0] M 2 N 2 m 1 n 1 w 3 s 1 > X range of indices: 0 2, Y range of indices: 0 2 > Vector Object:UVec 1 MPI processes > type: seq > 0 > 10 > 1 > 11 > Vector Object:VVec 1 MPI processes > type: seq > 0 > -100 > -200 > -10 > -110 > -210 > -1 > -101 > -201 > -11 > -111 > -211 > > Blaise > > > > -- > Department of Mathematics and Center for Computation & Technology > Louisiana State University, Baton Rouge, LA 70803, USA > Tel. +1 (225) 578 1612, Fax +1 (225) 578 4276 http://www.math.lsu.edu/~bourdin > > > > > > > > > -- Department of Mathematics and Center for Computation & Technology Louisiana State University, Baton Rouge, LA 70803, USA Tel. +1 (225) 578 1612, Fax +1 (225) 578 4276 http://www.math.lsu.edu/~bourdin

1 0

Re: [petsc-dev] Problem with VecLoad and .info files
by Jed Brown 10 Jun '12

10 Jun '12

Can you set a different prefix for each vector? (I'm surprised that more people don't run into this problem of incompatibility, but don't know an easy generic way to fix the problem.) $ cat vector.dat.info -vecload_block_size 1 -vecload_block_size 3 On Sun, Jun 10, 2012 at 8:37 AM, Blaise Bourdin <bourdin(a)lsu.edu> wrote: > Hi, > > I think I found another bug related to Vec blocksize... > > It looks like loading vectors with different block sizes from a binary > fine is broken in petsc-3.3: > try the attached example: > > iMac:Misc blaise$ ./TestVecLoad-dev -ndofU 1 -ndofV 3 -nx 2 -ny 2 -mode 0 > Vector Object:UVec 1 MPI processes > type: mpi > Process [0] > 0 > 10 > 1 > 11 > Vector Object:VVec 1 MPI processes > type: mpi > Process [0] > 0 > -100 > -200 > -10 > -110 > -210 > -1 > -101 > -201 > -11 > -111 > -211 > Writing vectors to file > iMac:Misc blaise$ ./TestVecLoad-dev -ndofU 1 -ndofV 3 -nx 2 -ny 2 -mode 1 > Reading vectors from file > Processor [0] M 2 N 2 m 1 n 1 w 1 s 1 > X range of indices: 0 2, Y range of indices: 0 2 > Processor [0] M 2 N 2 m 1 n 1 w 3 s 1 > X range of indices: 0 2, Y range of indices: 0 2 > [0]PETSC ERROR: --------------------- Error Message > ------------------------------------ > [0]PETSC ERROR: Arguments are incompatible! > [0]PETSC ERROR: Local size 4 not compatible with block size 3! > [0]PETSC ERROR: > ------------------------------------------------------------------------ > [0]PETSC ERROR: Petsc Release Version 3.3.0, Patch 0, unknown > [0]PETSC ERROR: See docs/changes/index.html for recent updates. > [0]PETSC ERROR: See docs/faq.html for hints about trouble shooting. > [0]PETSC ERROR: See docs/index.html for manual pages. > [0]PETSC ERROR: > ------------------------------------------------------------------------ > [0]PETSC ERROR: ./TestVecLoad-dev on a Darwin-gc named iMac.local by > blaise Sun Jun 10 08:35:07 2012 > [0]PETSC ERROR: Libraries linked from > /opt/HPC/petsc-3.3/Darwin-gcc4.2-g/lib > [0]PETSC ERROR: Configure run at Fri Jun 8 22:01:00 2012 > [0]PETSC ERROR: Configure options CFLAGS=-Wno-unused --download-hdf5=1 > --download-metis=1 --download-parmetis=1 --download-sowing=1 > --download-triangle=1 --with-cmake=cmake --with-debugging=1 > --with-gnu-compilers=1 --with-mpi-dir=/opt/HPC/mpich2-1.4.1p1-gcc4.2 > --with-pic --with-shared-libraries=1 --with-x11=1 > [0]PETSC ERROR: > ------------------------------------------------------------------------ > [0]PETSC ERROR: PetscLayoutSetBlockSize() line 459 in > /opt/HPC/petsc-3.3/src/vec/vec/impls/mpi/pmap.c > [0]PETSC ERROR: VecSetBlockSize() line 1526 in > /opt/HPC/petsc-3.3/src/vec/vec/interface/vector.c > [0]PETSC ERROR: VecLoad_Binary() line 104 in > /opt/HPC/petsc-3.3/src/vec/vec/utils/vecio.c > [0]PETSC ERROR: VecLoad_Default() line 348 in > /opt/HPC/petsc-3.3/src/vec/vec/utils/vecio.c > [0]PETSC ERROR: VecLoad() line 1111 in > /opt/HPC/petsc-3.3/src/vec/vec/interface/vector.c > [0]PETSC ERROR: main() line 86 in TestVecLoad-dev.c > application called MPI_Abort(MPI_COMM_WORLD, 75) - process 0 > [unset]: aborting job: > application called MPI_Abort(MPI_COMM_WORLD, 75) - process 0 > > but if one destroys the vector.dat.info file, everything seems fine: > > iMac:Misc blaise$ ./TestVecLoad-dev -ndofU 1 -ndofV 3 -nx 2 -ny 2 -mode 0 > Vector Object:UVec 1 MPI processes > type: mpi > Process [0] > 0 > 10 > 1 > 11 > Vector Object:VVec 1 MPI processes > type: mpi > Process [0] > 0 > -100 > -200 > -10 > -110 > -210 > -1 > -101 > -201 > -11 > -111 > -211 > Writing vectors to file > iMac:Misc blaise$ rm vector.dat.info > iMac:Misc blaise$ ./TestVecLoad-dev -ndofU 1 -ndofV 3 -nx 2 -ny 2 -mode 1 > Reading vectors from file > Processor [0] M 2 N 2 m 1 n 1 w 1 s 1 > X range of indices: 0 2, Y range of indices: 0 2 > Processor [0] M 2 N 2 m 1 n 1 w 3 s 1 > X range of indices: 0 2, Y range of indices: 0 2 > Vector Object:UVec 1 MPI processes > type: seq > 0 > 10 > 1 > 11 > Vector Object:VVec 1 MPI processes > type: seq > 0 > -100 > -200 > -10 > -110 > -210 > -1 > -101 > -201 > -11 > -111 > -211 > > Blaise > > > > -- > Department of Mathematics and Center for Computation & Technology > Louisiana State University, Baton Rouge, LA 70803, USA > Tel. +1 (225) 578 1612, Fax +1 (225) 578 4276 > http://www.math.lsu.edu/~bourdin > > > > > > > > >

1 0

Problem with VecLoad and .info files
by Blaise Bourdin 10 Jun '12

10 Jun '12

Hi, I think I found another bug related to Vec blocksize... It looks like loading vectors with different block sizes from a binary fine is broken in petsc-3.3: try the attached example: iMac:Misc blaise$ ./TestVecLoad-dev -ndofU 1 -ndofV 3 -nx 2 -ny 2 -mode 0 Vector Object:UVec 1 MPI processes type: mpi Process [0] 0 10 1 11 Vector Object:VVec 1 MPI processes type: mpi Process [0] 0 -100 -200 -10 -110 -210 -1 -101 -201 -11 -111 -211 Writing vectors to file iMac:Misc blaise$ ./TestVecLoad-dev -ndofU 1 -ndofV 3 -nx 2 -ny 2 -mode 1 Reading vectors from file Processor [0] M 2 N 2 m 1 n 1 w 1 s 1 X range of indices: 0 2, Y range of indices: 0 2 Processor [0] M 2 N 2 m 1 n 1 w 3 s 1 X range of indices: 0 2, Y range of indices: 0 2 [0]PETSC ERROR: --------------------- Error Message ------------------------------------ [0]PETSC ERROR: Arguments are incompatible! [0]PETSC ERROR: Local size 4 not compatible with block size 3! [0]PETSC ERROR: ------------------------------------------------------------------------ [0]PETSC ERROR: Petsc Release Version 3.3.0, Patch 0, unknown [0]PETSC ERROR: See docs/changes/index.html for recent updates. [0]PETSC ERROR: See docs/faq.html for hints about trouble shooting. [0]PETSC ERROR: See docs/index.html for manual pages. [0]PETSC ERROR: ------------------------------------------------------------------------ [0]PETSC ERROR: ./TestVecLoad-dev on a Darwin-gc named iMac.local by blaise Sun Jun 10 08:35:07 2012 [0]PETSC ERROR: Libraries linked from /opt/HPC/petsc-3.3/Darwin-gcc4.2-g/lib [0]PETSC ERROR: Configure run at Fri Jun 8 22:01:00 2012 [0]PETSC ERROR: Configure options CFLAGS=-Wno-unused --download-hdf5=1 --download-metis=1 --download-parmetis=1 --download-sowing=1 --download-triangle=1 --with-cmake=cmake --with-debugging=1 --with-gnu-compilers=1 --with-mpi-dir=/opt/HPC/mpich2-1.4.1p1-gcc4.2 --with-pic --with-shared-libraries=1 --with-x11=1 [0]PETSC ERROR: ------------------------------------------------------------------------ [0]PETSC ERROR: PetscLayoutSetBlockSize() line 459 in /opt/HPC/petsc-3.3/src/vec/vec/impls/mpi/pmap.c [0]PETSC ERROR: VecSetBlockSize() line 1526 in /opt/HPC/petsc-3.3/src/vec/vec/interface/vector.c [0]PETSC ERROR: VecLoad_Binary() line 104 in /opt/HPC/petsc-3.3/src/vec/vec/utils/vecio.c [0]PETSC ERROR: VecLoad_Default() line 348 in /opt/HPC/petsc-3.3/src/vec/vec/utils/vecio.c [0]PETSC ERROR: VecLoad() line 1111 in /opt/HPC/petsc-3.3/src/vec/vec/interface/vector.c [0]PETSC ERROR: main() line 86 in TestVecLoad-dev.c application called MPI_Abort(MPI_COMM_WORLD, 75) - process 0 [unset]: aborting job: application called MPI_Abort(MPI_COMM_WORLD, 75) - process 0 but if one destroys the vector.dat.info file, everything seems fine: iMac:Misc blaise$ ./TestVecLoad-dev -ndofU 1 -ndofV 3 -nx 2 -ny 2 -mode 0 Vector Object:UVec 1 MPI processes type: mpi Process [0] 0 10 1 11 Vector Object:VVec 1 MPI processes type: mpi Process [0] 0 -100 -200 -10 -110 -210 -1 -101 -201 -11 -111 -211 Writing vectors to file iMac:Misc blaise$ rm vector.dat.info iMac:Misc blaise$ ./TestVecLoad-dev -ndofU 1 -ndofV 3 -nx 2 -ny 2 -mode 1 Reading vectors from file Processor [0] M 2 N 2 m 1 n 1 w 1 s 1 X range of indices: 0 2, Y range of indices: 0 2 Processor [0] M 2 N 2 m 1 n 1 w 3 s 1 X range of indices: 0 2, Y range of indices: 0 2 Vector Object:UVec 1 MPI processes type: seq 0 10 1 11 Vector Object:VVec 1 MPI processes type: seq 0 -100 -200 -10 -110 -210 -1 -101 -201 -11 -111 -211 Blaise -- Department of Mathematics and Center for Computation & Technology Louisiana State University, Baton Rouge, LA 70803, USA Tel. +1 (225) 578 1612, Fax +1 (225) 578 4276 http://www.math.lsu.edu/~bourdin

1 0

Something to keep in mind when discussing suggestions on petsc-dev :-)
by Barry Smith 10 Jun '12

10 Jun '12

http://www.sciencedaily.com/releases/2012/06/120607142246.htm Barry

1 0

Re: [petsc-dev] Parallel Gauss - Seidel
by Jed Brown 09 Jun '12

09 Jun '12

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(a)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(a)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(a)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(a)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(a)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(a)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. > > > > > > > > > > > > > > > >

1 0

Re: [petsc-dev] Parallel Gauss - Seidel
by Barry Smith 09 Jun '12

09 Jun '12

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(a)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(a)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(a)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(a)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(a)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. > > > > > > > > > >

1 0

Re: [petsc-dev] Parallel Gauss - Seidel
by Jed Brown 09 Jun '12

09 Jun '12

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. On Jun 9, 2012 8:21 PM, "Barry Smith" <bsmith(a)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(a)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(a)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(a)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(a)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. > > > > > > > > > > >

1 0