February 2013 - petsc-dev - mailman3.cels.anl.gov

Re: [petsc-dev] Deflated Krylov solvers for PETSc
by Jie Chen 28 Feb '13

28 Feb '13

I guess viewing deflation as a preconditioner may destroy the symmetry of the system if it is originally symmetric? I wrote a one page technical details about deflation and it may help clarify my thoughts. http://www.mcs.anl.gov/~jiechen/tmp/deflate.pdf Jie ----- Original Message ----- From: "Barry Smith" <bsmith(a)mcs.anl.gov> To: "For users of the development version of PETSc" <petsc-dev(a)mcs.anl.gov> Sent: Wednesday, February 27, 2013 2:32:38 PM Subject: Re: [petsc-dev] Deflated Krylov solvers for PETSc The PETSc Krylov methods are supposed to be implemented using KSP_MatMult(), KSP_PCApply() etc instead of directly using MatMult() and PCApply() etc. This means, in theory, one can add support for techniques that modify the application of the multiply or preconditioner in a single place and does not mean changing code in each method implementation. I've always had trouble wrapping my head around if deflation techniques should be viewed as modifications to the Krylov methods or modification (via composing) of the preconditioners. Having this completely clarified would help us understand how deflation could be added to PETSc in a systematic (and hence maintainable) way. Barry On Feb 26, 2013, at 10:58 PM, Jie Chen <jiechen(a)mcs.anl.gov> wrote: > To whomever is interested in deflated Krylov solvers: > > Deflation is very useful technique for accelerating the convergence of Krylov iterations. It seems to me that deflation is not supported in PETSc, except for the special case DGMRES, where eigenvectors are used to deflate GMRES. As a matter of fact, one does not need to use eigenvectors for deflation, nor deflation is dependent on a specific choice of the iterative method. > > Coincidentally, Kees Vuik pointed to me an earlier paper of his: > > http://ta.twi.tudelft.nl/nw/users/vuik/papers/Vui01F.pdf > > which claims that the convergence of using block Jacobi/ILU(0) [the default preconditioner for PETSc] becomes in general better if deflation is included. In his paper, deflation means "coarse grid acceleration"; see page 935, the paragraph of equation (2). > > So I think (including my personal experience), deflation is a general enough technique that may benefit every user of the KSP solver. Implementation-wise, my personal take is that the library needs to change the A-multiply to (A-AZ(Z'AZ)^{-1}Z'A)-multiply, where Z is the deflation matrix: tall, skinny, and with a small number of columns. PETSc can ask the user to define Z, and it can also predefine several choices of Z, one example choice is from the above mentioned paper. In principle, deflation should be independent of the KSP types, but I don't see how exactly it can be implemented without touching the code of every krylov solvers. Perhaps adding one more layer of abstraction is needed. > > So for anyone who is interested in implementing deflation, I can provide a more detailed tutorial. The design of PETSc is becoming more and more complicated for me to digest, thus I think it may be unrealistic for me to implement deflation myself for the moment. > > Jie > > > > -- > Jie Chen > Mathematics and Computer Science Division > Argonne National Laboratory > Address: 9700 S Cass Ave, Bldg 240, Lemont, IL 60439 > Phone: (630) 252-3313 > Email: jiechen(a)mcs.anl.gov > Homepage: http://www.mcs.anl.gov/~jiechen

1 0

Re: [petsc-dev] [petsc-users] MATNEST with shell matrices
by Jed Brown 28 Feb '13

28 Feb '13

On Thu, Feb 28, 2013 at 12:58 PM, Amneet Bhalla <mail2amneet(a)gmail.com>wrote: > > Hi Jed, > > I have registered MatShell that operates on VecNest (which is composed of > shell vectors) with KSP. > When I am solving sytem of equations using this KSP, I am getting some > errors. The sequence of > events from gdb are: > > #0 0x0000000000000000 in ?? () > > #1 0x00000000009856fd in VecGetArrayRead (x=0x3422a30, a=0x7fff8b8c6620) > at > /home/amneet/SOFTWARES/LIBRARIESIBAMR/PETSc-3.3-Lite/include/petsc-private/vecimpl.h:350 > > #2 0x000000000098f57f in VecGetArray_Nest (X=0x34291d0, x=0x7fff8b8c6720) > at > /home/amneet/SOFTWARES/LIBRARIESIBAMR/PETSc-3.3-Lite/src/vec/vec/impls/nest/vecnest.c:719 > > #3 0x0000000000a3f3df in VecGetArrayRead (x=0x34291d0, a=0x7fff8b8c6720) > at > /home/amneet/SOFTWARES/LIBRARIESIBAMR/PETSc-3.3-Lite/include/petsc-private/vecimpl.h:350 > This is calling VecValidValues() checks for vec->ops->getarray. Since VecGetArray_Nest() is implemented, it calls the routine. I'm inclined to add a vec->ops->validvalues() and _not_ try to hack it in terms of VecGetArray(). > #4 0x0000000000a51c83 in MatMult (mat=0x32c34f0, x=0x336a1c0, y=0x34291d0) > at > /home/amneet/SOFTWARES/LIBRARIESIBAMR/PETSc-3.3-Lite/src/mat/interface/matrix.c:2092 > > #5 0x00000000012b2a78 in KSPInitialResidual (ksp=0x334bcb0, > vsoln=0x336a1c0, > vt1=0x34291d0, vt2=0x343c150, vres=0x344f220, vb=0x337d810) > at > /home/amneet/SOFTWARES/LIBRARIESIBAMR/PETSc-3.3-Lite/src/ksp/ksp/interface/itres.c:54 > > #6 0x0000000000da2eaa in KSPSolve_GMRES (ksp=0x334bcb0) > at > /home/amneet/SOFTWARES/LIBRARIESIBAMR/PETSc-3.3-Lite/src/ksp/ksp/impls/gmres/gmres.c:230 > > #7 0x0000000000d1ad6b in KSPSolve (ksp=0x334bcb0, b=0x337d810, > x=0x336a1c0) > at > /home/amneet/SOFTWARES/LIBRARIESIBAMR/PETSc-3.3-Lite/src/ksp/ksp/interface/itfunc.c:446 > > #8 0x00000000006b8a08 in IBTK::PETScBlockKrylovLinearSolver::solveSystem ( > this=0x3249380, x=0x32f3e20, b=0x330c6d0) > at PETScBlockKrylovLinearSolver.C:384 > > #9 0x00000000005777e3 in main (argc=6, argv=0x7fff8b8ca988) at main.C:274 > > In frame 4, the user defined matrix-vector multiplication is successful ( > I printed out a statement). > This would imply in the routine MatVecMult(Mat A, Vec x, Vec y), Vec y > would be a duplicated > vector. This in turn would mean that user-defined routine to duplicate > vectors MatGetVecs(Mat A, Vec *left, Vec* right) > should have been called. But this is not happening. It is in this routine, > that I am setting > > *left->ops->getarray = PETSC_NULL; > > Any thoughts on what might be going on? > > Thanks, > Amneet >

1 0

should this be moved now?
by Barry Smith 28 Feb '13

28 Feb '13

With the introduction of PetscPowRealInt() should this be moved to petscmath.h and renamed ? /* C++ does not promote int64_t to scalar or int32_t for std::pow() */ static PetscScalar Pow(PetscScalar b,PetscInt p) { return PetscPowScalar(b,(int)p); }

1 0

Re: [petsc-dev] [petsc-users] MATNEST with shell matrices
by Amneet Bhalla 28 Feb '13

28 Feb '13

Hi Jed, I have registered MatShell that operates on VecNest (which is composed of shell vectors) with KSP. When I am solving sytem of equations using this KSP, I am getting some errors. The sequence of events from gdb are: #0 0x0000000000000000 in ?? () #1 0x00000000009856fd in VecGetArrayRead (x=0x3422a30, a=0x7fff8b8c6620) at /home/amneet/SOFTWARES/LIBRARIESIBAMR/PETSc-3.3-Lite/include/petsc-private/vecimpl.h:350 #2 0x000000000098f57f in VecGetArray_Nest (X=0x34291d0, x=0x7fff8b8c6720) at /home/amneet/SOFTWARES/LIBRARIESIBAMR/PETSc-3.3-Lite/src/vec/vec/impls/nest/vecnest.c:719 #3 0x0000000000a3f3df in VecGetArrayRead (x=0x34291d0, a=0x7fff8b8c6720) at /home/amneet/SOFTWARES/LIBRARIESIBAMR/PETSc-3.3-Lite/include/petsc-private/vecimpl.h:350 #4 0x0000000000a51c83 in MatMult (mat=0x32c34f0, x=0x336a1c0, y=0x34291d0) at /home/amneet/SOFTWARES/LIBRARIESIBAMR/PETSc-3.3-Lite/src/mat/interface/matrix.c:2092 #5 0x00000000012b2a78 in KSPInitialResidual (ksp=0x334bcb0, vsoln=0x336a1c0, vt1=0x34291d0, vt2=0x343c150, vres=0x344f220, vb=0x337d810) at /home/amneet/SOFTWARES/LIBRARIESIBAMR/PETSc-3.3-Lite/src/ksp/ksp/interface/itres.c:54 #6 0x0000000000da2eaa in KSPSolve_GMRES (ksp=0x334bcb0) at /home/amneet/SOFTWARES/LIBRARIESIBAMR/PETSc-3.3-Lite/src/ksp/ksp/impls/gmres/gmres.c:230 #7 0x0000000000d1ad6b in KSPSolve (ksp=0x334bcb0, b=0x337d810, x=0x336a1c0) at /home/amneet/SOFTWARES/LIBRARIESIBAMR/PETSc-3.3-Lite/src/ksp/ksp/interface/itfunc.c:446 #8 0x00000000006b8a08 in IBTK::PETScBlockKrylovLinearSolver::solveSystem ( this=0x3249380, x=0x32f3e20, b=0x330c6d0) at PETScBlockKrylovLinearSolver.C:384 #9 0x00000000005777e3 in main (argc=6, argv=0x7fff8b8ca988) at main.C:274 In frame 4, the user defined matrix-vector multiplication is successful ( I printed out a statement). This would imply in the routine MatVecMult(Mat A, Vec x, Vec y), Vec y would be a duplicated vector. This in turn would mean that user-defined routine to duplicate vectors MatGetVecs(Mat A, Vec *left, Vec* right) should have been called. But this is not happening. It is in this routine, that I am setting *left->ops->getarray = PETSC_NULL; Any thoughts on what might be going on? Thanks, Amneet

1 0

Re: [petsc-dev] Vector operations on ghosted vectors
by Barry Smith 28 Feb '13

28 Feb '13

On Feb 28, 2013, at 7:30 AM, John Mousel <elafint.john(a)gmail.com> wrote: > Let me rephrase. Is a functionality like this something that is worth adding to the PETSc source, or should I just write and use my own wrapper. Just write and use your own wrapper. Barry > > > On Wed, Feb 27, 2013 at 4:51 PM, Matthew Knepley <knepley(a)gmail.com> wrote: > On Wed, Feb 27, 2013 at 4:56 PM, John Mousel <john.mousel(a)gmail.com> wrote: > Is there a possibility of adding a wrapper function around a few basic vector operations such as VecCopy, VecAXPY, VECAXPYPZ... to operate on ghosted vectors? I perform a lot of vector operations including the ghost region to avoid communication. It really clutters code to see stuff like > > CALL VecGhostGetLocalForm(solver%u,uk,ierr) > CALL VecGhostGetLocalForm(solver%ukm1,ukm1,ierr) > CALL VecGhostGetLocalForm(solver%p,pk,ierr) > CALL VecGhostGetLocalForm(solver%pkm1,pkm1,ierr) > CALL VecCopy(uk,ukm1,ierr) > CALL VecCopy(pk,pkm1,ierr) > CALL VecGhostRestoreLocalForm(solver%u,uk,ierr) > CALL VecGhostRestoreLocalForm(solver%ukm1,ukm1,ierr) > CALL VecGhostRestoreLocalForm(solver%p,pk,ierr) > CALL VecGhostRestoreLocalForm(solver%pkm1,pkm1,ierr) > > Maybe there's a good reason for this not already existing, or maybe it does exist and I haven't found it, but I'd think something like VecCopyGhosted would be fairly useful. > > Is there something we could do that would be better than you defining a function VecCopyGhosted()? > > Thanks, > > Matt > > John > -- > What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead. > -- Norbert Wiener >

1 0

Re: [petsc-dev] OpenMP in PETSc when calling from Fortran?
by Barry Smith 28 Feb '13

28 Feb '13

1) Work with petsc-dev http://www.mcs.anl.gov/petsc/developers/index.html the threaded stuff is more mature. 2) You should not need petscthreadcomm.h in Fortran. Simply using OpenMP progmas in your portion of the code. 3) See http://www.mcs.anl.gov/petsc/features/threads.html for a little bit of detail and feel free to post additional questions on this mailing list. Barry On Feb 28, 2013, at 4:35 AM, Åsmund Ervik <asmund.ervik(a)ntnu.no> wrote: > Hi, > > Having read that petsc-dev has the possibility of using OpenMP to run PC and KSP routines in parallel, and that people are getting nice speedup with this, I wanted to make use of this feature in our code. It is an incompressible two-phase flow solver programmed in Fortran which used PETSc for solving the pressure Poisson equation. > > I am, however, unable to use PETSC with OpenMP from Fortran. In particular, I can't seem to find a Fortran header file corresponding to petscthreadcomm.h (which is for C/C++). Is this not implemented yet? If so, is only the header file missing, or is more work required? > > Regards, > Åsmund Ervik

1 0

Re: [petsc-dev] Creating a dmplex mesh from a cell list
by Garnet Vaz 28 Feb '13

28 Feb '13

Hi Matt, Yes. Thank you. On Thu, Feb 28, 2013 at 6:53 AM, Matthew Knepley <knepley(a)gmail.com> wrote: > On Thu, Feb 28, 2013 at 8:24 AM, Garnet Vaz <garnet.vaz(a)gmail.com> wrote: >> >> Hi, >> >> I would like to know what some of the arguments for >> DMPlexCreateFromCellList are supposed to be for a triangular mesh. >> >> The function is: >> DMPlexCreateFromCellList(MPI_Comm comm, PetscInt dim, PetscInt >> numCells, PetscInt numVertices, PetscInt numCorners, PetscBool >> interpolate, const int cells[], PetscInt spaceDim, const double >> vertexCoords[], DM *dm) >> >> dim = 2 >> numCells = nCells (not the entire length of the cells array; instead >> length/3) >> numVertices = nVert (length of vertexCoords/2) >> numCorners = 3 >> spaceDim = 2 >> >> What is wrong with my understanding about the following variables? >> The error I get is the following: >> >> [0]PETSC ERROR: Argument out of range! >> [0]PETSC ERROR: Cone point 43 is not in the valid range [0, 43)! > > > I have created a proper webpage for this function which will be on the web > at midnight. Here it is > > /*@C > DMPlexCreateFromCellList - This takes as input common mesh generator > output, a list of the vertices for each cell, and produces a DM > > Input Parameters: > + comm - The communicator > . dim - The topological dimension of the mesh > . numCells - The number of cells > . numVertices - The number of vertices > . numCorners - The number of vertices for each cell > . interpolate - Flag indicating that intermediate mesh entities (faces, > edges) should be created automatically > . cells - An array of numCells*numCorners numbers, the vertices for each > cell > . spaceDim - The spatial dimension used for coordinates > - vertexCoords - An array of numVertices*spaceDim numbers, the coordinates > of each vertex > > Output Parameter: > . dm - The DM > > Note: Two triangles sharing a face > $ > $ 2 > $ / | \ > $ / | \ > $ / | \ > $ 0 0 | 1 3 > $ \ | / > $ \ | / > $ \ | / > $ 1 > would have input > $ numCells = 2, numVertices = 4 > $ cells = [0 1 2 1 3 2] > $ > which would result in the DMPlex > $ > $ 4 > $ / | \ > $ / | \ > $ / | \ > $ 2 0 | 1 5 > $ \ | / > $ \ | / > $ \ | / > $ 3 > > Level: beginner > > .seealso: DMPlexCreate() > @*/ > > Does this help? > > Thanks, > > Matt > >> >> >> Thank You. >> >> -- >> Regards, >> Garnet > > > > > -- > What most experimenters take for granted before they begin their experiments > is infinitely more interesting than any results to which their experiments > lead. > -- Norbert Wiener -- Regards, Garnet

1 0

Re: [petsc-dev] Creating a dmplex mesh from a cell list
by Matthew Knepley 28 Feb '13

28 Feb '13

On Thu, Feb 28, 2013 at 8:24 AM, Garnet Vaz <garnet.vaz(a)gmail.com> wrote: > Hi, > > I would like to know what some of the arguments for > DMPlexCreateFromCellList are supposed to be for a triangular mesh. > > The function is: > DMPlexCreateFromCellList(MPI_Comm comm, PetscInt dim, PetscInt > numCells, PetscInt numVertices, PetscInt numCorners, PetscBool > interpolate, const int cells[], PetscInt spaceDim, const double > vertexCoords[], DM *dm) > > dim = 2 > numCells = nCells (not the entire length of the cells array; instead > length/3) > numVertices = nVert (length of vertexCoords/2) > numCorners = 3 > spaceDim = 2 > > What is wrong with my understanding about the following variables? > The error I get is the following: > > [0]PETSC ERROR: Argument out of range! > [0]PETSC ERROR: Cone point 43 is not in the valid range [0, 43)! > I have created a proper webpage for this function which will be on the web at midnight. Here it is /*@C DMPlexCreateFromCellList - This takes as input common mesh generator output, a list of the vertices for each cell, and produces a DM Input Parameters: + comm - The communicator . dim - The topological dimension of the mesh . numCells - The number of cells . numVertices - The number of vertices . numCorners - The number of vertices for each cell . interpolate - Flag indicating that intermediate mesh entities (faces, edges) should be created automatically . cells - An array of numCells*numCorners numbers, the vertices for each cell . spaceDim - The spatial dimension used for coordinates - vertexCoords - An array of numVertices*spaceDim numbers, the coordinates of each vertex Output Parameter: . dm - The DM Note: Two triangles sharing a face $ $ 2 $ / | \ $ / | \ $ / | \ $ 0 0 | 1 3 $ \ | / $ \ | / $ \ | / $ 1 would have input $ numCells = 2, numVertices = 4 $ cells = [0 1 2 1 3 2] $ which would result in the DMPlex $ $ 4 $ / | \ $ / | \ $ / | \ $ 2 0 | 1 5 $ \ | / $ \ | / $ \ | / $ 3 Level: beginner .seealso: DMPlexCreate() @*/ Does this help? Thanks, Matt > > Thank You. > > -- > Regards, > Garnet > -- What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead. -- Norbert Wiener

1 0

Re: [petsc-dev] Vector operations on ghosted vectors
by John Mousel 28 Feb '13

28 Feb '13

Let me rephrase. Is a functionality like this something that is worth adding to the PETSc source, or should I just write and use my own wrapper. On Wed, Feb 27, 2013 at 4:51 PM, Matthew Knepley <knepley(a)gmail.com> wrote: > On Wed, Feb 27, 2013 at 4:56 PM, John Mousel <john.mousel(a)gmail.com>wrote: > >> Is there a possibility of adding a wrapper function around a few basic >> vector operations such as VecCopy, VecAXPY, VECAXPYPZ... to operate on >> ghosted vectors? I perform a lot of vector operations including the ghost >> region to avoid communication. It really clutters code to see stuff like >> >> CALL VecGhostGetLocalForm(solver%u,uk,ierr) >> CALL VecGhostGetLocalForm(solver%ukm1,ukm1,ierr) >> CALL VecGhostGetLocalForm(solver%p,pk,ierr) >> CALL VecGhostGetLocalForm(solver%pkm1,pkm1,ierr) >> CALL VecCopy(uk,ukm1,ierr) >> CALL VecCopy(pk,pkm1,ierr) >> CALL VecGhostRestoreLocalForm(solver%u,uk,ierr) >> CALL VecGhostRestoreLocalForm(solver%ukm1,ukm1,ierr) >> CALL VecGhostRestoreLocalForm(solver%p,pk,ierr) >> CALL VecGhostRestoreLocalForm(solver%pkm1,pkm1,ierr) >> >> Maybe there's a good reason for this not already existing, or maybe it >> does exist and I haven't found it, but I'd think something like >> VecCopyGhosted would be fairly useful. >> > > Is there something we could do that would be better than you defining a > function VecCopyGhosted()? > > Thanks, > > Matt > > >> John >> > -- > What most experimenters take for granted before they begin their > experiments is infinitely more interesting than any results to which their > experiments lead. > -- Norbert Wiener >

1 0

Creating a dmplex mesh from a cell list
by Garnet Vaz 28 Feb '13

28 Feb '13

Hi, I would like to know what some of the arguments for DMPlexCreateFromCellList are supposed to be for a triangular mesh. The function is: DMPlexCreateFromCellList(MPI_Comm comm, PetscInt dim, PetscInt numCells, PetscInt numVertices, PetscInt numCorners, PetscBool interpolate, const int cells[], PetscInt spaceDim, const double vertexCoords[], DM *dm) dim = 2 numCells = nCells (not the entire length of the cells array; instead length/3) numVertices = nVert (length of vertexCoords/2) numCorners = 3 spaceDim = 2 What is wrong with my understanding about the following variables? The error I get is the following: [0]PETSC ERROR: Argument out of range! [0]PETSC ERROR: Cone point 43 is not in the valid range [0, 43)! Thank You. -- Regards, Garnet

1 0