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

MatWrapML_SHELL
by John Fettig 09 May '12

09 May '12

In the function MatWrapML_SHELL, there is a return without creating the shell matrix. I have two questions about it: m = mlmat->outvec_leng; n = mlmat->invec_leng; if (!m || !n){ newmat = PETSC_NULL; PetscFunctionReturn(0); } 1) Should there be an error generated here? 2) If not, why does the function return without generating the matrix? This can result in deadlock if some processes return in this clause and some do not. The ones that do not get caught in the MatCreateShell() call. Regards, John

1 0

Re: [petsc-dev] API changes in MatIS
by Jed Brown 09 May '12

09 May '12

On Wed, May 9, 2012 at 1:29 PM, Stefano Zampini <stefano.zampini(a)gmail.com>wrote: > Finally I have some time to work on BDDC. I'm thinking to restyle some of > my BDDC code. In particular in setting up the coarse environment, I want > to implement a function (mostly reusing some BDDC code I already wrote) > > MatiSSubassemble(A, IS ISSubdomains,*B,*intercomm,*partial_scatter) > > which, given subdomain indices, creates a new MATIS matrix on a subcomm of > the communicator of A by subassembling. In particular, ISSubdomains should > come from a call to MatPartitioningApply on the adjacency matrix of the > subdomains. > > The communicator of B is the subcomm of A of processes which will receive > some of the local matrices > intercomm will be the communicator which can be used between each > receiving process and their sending "friends" > *partial_scatter will be the scatter context associated to the intercomm > > Jed, do you already wrote something similar for your PA? What's the best > logic for you? Should I use MatCoarsen? > MatCoarsen is for graph aggregation/matching, which is different. It seems to me that subassembly doesn't coarsen at all. If MatGetSubMatrices() could extract parallel matrices directly, then it looks more like that's what you're after. I don't think the implementation should use an intercomm (too many implementations are buggy and it requires MPI-2 which PETSc tries to avoid depending on because enough users still have antique systems that don't support it). In any case, I think that creating a partition should be separate from transforming the matrix. For MatIS, I would be fine with creating a new MatIS using the same local matrices on a subcommunicator, then MatConvert'ing it to AIJ (that will do assembly), and discarding the temporary subcomm MatIS.

1 0

Re: [petsc-dev] API changes in MatIS
by Stefano Zampini 09 May '12

09 May '12

Finally I have some time to work on BDDC. I'm thinking to restyle some of my BDDC code. In particular in setting up the coarse environment, I want to implement a function (mostly reusing some BDDC code I already wrote) MatiSSubassemble(A, IS ISSubdomains,*B,*intercomm,*partial_scatter) which, given subdomain indices, creates a new MATIS matrix on a subcomm of the communicator of A by subassembling. In particular, ISSubdomains should come from a call to MatPartitioningApply on the adjacency matrix of the subdomains. The communicator of B is the subcomm of A of processes which will receive some of the local matrices intercomm will be the communicator which can be used between each receiving process and their sending "friends" *partial_scatter will be the scatter context associated to the intercomm Jed, do you already wrote something similar for your PA? What's the best logic for you? Should I use MatCoarsen?

1 0

Re: [petsc-dev] OpenMP Support
by Gerard Gorman 09 May '12

09 May '12

Hi Dave OpenMP support exists for vec and mat (only AIJ so far). There is a big difference in performance depending on available memory bandwidth and the compiler OpenMP implementation. In application codes (such as Fluidity which is our main target code for this work) there are other significant costs such as matrix assembly. So in general you have to consider how easy it will be to thread to the other computationally expensive sections of your code, otherwise the overall speed-up of your application will be modest. Cheers Gerard Dave Nystrom emailed the following on 09/05/12 04:29: > Is the pthreads support further along than the OpenMP support? I have not > tried the pthreads support yet. Does either the pthreads support or the > OpenMP support implement the matvec or do they just do vector type > operations? > > Jed Brown writes: > > On Tue, May 8, 2012 at 9:23 PM, Dave Nystrom <dnystrom1(a)comcast.net> wrote: > > > > > I see that petsc-dev now has some OpenMP support. Would a serial, non-mpi > > > code that uses petsc-dev be able to gain much performance improvement from > > > it > > > now for the case of doing sparse linear solve with cg and jacobi > > > preconditioning? > > > > > > > The kernels are being transitioned to use the threadcomm, which enables > > OpenMP and other threading models. > > > > We anticipate that pthreads will provide the best performance because > > operations can be less synchronous than with OpenMP (for which a parallel > > region implies barrier semantics). But if other parts of an application > > are using OpenMP, it would be preferable for PETSc to also use OpenMP so > > that it can share the same thread pool. The same applies to TBB.

1 0

Re: [petsc-dev] OpenMP Support
by Jed Brown 08 May '12

08 May '12

On Tue, May 8, 2012 at 10:29 PM, Dave Nystrom <Dave.Nystrom(a)tachyonlogic.com > wrote: > Is the pthreads support further along than the OpenMP support? I have not > tried the pthreads support yet. Does either the pthreads support or the > OpenMP support implement the matvec or do they just do vector type > operations? > The pthreads stuff has been available for a while, including some matrix kernels, but it's being transformed now, so it will get better.

1 0

Re: [petsc-dev] OpenMP Support
by Dave Nystrom 08 May '12

08 May '12

Is the pthreads support further along than the OpenMP support? I have not tried the pthreads support yet. Does either the pthreads support or the OpenMP support implement the matvec or do they just do vector type operations? Jed Brown writes: > On Tue, May 8, 2012 at 9:23 PM, Dave Nystrom <dnystrom1(a)comcast.net> wrote: > > > I see that petsc-dev now has some OpenMP support. Would a serial, non-mpi > > code that uses petsc-dev be able to gain much performance improvement from > > it > > now for the case of doing sparse linear solve with cg and jacobi > > preconditioning? > > > > The kernels are being transitioned to use the threadcomm, which enables > OpenMP and other threading models. > > We anticipate that pthreads will provide the best performance because > operations can be less synchronous than with OpenMP (for which a parallel > region implies barrier semantics). But if other parts of an application > are using OpenMP, it would be preferable for PETSc to also use OpenMP so > that it can share the same thread pool. The same applies to TBB.

1 0

Re: [petsc-dev] OpenMP Support
by Jed Brown 08 May '12

08 May '12

On Tue, May 8, 2012 at 9:23 PM, Dave Nystrom <dnystrom1(a)comcast.net> wrote: > I see that petsc-dev now has some OpenMP support. Would a serial, non-mpi > code that uses petsc-dev be able to gain much performance improvement from > it > now for the case of doing sparse linear solve with cg and jacobi > preconditioning? > The kernels are being transitioned to use the threadcomm, which enables OpenMP and other threading models. We anticipate that pthreads will provide the best performance because operations can be less synchronous than with OpenMP (for which a parallel region implies barrier semantics). But if other parts of an application are using OpenMP, it would be preferable for PETSc to also use OpenMP so that it can share the same thread pool. The same applies to TBB.

1 0

OpenMP Support
by Dave Nystrom 08 May '12

08 May '12

I see that petsc-dev now has some OpenMP support. Would a serial, non-mpi code that uses petsc-dev be able to gain much performance improvement from it now for the case of doing sparse linear solve with cg and jacobi preconditioning? Thanks, Dave

1 0

Re: [petsc-dev] API changes in MatIS
by Jed Brown 08 May '12

08 May '12

On Tue, May 8, 2012 at 11:05 AM, Stefano Zampini <stefano.zampini(a)gmail.com>wrote: > Not since now. I think it is were such an approach belongs, apart from the > fact that PCREDISTRIBUTE is not so visible from a standard user's > perspective (as me). I think I can put my code in PCSetUp_REDISTRIBUTE, > Adding . If I understood pcredistribute, I should do from the higher level > something like > > MatCreateIS(A) > KSPCreate(true_ksp) > KSPSetOperators(true_ksp,A,A) > KSPSetType(true_ksp,KSPPREONLY) > KSPGetPC(true_ksp,red_pc) > PCSetType(red_pc,PC_REDISTRIBUTE) > KSPSetup(true_ksp) > This KSPSetUp() need not be called explicitly. > > KSPSolve(true_ksp,b,x) > > I didn't understand the "diagonal portion" stuff. > It's common for people to just "zero" rows that correspond to Dirichlet boundary conditions. This eliminates those trivial equations. It is the place that I think we should put other redistribution (we can have variants that are more or less aggressive). > Load balancing is assured if you come from a Mesh partitioner which does > it for you. For communications costs, optimizing the global ordering of a > MATIS object is not so simple. In principle you have more local dofs than > that stored in the local part of the global vector you own. With my code > (I've tested it on structured grids) I've always observed more than 80% of > the whole number of global dofs assigned to the process to which they > belong in the local representation of MATIS (mostly INTERNAL nodes). As an > example, a usual Natural ordering for structured grids counts to about 1%. > I thought a fairly natural thing is to partition elements and then distribute vertex ownership to the owner of an element that touches, for approximate vertex balance. Multi-objective partitioning is also possible. I just think that for good assembly efficiency, this is something that should be done earlier. It's fine to also support doing it later, but the semantics have to be consistent. > I didn't used DM so far. How can you set up a DM for a MATIS object? > DM manages geometric and problem-specific information. The most commonly used implementation of DM is DMDA which does logically structured grids. To support MATIS, you would look at DMCreateMatrix_DA(). > So, for the smoother something like the hybrid approach proposed here > http://www.cs.nyu.edu/cs/faculty/widlund/ijnme_submit_revision_ijnme.pdf ? It could use that, but I think/hope that less will be required since the only thing lacking is accounting for error in the multilevel decomposition, not the spatial decomposition.

1 0

Re: [petsc-dev] API changes in MatIS
by Stefano Zampini 08 May '12

08 May '12

Asking the first question: >> The function works by first reordering local nodes with INTERNAL and >> INTERFACE nodes (if the second argument LocalIS is not PETSC_NULL), in >> order to be able to use memcpy when using local scatters of PCIS (e.g. >> pcis->N_to_B). Then, it reorders the global numbering (if GlobalIS is not >> PETSC_NULL) by exploiting data locality in the global vector and reduce the >> amount of communications when scattering from global to local and >> viceversa; in practice, a numbering of global nodes is constructed in such >> a way that the greatest number of subdomain INTERNAL nodes will belong to >> the local part of the global vector. Then, remaining INTERNAL nodes (notes >> that if the parallel problem is well balanced, this will not occur) and >> INTERFACE nodes are placed in the remaining global nodes, keeping them >> sorted as they were in the previous global numbering. I think the latter >> part can be further modified. When I'll add to matis.c, maybe you can take >> a look at the code and give me your suggestions. >> > > Okay, this is different from what I thought you were describing. Moving > the global dofs is a tricky thing because the vectors we want to solve > against and interpret results for come from a higher level. Changing the > global ordering really produces a different matrix (MatMult is different), > so it's not "optimization". If the linear operator stayed the same, but you > just changed the format, perhaps we could call it "optimize". > > Are you aware of PCRedistribute? It could be made to perform the > transformation you are thinking about for MATIS, solving the problem in the > "optimized" ordering, and mapping the solution back to the original > ordering. > > Not since now. I think it is were such an approach belongs, apart from the fact that PCREDISTRIBUTE is not so visible from a standard user's perspective (as me). I think I can put my code in PCSetUp_REDISTRIBUTE, Adding . If I understood pcredistribute, I should do from the higher level something like MatCreateIS(A) KSPCreate(true_ksp) KSPSetOperators(true_ksp,A,A) KSPSetType(true_ksp,KSPPREONLY) KSPGetPC(true_ksp,red_pc) PCSetType(red_pc,PC_REDISTRIBUTE) KSPSetup(true_ksp) KSPSolve(true_ksp,b,x) I didn't understand the "diagonal portion" stuff. > But why is this needed in the first place? Don't people who use MATIS also > create a reasonable global ordering at setup time? It's bad for other parts > of the algorithm to have a global numbering that doesn't line up well with > the partition of elements. > > Load balancing is assured if you come from a Mesh partitioner which does it for you. For communications costs, optimizing the global ordering of a MATIS object is not so simple. In principle you have more local dofs than that stored in the local part of the global vector you own. With my code (I've tested it on structured grids) I've always observed more than 80% of the whole number of global dofs assigned to the process to which they belong in the local representation of MATIS (mostly INTERNAL nodes). As an example, a usual Natural ordering for structured grids counts to about 1%. > >> Regarding the other question: >> With MATIS (and in principle every subassembled matrix) you are forced to >> use the local matrices assembled by the user or an assembled group of them >> (as you surely know since you are developing a PA). You cannot safely (in a >> theoretical way) move subdomain boundaries since you don't know (at least >> in principle) how to reproduce the matrix values on the local interface >> since they are neumann boundaries for the local PDE. The main question is: >> how can I split the assembled values on the nodes which previously were >> internal to the subdomain? >> > > This can't be done, at least without more information or making > assumptions about the discretization. That's why I asked. > > >> (I think this is the main question when thinking to a MatConvert >> implementation which converts from a matrix obtained by DMDAGetMatrix to a >> MATIS). >> > > There is no reason why DMCreateMatrix() cannot return a MATIS. In any > case, that's a better model than assembling in one format and trying to > convert, especially since conversion really can't done in this direction (a > MATIS can be converted to MATAIJ, not the other way). > > I didn't used DM so far. How can you set up a DM for a MATIS object? > >> >>> >>> For a different, but related algorithm, I'm considering writing a MatPA >>> (PA = partially assembled) which would choose some points to assemble >>> internally. The reason for storing in partially assembled form is that I >>> will want to solve a number of pinned Neumann problems to determine the >>> coarse space and I'd rather not have to take unnecessary submatrices. Do >>> you think this is a bad idea and I should just use MatIS? (I also want >>> support for multiple subdomains per process and multiple processes per >>> subdomain, though that could be added to MatIS.) >>> >> >> Cool. Are you implementing FETI-DP? I'm planning to implement it in the >> very next future for isogeometric analysis. Asking to your question, it >> depends on how you need to solve these "pinned Neumann problems": you need >> to solve unrelated local problems defined by the local matrices or solve a >> parallel Neumann problem? In the first case, MatIS is your choice. In the >> other case, I think you need to resort to a multigrid for a new MATIS >> object defined by the subdomains of interest. >> However, I will be glad to help you with PA if you want. > > > I'm not that interested in FETI-DP since it's basically the same as BDDC > and doesn't fix the exponential dependence on the number of levels. Also, I > hate that it doesn't work well for irregular coefficients on subdomain > boundaries. So I have a different algorithm in mind, more naturally viewed > as a multigrid or multilevel Schwarz, but with aggressive coarsening > enabled by the subdomains. One part of the method is compatible relaxation, > for which pinned Neumann problems are used. The idea is to start with just > enough pinning to make the subdomains non-singular, then enrich the coarse > space by adaptively selecting vertices or edge/face moments that are not > well controlled (as determined by compatible relaxation). Once the coarse > injection is found, we construct the coarse basis functions by harmonic > extension (again using pinned Neumann problems). With the multigrid > hierarchy is constructed, the smoother would be based on overlapping > subdomains or a multi-stage smoother involving both pinned Neumann problems > and surface strips (or perhaps just point-block Jacobi). The problem with > only using the pinned Neumann problems in the smoother seems to be that we > can't remove the exponential dependence without allowing the smoother to > move the injection (compatible relaxation is good for telling us whether a > coarse space is good, but not actually the right thing for multigrid > convergence). > So, for the smoother something like the hybrid approach proposed here http://www.cs.nyu.edu/cs/faculty/widlund/ijnme_submit_revision_ijnme.pdf ? -- Stefano

1 0