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

Re: [petsc-dev] the complex dot business
by Matthew Knepley 27 Feb '12

27 Feb '12

On Mon, Feb 27, 2012 at 10:55 AM, Barry Smith <bsmith(a)mcs.anl.gov> wrote: > > I'm afraid the sudden appearance of the blas dot for complex numbers > problem is my fault. > > Traditionally since 199x we never used the complex dot product from blas > because of the return complex type problem. VecDot_Seq() used #ifdef for > complex and a simple loop. > > Since some blas norm suck I off-hand told Shri to use dot for norm and > wasn't thinking the complex case. > > Since then, of course, the complex dot problem has been coming up. > > Maybe it is best to continue to never use complex blas dot and just > change the VecNorm_ routines to use norm for complex and not dot.? Just > throw away this silly problem and time sink? Why is this such a problem? Just use PetscScalar for the BLASdot_() return type and its fine. Matt > > Barry > > -- 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] the complex dot business
by Jed Brown 27 Feb '12

27 Feb '12

On Mon, Feb 27, 2012 at 10:55, Barry Smith <bsmith(a)mcs.anl.gov> wrote: > > I'm afraid the sudden appearance of the blas dot for complex numbers > problem is my fault. > > Traditionally since 199x we never used the complex dot product from blas > because of the return complex type problem. VecDot_Seq() used #ifdef for > complex and a simple loop. > > Since some blas norm suck I off-hand told Shri to use dot for norm and > wasn't thinking the complex case. > Most BLAS implementations don't have great code for BLAS level 1 anyway and it tends to be memory bandwidth limited, so I don't have a problem with using a simple loop. > > Since then, of course, the complex dot problem has been coming up. > > Maybe it is best to continue to never use complex blas dot and just > change the VecNorm_ routines to use norm for complex and not dot.? Just > throw away this silly problem and time sink? I think we should use norm for norm.

1 0

Re: [petsc-dev] SNES, Coloring, DM, and TS
by Matthew Knepley 27 Feb '12

27 Feb '12

On Mon, Feb 27, 2012 at 10:48 AM, Jed Brown <jedbrown(a)mcs.anl.gov> wrote: > I'm getting lost with too many ways to do the same thing. The DM has only > one slot for a residual function, I'm not wild about having TS overwrite > that. So TS uses SNESSetFunction(), but now I need it on levels and I want > coloring to work. But the coloring setup code is inside DMComputeJacobian() > and I'd rather not duplicate it in SNES > > if (snes->ops->computefunction && !snes->ops->computejacobian && snes->dm > && !snes->dm->ops->function) { > /* set up coloring using DM, but apply using snes->ops->computefunction > */ > } > > or some such nastiness. > > I would like to either _always_ use DM (possibly created implicitly by > TS/SNES/KSP and just used to hold function pointers) or _never_ use DM. > > The problem with never using DM is that we need a way for the user to > decide to implement different physics on coarse levels. They currently do > this by implementing DMCoarsen() or using DMCoarsenHookAdd(). If the > function pointers are stored in SNES and TS, the hook would have to be able > to access them somehow. > > The problem with always using DM is that we effectively create a > dependency loop because the user is telling DM about a nonlinear or > transient problem even though that concept is not part of the DM interface. > I don't see why this is a dependency loop. > As a related matter, I'm stuck with a reference counting loop where I want > to attach special Vecs to a DM. I can wrap them in a PetscContainer to > circumvent the code that checks for reference cycles, but that makes > freeing them delicate. We could make a named cache similar to DM work > vectors, but where restored Vecs are not available for reuse. > I know that rationale for having Vecs know about a DM (output), which I don't like. What is the rationale for having DMs know about a certain Vec? Matt -- 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

the complex dot business
by Barry Smith 27 Feb '12

27 Feb '12

I'm afraid the sudden appearance of the blas dot for complex numbers problem is my fault. Traditionally since 199x we never used the complex dot product from blas because of the return complex type problem. VecDot_Seq() used #ifdef for complex and a simple loop. Since some blas norm suck I off-hand told Shri to use dot for norm and wasn't thinking the complex case. Since then, of course, the complex dot problem has been coming up. Maybe it is best to continue to never use complex blas dot and just change the VecNorm_ routines to use norm for complex and not dot.? Just throw away this silly problem and time sink? Barry

1 0

Re: [petsc-dev] SNES, Coloring, DM, and TS
by Dmitry Karpeev 27 Feb '12

27 Feb '12

On Mon, Feb 27, 2012 at 10:48 AM, Jed Brown <jedbrown(a)mcs.anl.gov> wrote: > I'm getting lost with too many ways to do the same thing. The DM has only > one slot for a residual function, I'm not wild about having TS overwrite > that. So TS uses SNESSetFunction(), but now I need it on levels and I want > coloring to work. But the coloring setup code is inside DMComputeJacobian() > and I'd rather not duplicate it in SNES > > if (snes->ops->computefunction && !snes->ops->computejacobian && snes->dm > && !snes->dm->ops->function) { > /* set up coloring using DM, but apply using snes->ops->computefunction > */ > } > > or some such nastiness. > > I would like to either _always_ use DM (possibly created implicitly by > TS/SNES/KSP and just used to hold function pointers) or _never_ use DM. > > The problem with never using DM is that we need a way for the user to > decide to implement different physics on coarse levels. They currently do > this by implementing DMCoarsen() or using DMCoarsenHookAdd(). If the > function pointers are stored in SNES and TS, the hook would have to be able > to access them somehow. > > The problem with always using DM is that we effectively create a > dependency loop because the user is telling DM about a nonlinear or > transient problem even though that concept is not part of the DM interface. > Why is a nonlinear problem not part of the DM interface? In my mind a DM encapsulates a problem, that also knows how to linearize/split/coarsen/refine itself by returning the corresponding subDMs. Dmitry. > > > As a related matter, I'm stuck with a reference counting loop where I want > to attach special Vecs to a DM. I can wrap them in a PetscContainer to > circumvent the code that checks for reference cycles, but that makes > freeing them delicate. We could make a named cache similar to DM work > vectors, but where restored Vecs are not available for reuse. >

1 0

SNES, Coloring, DM, and TS
by Jed Brown 27 Feb '12

27 Feb '12

I'm getting lost with too many ways to do the same thing. The DM has only one slot for a residual function, I'm not wild about having TS overwrite that. So TS uses SNESSetFunction(), but now I need it on levels and I want coloring to work. But the coloring setup code is inside DMComputeJacobian() and I'd rather not duplicate it in SNES if (snes->ops->computefunction && !snes->ops->computejacobian && snes->dm && !snes->dm->ops->function) { /* set up coloring using DM, but apply using snes->ops->computefunction */ } or some such nastiness. I would like to either _always_ use DM (possibly created implicitly by TS/SNES/KSP and just used to hold function pointers) or _never_ use DM. The problem with never using DM is that we need a way for the user to decide to implement different physics on coarse levels. They currently do this by implementing DMCoarsen() or using DMCoarsenHookAdd(). If the function pointers are stored in SNES and TS, the hook would have to be able to access them somehow. The problem with always using DM is that we effectively create a dependency loop because the user is telling DM about a nonlinear or transient problem even though that concept is not part of the DM interface. As a related matter, I'm stuck with a reference counting loop where I want to attach special Vecs to a DM. I can wrap them in a PetscContainer to circumvent the code that checks for reference cycles, but that makes freeing them delicate. We could make a named cache similar to DM work vectors, but where restored Vecs are not available for reuse.

1 0

Re: [petsc-dev] OpenMP/Vec
by Jed Brown 27 Feb '12

27 Feb '12

On Sun, Feb 26, 2012 at 22:54, recrusader <recrusader(a)gmail.com> wrote: > For multithread support in PETSc, my question is whether KSP and/or PC > work when Vec and Mat use multithread mode. > Yes

1 0

Re: [petsc-dev] OpenMP/Vec
by recrusader 26 Feb '12

26 Feb '12

Dear Jed, For multithread support in PETSc, my question is whether KSP and/or PC work when Vec and Mat use multithread mode. Thanks, Yujie On Sun, Feb 26, 2012 at 6:39 PM, Jed Brown <jedbrown(a)mcs.anl.gov> wrote: > On Sun, Feb 26, 2012 at 04:07, Gerard Gorman <g.gorman(a)imperial.ac.uk>wrote: > >> > Did you post a repository yet? I'd like to have a look at the code. >> >> It's on Barry's favourite collaborative software development site of >> course ;-) >> >> https://bitbucket.org/wence/petsc-dev-omp/overview >> > > I looked through the code and I'm concerned that all the OpenMP code is > inlined into vec/impls/seq and mat/impls/aij/seq with, as far as I can > tell, no way to use OpenMP for some objects in a simulation, but not > others. I think that all the pragmas should have num_threads(x->nthreads) > clauses. We can compute the correct number of threads based on sizes when > memory is allocated (or specified through command line options, inherited > from related objects, etc). > > I don't think we can get away from OpenMP schedule overhead (several > hundred cycles) even for those objects that we choose not to use threads > for, but (at least with my gcc tests), that overhead is only about a third > of the cost of actually starting a parallel region. > > It's really not acceptable to insert unguarded > > #pragma omp ... > > into the code because this will generate tons of warnings or errors with > compilers that don't know about OpenMP. It would be better to test for > _Pragma and use > > #define PetscPragmatize(x) _Pragma(#x) > #if defined(PETSC_USE_OPENMP) > # define PetscPragmaOMP(x) PetscPragmatize(omp x) > #else > # define PetscPragmaOMP(x) > #endif > > then use > > PetscPragmaOMP(parallel for ...) > > We should probably use a variant for object-based threading > > #define PetscPragmaOMPObject(obj,x) PetscPragmaOMP(x > num_threads((obj)->nthreads)) > > In the case of multiple objects, I think you usually want the object being > modified to control the number of threads. > > In many cases, I would prefer more control over the partition of the loop. > For example, in many cases, I'd be willing to tolerate a slight > computational imbalance between threads in exchange for working exclusively > within my page. Note that the arithmetic to compute such things is orders > of magnitude less expensive than the schedule/distribution to threads. I > don't know how to do that except to > > PragmaOMP(parallel) { > int nthreads = omp_get_num_threads(); > int tnum = omp_get_thread_num(); > int start,end; > // compute start and end > for (int i=start; i<end; i++) { > // the work > } > } > > We could perhaps capture some of this common logic in a macro: > > #define VecOMPParallelBegin(X,args) do { \ > PragmaOMPObject(X,parallel args) { \ > PetscInt _start, _end; \ > VecOMPGetThreadLocalPart(X,&_start,&_end); \ > { do {} while(0) > > #define VecOMPParallelEnd() }}} while(0) > > VecOMPParallelBegin(X, shared/private ...); > { > PetscInt i; > for (i=_start; i<_end; i++) { > // the work > } > } > VecOMPParallelEnd(); > > That should reasonably give us complete run-time control of the number of > parallel threads per object and their distribution, within the constraints > of contiguous thread partition. That also leaves open the possibility of > using libnuma to query and migrate pages. (For example, a short vector that > needs to be accessed from multiple NUMA nodes might intentionally be > faulted with pages spread apart even though other vectors of similar size > might be accessed from within one NUMA nodes and thus not use threads at > all. (One 4 KiB page is only 512 doubles, but if the memory is local to a > single NUMA node, we wouldn't use threads until the vector length was 4 to > 8 times larger.) >

1 0

Re: [petsc-dev] DM field interfaces
by Jed Brown 26 Feb '12

26 Feb '12

On Sun, Feb 26, 2012 at 21:56, Barry Smith <bsmith(a)mcs.anl.gov> wrote: > Presumably a "new" DM that represents the new set of equations ... and > to start we will need to do very specific examples before coming up with a > general paradigm. > > But it seems best to me to make this orthogonal to PCFieldSplit. > Okay, sounds good to me.

1 0

Re: [petsc-dev] DM field interfaces
by Barry Smith 26 Feb '12

26 Feb '12

See bottom. On Feb 26, 2012, at 9:33 PM, Jed Brown wrote: > On Sat, Feb 25, 2012 at 15:27, Barry Smith <bsmith(a)mcs.anl.gov> wrote: > > Let's back up a step before everyone gets totally confused (and worry less about exact syntax to use instead focus on concepts) > > Proposal 1: Jed > > // Conservative gas dynamics with fields ordered as [rho, rho*u, rho*v, rho*w, E] > DMGetFieldSplitting(dm,"conservative-vector",&fscons); // creates if it doesn't exist > DMFieldSplittingAddField(dm,fscons,"rho",{0}); > DMFieldSplittingAddField(dm,fscons,"momentum",{1,2,3}); > > This seems to be premised on having an underlying "canonical" set of fields that one can define unions of to make new fields. For example with DMDA the underlying > fields are simply the 0th, 1st, 2nd etc DOF at each node. For a simple staggered grid 0 canonical may be the pressure centers, 1 the x velocities, 2 the y velocities > > Alright, lets abandon this "canonical field" idea, it was broken to begin with and less flexible, as you say. > > > Proposal 2: Barry > > Use IS's to define the fields, not unions of canonical fields. Reason: more general than proposal 1 > > Great, my concern is that we would like DM to be able to declare multiple named partitions. For example, so that we can have a conservative-scalar splitting, a conservative vector splitting (keeps vector components together), primitive variable splittings, local-characteristic splittings (analytically diagonalize the flux Jacobian for hyperbolic problems), etc. > > We can think of these transformations as symmetric preconditioning. Suppose that we want to solve with J, operating on conservative variables. Suppose P converts conservative to primitive variables and C = P^{-1} converts primitive to conservative. I'm using matrix notation for because this would often be applied to perturbations. We need to solve a linear system with J_c, the Jacobian acting on conservative variables, but the preconditioner wants to operate in primitive variables, so the user assembles > > J_p = P J_c C. > > (Usually they do this only approximately and it's done without looking at J_c.) The conservative system is > > J_c = C J_p P > > At the algebraic level, I think this logic belongs in a new preconditioner that applies a user-defined pre- and post-transformation, the concern is how to forward DM-based information into PCFieldSplit and PCMG operating on J_p. > > > Because Jed's proposal 1 violates the use of IS for indexing paradigm and is less general than Proposal 2 I'd like to drop that Proposal 1; but it does bring up an issue. I think we need another IS implementation that is more "general" than stride for handling sets of DOFs of freedom associated with those canonical sets of fields. So it is easy to access canonical fields that are not simply stride. For example ISCreateCanonical(MPI_Comm,PetscInt ncanonical,PetscInt *fields), we can even have "predefined" ones available for use such as IS_CANONICAL_(MPI_Comm,3,{1,2,3}) Thus we can preserve some of the simplicity of proposal 1. Jed's (again not worrying about syntax) > > DMFieldSplittingAddField(dm,fscons,"momentum",{1,2,3}); becomes > > DMFieldSplittingAddField(dm,fscons,"momentum",IS_CANONICAL_(comm,3,{1,2,3}) > > Note that for a Vec with a given block size the canonical fields are simply the stride fields so IS_CANONICAL_(comm,1,{start}) is pretty much ISStride(). > > Ignoring syntax, I think this is fine. > > Proposal 3: Jed > > Observation: Using IS does not give all desired generality since the new field may be a linear combination of DOFs. Or, yikes, a nonlinear transformation of DOF. > Thus we may need something like DMTransformVec() to do these transformations. > > > So this moves us beyond the traditional block linear algebra methods but if we could solve it we'd have some very nice new powerful capabilities on our hands. > > Question 1: Do we need to solve this for the current generation of PCFieldsplit? > Question 2: Do we want to solve it in PCFieldSplit or should that come later in a different PC, PCBasisChangeThingy or something? > > Yes, but what DM should PCBasisChangeThingy forward to the inner PC? Presumably a "new" DM that represents the new set of equations ... and to start we will need to do very specific examples before coming up with a general paradigm. But it seems best to me to make this orthogonal to PCFieldSplit. Barry

1 0