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

Re: [petsc-dev] build on cygwin for cygwin
by Satish Balay 28 Feb '17

28 Feb '17

On Tue, 28 Feb 2017, Marco Atzeri wrote: > On 27/02/2017 19:56, Satish Balay wrote: > > On Mon, 27 Feb 2017, Marco Atzeri wrote: > > > > > On 27/02/2017 19:15, Satish Balay wrote: > > > > Ah - previously I didn't check the error log closely. > > > > > > > > > Hm - the test for sequential build failed? Perhaps you had some mixup > > > > of files with the previous openmpi build. > > > > When changing builds (like switching mpi) - its best to remove build > > > > files in PETSC_ARCH - and install files in prefix-dir - before the > > > > next build. Or use a different PETSC_ARCH/prefix values for the > > > > second build. > > > > > > I am using in general different directories. > > > I will try again just in case I missed something. > > > > > > > I tried both builds - and they ran fine for me. [with maint/3.7] > > static build works fine, including the stream test. I was refering to the dll build [and with-mpi=0]. The C examples run fine for me. > > > > > > > > > Either way - fortran examples don't work with dlls - so its best to > > > > use --with-shared-libraries=0 > > > > > > don't work in general or just for cygwin ? > > > > PETSc fortran interface uses fortran common blocks - and DLLs (with > > common blocks) require some compiler directives (perhaps compiler > > specific wrt ifort vs gfortran) to handle this. > > > > We don't yet have that in petsc. > > I assume it is related to you inquiry on the cygwin mailing list > https://cygwin.com/ml/cygwin/2015-06/msg00352.html > yes > Had you any progress or the matter or is it still there the blocking point ? > Just to know what to investigate as simple test case. This is where it stands. If you are able to make progress - you can submit a PR thanks, Satish > > > Satish > > Regards > Marco > > >

1 0

Re: [petsc-dev] [petsc-users] Multigrid with defect correction
by Barry Smith 28 Feb '17

28 Feb '17

> On Feb 28, 2017, at 11:16 AM, Dave May <dave.mayhem23(a)gmail.com> wrote: > > > > PCMG can actually do this for you if you tell it to use Galerkin so you shouldn't need to do this (with both the "standard" PCMG interface and its use of DM). > > > Sure. I realize if there is a DM on the fine grid, and the DM has the method DMCoarsen(), then galerkin operators can/will be constructed by PCMG if the galerkin flag is set. However, my understanding is that PCMG will create Galerkin operators for on ALL levels (except the finest). > > I wanted to have control over exactly which levels in the hierarchy defined Galerkin operators. I also wanted more control over how the mesh hierarchy was constructed than DMCoarsen() provided, so i found it more convenient to build the mesh hierarchy myself. I understand and that is fine. We have these two "interfaces" for PCMG, the "loop over the levels and do what you want on each PCMG interface" that has the most flexibility and the "provide a callback that can set some things for each level (provide the DM) interface". This second one has more limitations, and likely will always have more limitations but I would like to make sure it does not have unneeded limitations. > > Cheers > Dave > > > Barry > > > > > > > > > >> I would ask, why just two discretizations? I've always thought a better > > >> interface would be for Mats (and any other operators/functions) to have > > >> optional approximations or supplementary data attached. We can do this > > >> with PetscObjectCompose, but that's hard to work with in a structured > > >> way. Anyway, I think I would rather just have the Amat with ability to > > >> attach one or more Pmats. > > > > > > Understood. I don't like the API as "attach one or more Pmats to a > > > Mat" but instead prefer to think of it as a linear operator that has > > > one or more "levels of approximation" > > > > Should they be "levels" (numbered?) or named approximations? > > > > > that can optionally be indicated; which one you wish to use is up to > > > the algorithm using the linear operator. Changing this in PETSc > > > would require some conceptional changes and refactoring in > > > PETSc. Would simplify things like SNESSetJacobian() because we would > > > not need to track mat and pmat. But would not be terribly difficult > > > for the PETSc developers or user. > > > > > > If we used this model, instead of adding a new DM to the interface > > > do we keep the one DM but have a way to "attach" other "lower order" > > > DMs to the that one DM that can be accessed? > > > > Presumably. Is a DM a function space, an operator, or both? > > > > > DMCreateMatrix(dm,&mat); > > > pmat = mat; > > > DMGetLowerDM(dm,&ldm); > > > if (ldm) { > > > DMCreateMatrix(ldm,&pmat); > > > } > > > KSP/SNES/TSSetIJacobian(ts,mat,pmat,...); > > > > > > more generally DMGetLowerDM(dm,PetscOrderIndicater i,&ldm); Similarly MatGetLowerMat(mat,PetscOrderIndicater i,&lmat)? > > > > > > ----- an aside > > > BTW: this introduces a possible "missing" object in PETSc. We have > > > as a fundamental object the Mat which is a linear operator, we do > > > not have a fundamental object for a "general" operator, instead we > > > accept user provided functions directly (SNESSetFunction(), > > > TSSetIFunction() etc). I've resisted the PetscOperator because I > > > felt it was overly abstract (and thus confusing to most users) and > > > is more difficult to implement cleanly since the operator may depend > > > on more than one input vector (TSSetIFunction depends on u and udot; > > > would we have a PetscOperator2 that had two input vectors?). > > > > DMSNES and DMTS effectively provide this, it's just intimately coupled > > to the DM. I advocated before for separating these because function > > spaces have meaning independent of operators defined on them, but that's > > more of a performance consideration that assumes DM implementations > > don't have cheap sharing of common data. > > > > > I did try the PF object but it is very limited and rarely > > > used. Certainly I think for a wide class of users, especially Fortran, > > > SNESSetFunction() is understandable but > > > > > > PetscOperatorCreate(&operator) > > > PetscOperatorSetRawFunction(operator, function()) > > > SNESSetOperator(snes,operator) > > > > > > becomes gibberish and so we would lose many potential users. Of course we could introduce the PetscOperator and maintain the simplicity of the current interface by having things like SNESSetFunction() automatically handle the conversion to the "operator" interface and have a new SNESSetOperator() for the general case. I hate having this kind of "wrapper" code because it usually bit decays over time and I don't like have multiple equivalent interfaces for the same effect (how do you document them, how does the user decide which one to use? I hate perl :-), there should be one way to do each thing). > > > > SNESSetFunction calls DMSNESSetFunction. > > > > > Or does one every need "lower order operation evaluation"? Is lower > > > order linearization enough? > > > > No, dual-order discretizations are possibly more important for nonlinear > > problems. > > > > > Currently with mat and pmat we only support lower order linearization, > > > we don't for example allow passing a higher order and lower order > > > nonlinear function both to SNESSetFunction. If the lower order > > > operator is never needed, only the lower order linearization then > > > PetscOperator (which would allow managing several orders of general > > > operator) would not be needed so API for users is simpler. Or we could > > > introduce lower order linear operators now and not worry about lower > > > order function evaluations for years? --- end of the aside > > > > One common application of lower order discretizations is in in the > > function provided to SNESSetNGS. > > > > > I also don't disagree with Matt, should we just "complete" support for two levels of approximation by having proper two DMs for TS/KSP/SNES, see how that goes rather than jumping directly from a (far from perfect since users need to manually do ugly stuff since we only have one DM) two level support (as we currently have) to multilevel support? I won't change the current interface to do this instead I would just add TS/SNES/KSPSetDMs() > > > > Yup. > > > > > Imagine poor Ed's book if we suddenly change stuff like this in PETSc :-( So the question is what should/could we do now to improve flexibility without requiring Ed to start all over again? That is, how gradually we need to introduce this kind of overhaul, and with how much backward compatibility? > > > > I would say quite gradual and with at least a full release cycle of > > backward compatibility. And although I think it is ultimately better > > for extensibility, I'm not 100% convinced it is worth the cost to users > > to overhaul. One nice thing about the current interface is that it > > forces people to acknowledge that those two arguments could be > > different, so they may read or find examples. If we switch so our > > operators have a basket of approximations, lots of users may never learn > > about that basket. > > > >

1 0

Re: [petsc-dev] [petsc-users] Multigrid with defect correction
by Dave May 28 Feb '17

28 Feb '17

> > PCMG can actually do this for you if you tell it to use Galerkin so you > shouldn't need to do this (with both the "standard" PCMG interface and its > use of DM). > > Sure. I realize if there is a DM on the fine grid, and the DM has the method DMCoarsen(), then galerkin operators can/will be constructed by PCMG if the galerkin flag is set. However, my understanding is that PCMG will create Galerkin operators for on ALL levels (except the finest). I wanted to have control over exactly which levels in the hierarchy defined Galerkin operators. I also wanted more control over how the mesh hierarchy was constructed than DMCoarsen() provided, so i found it more convenient to build the mesh hierarchy myself. Cheers Dave > Barry > > > > > > > > > >> I would ask, why just two discretizations? I've always thought a > better > > >> interface would be for Mats (and any other operators/functions) to > have > > >> optional approximations or supplementary data attached. We can do > this > > >> with PetscObjectCompose, but that's hard to work with in a structured > > >> way. Anyway, I think I would rather just have the Amat with ability > to > > >> attach one or more Pmats. > > > > > > Understood. I don't like the API as "attach one or more Pmats to a > > > Mat" but instead prefer to think of it as a linear operator that has > > > one or more "levels of approximation" > > > > Should they be "levels" (numbered?) or named approximations? > > > > > that can optionally be indicated; which one you wish to use is up to > > > the algorithm using the linear operator. Changing this in PETSc > > > would require some conceptional changes and refactoring in > > > PETSc. Would simplify things like SNESSetJacobian() because we would > > > not need to track mat and pmat. But would not be terribly difficult > > > for the PETSc developers or user. > > > > > > If we used this model, instead of adding a new DM to the interface > > > do we keep the one DM but have a way to "attach" other "lower order" > > > DMs to the that one DM that can be accessed? > > > > Presumably. Is a DM a function space, an operator, or both? > > > > > DMCreateMatrix(dm,&mat); > > > pmat = mat; > > > DMGetLowerDM(dm,&ldm); > > > if (ldm) { > > > DMCreateMatrix(ldm,&pmat); > > > } > > > KSP/SNES/TSSetIJacobian(ts,mat,pmat,...); > > > > > > more generally DMGetLowerDM(dm,PetscOrderIndicater i,&ldm); > Similarly MatGetLowerMat(mat,PetscOrderIndicater i,&lmat)? > > > > > > ----- an aside > > > BTW: this introduces a possible "missing" object in PETSc. We have > > > as a fundamental object the Mat which is a linear operator, we do > > > not have a fundamental object for a "general" operator, instead we > > > accept user provided functions directly (SNESSetFunction(), > > > TSSetIFunction() etc). I've resisted the PetscOperator because I > > > felt it was overly abstract (and thus confusing to most users) and > > > is more difficult to implement cleanly since the operator may depend > > > on more than one input vector (TSSetIFunction depends on u and udot; > > > would we have a PetscOperator2 that had two input vectors?). > > > > DMSNES and DMTS effectively provide this, it's just intimately coupled > > to the DM. I advocated before for separating these because function > > spaces have meaning independent of operators defined on them, but that's > > more of a performance consideration that assumes DM implementations > > don't have cheap sharing of common data. > > > > > I did try the PF object but it is very limited and rarely > > > used. Certainly I think for a wide class of users, especially Fortran, > > > SNESSetFunction() is understandable but > > > > > > PetscOperatorCreate(&operator) > > > PetscOperatorSetRawFunction(operator, function()) > > > SNESSetOperator(snes,operator) > > > > > > becomes gibberish and so we would lose many potential users. Of course > we could introduce the PetscOperator and maintain the simplicity of the > current interface by having things like SNESSetFunction() automatically > handle the conversion to the "operator" interface and have a new > SNESSetOperator() for the general case. I hate having this kind of > "wrapper" code because it usually bit decays over time and I don't like > have multiple equivalent interfaces for the same effect (how do you > document them, how does the user decide which one to use? I hate perl :-), > there should be one way to do each thing). > > > > SNESSetFunction calls DMSNESSetFunction. > > > > > Or does one every need "lower order operation evaluation"? Is lower > > > order linearization enough? > > > > No, dual-order discretizations are possibly more important for nonlinear > > problems. > > > > > Currently with mat and pmat we only support lower order linearization, > > > we don't for example allow passing a higher order and lower order > > > nonlinear function both to SNESSetFunction. If the lower order > > > operator is never needed, only the lower order linearization then > > > PetscOperator (which would allow managing several orders of general > > > operator) would not be needed so API for users is simpler. Or we could > > > introduce lower order linear operators now and not worry about lower > > > order function evaluations for years? --- end of the aside > > > > One common application of lower order discretizations is in in the > > function provided to SNESSetNGS. > > > > > I also don't disagree with Matt, should we just "complete" support for > two levels of approximation by having proper two DMs for TS/KSP/SNES, see > how that goes rather than jumping directly from a (far from perfect since > users need to manually do ugly stuff since we only have one DM) two level > support (as we currently have) to multilevel support? I won't change the > current interface to do this instead I would just add TS/SNES/KSPSetDMs() > > > > Yup. > > > > > Imagine poor Ed's book if we suddenly change stuff like this in PETSc > :-( So the question is what should/could we do now to improve flexibility > without requiring Ed to start all over again? That is, how gradually we > need to introduce this kind of overhaul, and with how much backward > compatibility? > > > > I would say quite gradual and with at least a full release cycle of > > backward compatibility. And although I think it is ultimately better > > for extensibility, I'm not 100% convinced it is worth the cost to users > > to overhaul. One nice thing about the current interface is that it > > forces people to acknowledge that those two arguments could be > > different, so they may read or find examples. If we switch so our > > operators have a basket of approximations, lots of users may never learn > > about that basket. > > > >

1 0

Re: [petsc-dev] [petsc-users] Multigrid with defect correction
by Barry Smith 28 Feb '17

28 Feb '17

> On Feb 28, 2017, at 4:03 AM, Dave May <dave.mayhem23(a)gmail.com> wrote: > > > > On 26 February 2017 at 05:30, Jed Brown <jedbrown(a)mcs.anl.gov> wrote: > Barry Smith <bsmith(a)mcs.anl.gov> writes: > > >> On Feb 25, 2017, at 7:35 PM, Jed Brown <jedbrown(a)mcs.anl.gov> wrote: > >> > >> Barry Smith <bsmith(a)mcs.anl.gov> writes: > >>> Do you think this is a reasonable approach or am I missing > >>> something fundamental? I am assuming generally for the "higher > >>> order" DM the Mat it returns is a MATSHELL or a new matrix class > >>> built on "tensor contractions" and that kind of nonsense. I don't > >>> want to do all the coding and then have it turn out that it is > >>> totally useless for CEED etc. > >> > >> Well, this is exactly what we do in pTatin. > > > > So why didn't you fix it in PETSc at that time? > > I didn't know it "didn't work" because what we were doing did work. (I > think we used PCMG interfaces directly; we were mixing unassembled and > assembled representations, plus Galerkin on coarser levels.) > > > I can confirm that Jed's recollection is correct. > > The configuration used in ptatin was performed in the context of a non-linear problem. Essentially all the MG setup has conducted during my method for computeJacobian(). The setup / config for the different operators was a little delicate, but the standard PCMG API permitted mixing matrix-free, assembled and Galerkin coarse operators. Sure the PCMG API has the flexibility; it is the DM usage in PCMG etc that only has a single DM and no way to handle mat and pmat separately. > For the Galerkin case, I performed the call to MatPtAP myself in my computeJacobian() function and would stuff the result into the appropriate KSP defining the smoother. PCMG can actually do this for you if you tell it to use Galerkin so you shouldn't need to do this (with both the "standard" PCMG interface and its use of DM). Barry > > > > >> I would ask, why just two discretizations? I've always thought a better > >> interface would be for Mats (and any other operators/functions) to have > >> optional approximations or supplementary data attached. We can do this > >> with PetscObjectCompose, but that's hard to work with in a structured > >> way. Anyway, I think I would rather just have the Amat with ability to > >> attach one or more Pmats. > > > > Understood. I don't like the API as "attach one or more Pmats to a > > Mat" but instead prefer to think of it as a linear operator that has > > one or more "levels of approximation" > > Should they be "levels" (numbered?) or named approximations? > > > that can optionally be indicated; which one you wish to use is up to > > the algorithm using the linear operator. Changing this in PETSc > > would require some conceptional changes and refactoring in > > PETSc. Would simplify things like SNESSetJacobian() because we would > > not need to track mat and pmat. But would not be terribly difficult > > for the PETSc developers or user. > > > > If we used this model, instead of adding a new DM to the interface > > do we keep the one DM but have a way to "attach" other "lower order" > > DMs to the that one DM that can be accessed? > > Presumably. Is a DM a function space, an operator, or both? > > > DMCreateMatrix(dm,&mat); > > pmat = mat; > > DMGetLowerDM(dm,&ldm); > > if (ldm) { > > DMCreateMatrix(ldm,&pmat); > > } > > KSP/SNES/TSSetIJacobian(ts,mat,pmat,...); > > > > more generally DMGetLowerDM(dm,PetscOrderIndicater i,&ldm); Similarly MatGetLowerMat(mat,PetscOrderIndicater i,&lmat)? > > > > ----- an aside > > BTW: this introduces a possible "missing" object in PETSc. We have > > as a fundamental object the Mat which is a linear operator, we do > > not have a fundamental object for a "general" operator, instead we > > accept user provided functions directly (SNESSetFunction(), > > TSSetIFunction() etc). I've resisted the PetscOperator because I > > felt it was overly abstract (and thus confusing to most users) and > > is more difficult to implement cleanly since the operator may depend > > on more than one input vector (TSSetIFunction depends on u and udot; > > would we have a PetscOperator2 that had two input vectors?). > > DMSNES and DMTS effectively provide this, it's just intimately coupled > to the DM. I advocated before for separating these because function > spaces have meaning independent of operators defined on them, but that's > more of a performance consideration that assumes DM implementations > don't have cheap sharing of common data. > > > I did try the PF object but it is very limited and rarely > > used. Certainly I think for a wide class of users, especially Fortran, > > SNESSetFunction() is understandable but > > > > PetscOperatorCreate(&operator) > > PetscOperatorSetRawFunction(operator, function()) > > SNESSetOperator(snes,operator) > > > > becomes gibberish and so we would lose many potential users. Of course we could introduce the PetscOperator and maintain the simplicity of the current interface by having things like SNESSetFunction() automatically handle the conversion to the "operator" interface and have a new SNESSetOperator() for the general case. I hate having this kind of "wrapper" code because it usually bit decays over time and I don't like have multiple equivalent interfaces for the same effect (how do you document them, how does the user decide which one to use? I hate perl :-), there should be one way to do each thing). > > SNESSetFunction calls DMSNESSetFunction. > > > Or does one every need "lower order operation evaluation"? Is lower > > order linearization enough? > > No, dual-order discretizations are possibly more important for nonlinear > problems. > > > Currently with mat and pmat we only support lower order linearization, > > we don't for example allow passing a higher order and lower order > > nonlinear function both to SNESSetFunction. If the lower order > > operator is never needed, only the lower order linearization then > > PetscOperator (which would allow managing several orders of general > > operator) would not be needed so API for users is simpler. Or we could > > introduce lower order linear operators now and not worry about lower > > order function evaluations for years? --- end of the aside > > One common application of lower order discretizations is in in the > function provided to SNESSetNGS. > > > I also don't disagree with Matt, should we just "complete" support for two levels of approximation by having proper two DMs for TS/KSP/SNES, see how that goes rather than jumping directly from a (far from perfect since users need to manually do ugly stuff since we only have one DM) two level support (as we currently have) to multilevel support? I won't change the current interface to do this instead I would just add TS/SNES/KSPSetDMs() > > Yup. > > > Imagine poor Ed's book if we suddenly change stuff like this in PETSc :-( So the question is what should/could we do now to improve flexibility without requiring Ed to start all over again? That is, how gradually we need to introduce this kind of overhaul, and with how much backward compatibility? > > I would say quite gradual and with at least a full release cycle of > backward compatibility. And although I think it is ultimately better > for extensibility, I'm not 100% convinced it is worth the cost to users > to overhaul. One nice thing about the current interface is that it > forces people to acknowledge that those two arguments could be > different, so they may read or find examples. If we switch so our > operators have a basket of approximations, lots of users may never learn > about that basket. >

1 0

Re: [petsc-dev] no petsc on Edison
by Hong 28 Feb '17

28 Feb '17

This work is merged to master https://bitbucket.org/petsc/petsc/commits/57c33550d09512dabc9444743ef186e1d… Hong On Tue, Feb 28, 2017 at 8:04 AM, Mark Adams <mfadams(a)lbl.gov> wrote: > Yes, this looks great. Good to get some data on this. Hard problem. > > On Mon, Feb 27, 2017 at 3:31 PM, Barry Smith <bsmith(a)mcs.anl.gov> wrote: > >> >> Hong, >> >> Very nice, thanks! >> >> Barry >> >> >> > On Feb 27, 2017, at 2:06 PM, Hong <hzhang(a)mcs.anl.gov> wrote: >> > >> > Mark, >> > I implemented scalable MatPtAP and did comparisons of three >> implementations using ex56.c on alcf cetus machine (this machine has small >> memory, 1GB/core): >> > - nonscalable PtAP: use an array of length PN to do dense axpy >> > - scalable PtAP: do sparse axpy without use of PN array >> > - hypre PtAP. >> > >> > The results are attached. Summary: >> > - nonscalable PtAP is 2x faster than scalable, 8x faster than hypre PtAP >> > - scalable PtAP is 4x faster than hypre PtAP >> > - hypre uses less memory (see job.ne399.n63.np1000.sh) >> > >> > Based on above observation, I set the default PtAP algorithm as >> 'nonscalable'. >> > When PN > local estimated nonzero of C=PtAP, then switch default to >> 'scalable'. >> > User can overwrite default. >> > >> > For the case of np=8000, ne=599 (see job.ne599.n500.np8000.sh), I get >> > MatPtAP 3.6224e+01 (nonscalable for small mats, >> scalable for larger ones) >> > scalable MatPtAP 4.6129e+01 >> > hypre 1.9389e+02 >> > >> > I'm merging this work to next, then to master soon. >> > >> > Hong >> > >> > On Wed, Feb 8, 2017 at 7:25 PM, Mark Adams <mfadams(a)lbl.gov> wrote: >> > Thanks Hong, that sounds great. >> > >> > I am weary of silent optimizations like you suggest but 2x is big! and >> failure is very bad. So I would vote for your suggestion. >> > >> > ex56 is an elasticity problem. It would be nice, now that you have this >> experimental setup in place, to compare with hypre on a 3D scalar problem. >> Hypre might not spend much effort optimizing for block matrices. 3x better >> than hypre seems large to me. I have to suspect that hypre does not >> exploit blocked matrices as much as we do. >> > >> > Thanks again, >> > >> > >> > >> > On Wed, Feb 8, 2017 at 12:07 PM, Hong <hzhang(a)mcs.anl.gov> wrote: >> > I conducted tests on MatMatMult() and MatPtAP() using >> petsc/src/ksp/ksp/examples/tutorials/ex56.c (gamg) on a 8-core machine >> (petsc machine). The output file is attached. >> > >> > Summary: >> > 1) non-scalable MatMatMult() for mpiaij format is 2x faster than >> scalable version. The major difference between the two is dense-axpy vs. >> sparse-axpy. >> > >> > Currently, we set non-scalable as default, which leads to problem when >> running large problems. >> > How about setting default as >> > - non-scalable for small to medium size matrices >> > - scalable for larger ones, e.g. >> > >> > + ierr = PetscOptionsEList("-matmatmult_via","Algorithmic >> approach","MatMatMult",algTypes,nalg,algTypes[1],&alg,&flg); >> > >> > + if (!flg) { /* set default algorithm based on B->cmap->N */ >> > + PetscMPIInt size; >> > + ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); >> > + if ((PetscReal)(B->cmap->N)/size > 100000.0) alg = 0; /* >> scalable algorithm */ >> > + } >> > >> > i.e., if user does NOT pick an algorithm, when ave cols per process > >> 100k, use scalable implementation; otherwise, non-scalable version. >> > >> > 2) We do NOT have scalable implementation for MatPtAP() yet. >> > We have non-scalable PtAP and interface to Hypre's PtAP. Comparing the >> two, >> > Petsc MatPtAP() is approx 3x faster than Hypre's. >> > >> > I'm writing a scalable MatPtAP() now. >> > >> > Hong >> > >> > On Thu, Feb 2, 2017 at 2:54 PM, Stefano Zampini < >> stefano.zampini(a)gmail.com> wrote: >> > >> > >> > Il 02 Feb 2017 23:43, "Mark Adams" <mfadams(a)lbl.gov> ha scritto: >> > >> > >> > On Thu, Feb 2, 2017 at 12:02 PM, Stefano Zampini < >> stefano.zampini(a)gmail.com> wrote: >> > Mark, >> > >> > I saw your configuration has hypre. If you could run with master, you >> may try -matptap_via hypre. >> > >> > This is worth trying. Does this even work with GAMG? >> > >> > Yes, it should work, except that the block sizes, if any, are not >> propagated to the resulting matrix. I can add it if you need it. >> > >> > >> > >> > Treb: try hypre anyway. It has its own RAP code. >> > >> > >> > With that option, you will use hypre's RAP with MATAIJ >> > >> > >> > It uses BoomerAMGBuildCoarseOperator directly with the AIJ matrices. >> > >> > Stefano >> > >> >> On Feb 2, 2017, at 7:28 PM, Mark Adams <mfadams(a)lbl.gov> wrote: >> >> >> >> >> >> >> >> On Thu, Feb 2, 2017 at 11:13 AM, Hong <hzhang(a)mcs.anl.gov> wrote: >> >> Mark: >> >> Try '-matmatmult_via scalable' first. If this works, should we set it >> as default? >> >> >> >> If it is robust I would say yes unless it is noticeably slower (say >> >20%) small scale problems. >> >> >> > >> > >> > >> > >> > >> > >> > <out_ex56_cetus_short> >> >> >

1 0

Re: [petsc-dev] no petsc on Edison
by Mark Adams 28 Feb '17

28 Feb '17

Yes, this looks great. Good to get some data on this. Hard problem. On Mon, Feb 27, 2017 at 3:31 PM, Barry Smith <bsmith(a)mcs.anl.gov> wrote: > > Hong, > > Very nice, thanks! > > Barry > > > On Feb 27, 2017, at 2:06 PM, Hong <hzhang(a)mcs.anl.gov> wrote: > > > > Mark, > > I implemented scalable MatPtAP and did comparisons of three > implementations using ex56.c on alcf cetus machine (this machine has small > memory, 1GB/core): > > - nonscalable PtAP: use an array of length PN to do dense axpy > > - scalable PtAP: do sparse axpy without use of PN array > > - hypre PtAP. > > > > The results are attached. Summary: > > - nonscalable PtAP is 2x faster than scalable, 8x faster than hypre PtAP > > - scalable PtAP is 4x faster than hypre PtAP > > - hypre uses less memory (see job.ne399.n63.np1000.sh) > > > > Based on above observation, I set the default PtAP algorithm as > 'nonscalable'. > > When PN > local estimated nonzero of C=PtAP, then switch default to > 'scalable'. > > User can overwrite default. > > > > For the case of np=8000, ne=599 (see job.ne599.n500.np8000.sh), I get > > MatPtAP 3.6224e+01 (nonscalable for small mats, > scalable for larger ones) > > scalable MatPtAP 4.6129e+01 > > hypre 1.9389e+02 > > > > I'm merging this work to next, then to master soon. > > > > Hong > > > > On Wed, Feb 8, 2017 at 7:25 PM, Mark Adams <mfadams(a)lbl.gov> wrote: > > Thanks Hong, that sounds great. > > > > I am weary of silent optimizations like you suggest but 2x is big! and > failure is very bad. So I would vote for your suggestion. > > > > ex56 is an elasticity problem. It would be nice, now that you have this > experimental setup in place, to compare with hypre on a 3D scalar problem. > Hypre might not spend much effort optimizing for block matrices. 3x better > than hypre seems large to me. I have to suspect that hypre does not > exploit blocked matrices as much as we do. > > > > Thanks again, > > > > > > > > On Wed, Feb 8, 2017 at 12:07 PM, Hong <hzhang(a)mcs.anl.gov> wrote: > > I conducted tests on MatMatMult() and MatPtAP() using > petsc/src/ksp/ksp/examples/tutorials/ex56.c (gamg) on a 8-core machine > (petsc machine). The output file is attached. > > > > Summary: > > 1) non-scalable MatMatMult() for mpiaij format is 2x faster than > scalable version. The major difference between the two is dense-axpy vs. > sparse-axpy. > > > > Currently, we set non-scalable as default, which leads to problem when > running large problems. > > How about setting default as > > - non-scalable for small to medium size matrices > > - scalable for larger ones, e.g. > > > > + ierr = PetscOptionsEList("-matmatmult_via","Algorithmic > approach","MatMatMult",algTypes,nalg,algTypes[1],&alg,&flg); > > > > + if (!flg) { /* set default algorithm based on B->cmap->N */ > > + PetscMPIInt size; > > + ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); > > + if ((PetscReal)(B->cmap->N)/size > 100000.0) alg = 0; /* scalable > algorithm */ > > + } > > > > i.e., if user does NOT pick an algorithm, when ave cols per process > > 100k, use scalable implementation; otherwise, non-scalable version. > > > > 2) We do NOT have scalable implementation for MatPtAP() yet. > > We have non-scalable PtAP and interface to Hypre's PtAP. Comparing the > two, > > Petsc MatPtAP() is approx 3x faster than Hypre's. > > > > I'm writing a scalable MatPtAP() now. > > > > Hong > > > > On Thu, Feb 2, 2017 at 2:54 PM, Stefano Zampini < > stefano.zampini(a)gmail.com> wrote: > > > > > > Il 02 Feb 2017 23:43, "Mark Adams" <mfadams(a)lbl.gov> ha scritto: > > > > > > On Thu, Feb 2, 2017 at 12:02 PM, Stefano Zampini < > stefano.zampini(a)gmail.com> wrote: > > Mark, > > > > I saw your configuration has hypre. If you could run with master, you > may try -matptap_via hypre. > > > > This is worth trying. Does this even work with GAMG? > > > > Yes, it should work, except that the block sizes, if any, are not > propagated to the resulting matrix. I can add it if you need it. > > > > > > > > Treb: try hypre anyway. It has its own RAP code. > > > > > > With that option, you will use hypre's RAP with MATAIJ > > > > > > It uses BoomerAMGBuildCoarseOperator directly with the AIJ matrices. > > > > Stefano > > > >> On Feb 2, 2017, at 7:28 PM, Mark Adams <mfadams(a)lbl.gov> wrote: > >> > >> > >> > >> On Thu, Feb 2, 2017 at 11:13 AM, Hong <hzhang(a)mcs.anl.gov> wrote: > >> Mark: > >> Try '-matmatmult_via scalable' first. If this works, should we set it > as default? > >> > >> If it is robust I would say yes unless it is noticeably slower (say > >20%) small scale problems. > >> > > > > > > > > > > > > > > <out_ex56_cetus_short> > >

1 0

Re: [petsc-dev] [petsc-users] Multigrid with defect correction
by Dave May 28 Feb '17

28 Feb '17

On 26 February 2017 at 05:30, Jed Brown <jedbrown(a)mcs.anl.gov> wrote: > Barry Smith <bsmith(a)mcs.anl.gov> writes: > > >> On Feb 25, 2017, at 7:35 PM, Jed Brown <jedbrown(a)mcs.anl.gov> wrote: > >> > >> Barry Smith <bsmith(a)mcs.anl.gov> writes: > >>> Do you think this is a reasonable approach or am I missing > >>> something fundamental? I am assuming generally for the "higher > >>> order" DM the Mat it returns is a MATSHELL or a new matrix class > >>> built on "tensor contractions" and that kind of nonsense. I don't > >>> want to do all the coding and then have it turn out that it is > >>> totally useless for CEED etc. > >> > >> Well, this is exactly what we do in pTatin. > > > > So why didn't you fix it in PETSc at that time? > > I didn't know it "didn't work" because what we were doing did work. (I > think we used PCMG interfaces directly; we were mixing unassembled and > assembled representations, plus Galerkin on coarser levels.) > I can confirm that Jed's recollection is correct. The configuration used in ptatin was performed in the context of a non-linear problem. Essentially all the MG setup has conducted during my method for computeJacobian(). The setup / config for the different operators was a little delicate, but the standard PCMG API permitted mixing matrix-free, assembled and Galerkin coarse operators. For the Galerkin case, I performed the call to MatPtAP myself in my computeJacobian() function and would stuff the result into the appropriate KSP defining the smoother. > > >> I would ask, why just two discretizations? I've always thought a better > >> interface would be for Mats (and any other operators/functions) to have > >> optional approximations or supplementary data attached. We can do this > >> with PetscObjectCompose, but that's hard to work with in a structured > >> way. Anyway, I think I would rather just have the Amat with ability to > >> attach one or more Pmats. > > > > Understood. I don't like the API as "attach one or more Pmats to a > > Mat" but instead prefer to think of it as a linear operator that has > > one or more "levels of approximation" > > Should they be "levels" (numbered?) or named approximations? > > > that can optionally be indicated; which one you wish to use is up to > > the algorithm using the linear operator. Changing this in PETSc > > would require some conceptional changes and refactoring in > > PETSc. Would simplify things like SNESSetJacobian() because we would > > not need to track mat and pmat. But would not be terribly difficult > > for the PETSc developers or user. > > > > If we used this model, instead of adding a new DM to the interface > > do we keep the one DM but have a way to "attach" other "lower order" > > DMs to the that one DM that can be accessed? > > Presumably. Is a DM a function space, an operator, or both? > > > DMCreateMatrix(dm,&mat); > > pmat = mat; > > DMGetLowerDM(dm,&ldm); > > if (ldm) { > > DMCreateMatrix(ldm,&pmat); > > } > > KSP/SNES/TSSetIJacobian(ts,mat,pmat,...); > > > > more generally DMGetLowerDM(dm,PetscOrderIndicater i,&ldm); Similarly > MatGetLowerMat(mat,PetscOrderIndicater i,&lmat)? > > > > ----- an aside > > BTW: this introduces a possible "missing" object in PETSc. We have > > as a fundamental object the Mat which is a linear operator, we do > > not have a fundamental object for a "general" operator, instead we > > accept user provided functions directly (SNESSetFunction(), > > TSSetIFunction() etc). I've resisted the PetscOperator because I > > felt it was overly abstract (and thus confusing to most users) and > > is more difficult to implement cleanly since the operator may depend > > on more than one input vector (TSSetIFunction depends on u and udot; > > would we have a PetscOperator2 that had two input vectors?). > > DMSNES and DMTS effectively provide this, it's just intimately coupled > to the DM. I advocated before for separating these because function > spaces have meaning independent of operators defined on them, but that's > more of a performance consideration that assumes DM implementations > don't have cheap sharing of common data. > > > I did try the PF object but it is very limited and rarely > > used. Certainly I think for a wide class of users, especially Fortran, > > SNESSetFunction() is understandable but > > > > PetscOperatorCreate(&operator) > > PetscOperatorSetRawFunction(operator, function()) > > SNESSetOperator(snes,operator) > > > > becomes gibberish and so we would lose many potential users. Of course > we could introduce the PetscOperator and maintain the simplicity of the > current interface by having things like SNESSetFunction() automatically > handle the conversion to the "operator" interface and have a new > SNESSetOperator() for the general case. I hate having this kind of > "wrapper" code because it usually bit decays over time and I don't like > have multiple equivalent interfaces for the same effect (how do you > document them, how does the user decide which one to use? I hate perl :-), > there should be one way to do each thing). > > SNESSetFunction calls DMSNESSetFunction. > > > Or does one every need "lower order operation evaluation"? Is lower > > order linearization enough? > > No, dual-order discretizations are possibly more important for nonlinear > problems. > > > Currently with mat and pmat we only support lower order linearization, > > we don't for example allow passing a higher order and lower order > > nonlinear function both to SNESSetFunction. If the lower order > > operator is never needed, only the lower order linearization then > > PetscOperator (which would allow managing several orders of general > > operator) would not be needed so API for users is simpler. Or we could > > introduce lower order linear operators now and not worry about lower > > order function evaluations for years? --- end of the aside > > One common application of lower order discretizations is in in the > function provided to SNESSetNGS. > > > I also don't disagree with Matt, should we just "complete" support for > two levels of approximation by having proper two DMs for TS/KSP/SNES, see > how that goes rather than jumping directly from a (far from perfect since > users need to manually do ugly stuff since we only have one DM) two level > support (as we currently have) to multilevel support? I won't change the > current interface to do this instead I would just add TS/SNES/KSPSetDMs() > > Yup. > > > Imagine poor Ed's book if we suddenly change stuff like this in PETSc > :-( So the question is what should/could we do now to improve flexibility > without requiring Ed to start all over again? That is, how gradually we > need to introduce this kind of overhaul, and with how much backward > compatibility? > > I would say quite gradual and with at least a full release cycle of > backward compatibility. And although I think it is ultimately better > for extensibility, I'm not 100% convinced it is worth the cost to users > to overhaul. One nice thing about the current interface is that it > forces people to acknowledge that those two arguments could be > different, so they may read or find examples. If we switch so our > operators have a basket of approximations, lots of users may never learn > about that basket. >

1 0

Re: [petsc-dev] build on cygwin for cygwin
by Marco Atzeri 28 Feb '17

28 Feb '17

On 27/02/2017 19:56, Satish Balay wrote: > On Mon, 27 Feb 2017, Marco Atzeri wrote: > >> On 27/02/2017 19:15, Satish Balay wrote: >>> Ah - previously I didn't check the error log closely. >>> >>> Hm - the test for sequential build failed? Perhaps you had some mixup >>> of files with the previous openmpi build. >>> When changing builds (like switching mpi) - its best to remove build >>> files in PETSC_ARCH - and install files in prefix-dir - before the >>> next build. Or use a different PETSC_ARCH/prefix values for the >>> second build. >> >> I am using in general different directories. >> I will try again just in case I missed something. >> >>> I tried both builds - and they ran fine for me. [with maint/3.7] static build works fine, including the stream test. >>> >>> Either way - fortran examples don't work with dlls - so its best to >>> use --with-shared-libraries=0 >> >> don't work in general or just for cygwin ? > > PETSc fortran interface uses fortran common blocks - and DLLs (with > common blocks) require some compiler directives (perhaps compiler > specific wrt ifort vs gfortran) to handle this. > > We don't yet have that in petsc. I assume it is related to you inquiry on the cygwin mailing list https://cygwin.com/ml/cygwin/2015-06/msg00352.html Had you any progress or the matter or is it still there the blocking point ? Just to know what to investigate as simple test case. > Satish Regards Marco

1 0

Re: [petsc-dev] no petsc on Edison
by Barry Smith 27 Feb '17

27 Feb '17

Hong, Very nice, thanks! Barry > On Feb 27, 2017, at 2:06 PM, Hong <hzhang(a)mcs.anl.gov> wrote: > > Mark, > I implemented scalable MatPtAP and did comparisons of three implementations using ex56.c on alcf cetus machine (this machine has small memory, 1GB/core): > - nonscalable PtAP: use an array of length PN to do dense axpy > - scalable PtAP: do sparse axpy without use of PN array > - hypre PtAP. > > The results are attached. Summary: > - nonscalable PtAP is 2x faster than scalable, 8x faster than hypre PtAP > - scalable PtAP is 4x faster than hypre PtAP > - hypre uses less memory (see job.ne399.n63.np1000.sh) > > Based on above observation, I set the default PtAP algorithm as 'nonscalable'. > When PN > local estimated nonzero of C=PtAP, then switch default to 'scalable'. > User can overwrite default. > > For the case of np=8000, ne=599 (see job.ne599.n500.np8000.sh), I get > MatPtAP 3.6224e+01 (nonscalable for small mats, scalable for larger ones) > scalable MatPtAP 4.6129e+01 > hypre 1.9389e+02 > > I'm merging this work to next, then to master soon. > > Hong > > On Wed, Feb 8, 2017 at 7:25 PM, Mark Adams <mfadams(a)lbl.gov> wrote: > Thanks Hong, that sounds great. > > I am weary of silent optimizations like you suggest but 2x is big! and failure is very bad. So I would vote for your suggestion. > > ex56 is an elasticity problem. It would be nice, now that you have this experimental setup in place, to compare with hypre on a 3D scalar problem. Hypre might not spend much effort optimizing for block matrices. 3x better than hypre seems large to me. I have to suspect that hypre does not exploit blocked matrices as much as we do. > > Thanks again, > > > > On Wed, Feb 8, 2017 at 12:07 PM, Hong <hzhang(a)mcs.anl.gov> wrote: > I conducted tests on MatMatMult() and MatPtAP() using petsc/src/ksp/ksp/examples/tutorials/ex56.c (gamg) on a 8-core machine (petsc machine). The output file is attached. > > Summary: > 1) non-scalable MatMatMult() for mpiaij format is 2x faster than scalable version. The major difference between the two is dense-axpy vs. sparse-axpy. > > Currently, we set non-scalable as default, which leads to problem when running large problems. > How about setting default as > - non-scalable for small to medium size matrices > - scalable for larger ones, e.g. > > + ierr = PetscOptionsEList("-matmatmult_via","Algorithmic approach","MatMatMult",algTypes,nalg,algTypes[1],&alg,&flg); > > + if (!flg) { /* set default algorithm based on B->cmap->N */ > + PetscMPIInt size; > + ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); > + if ((PetscReal)(B->cmap->N)/size > 100000.0) alg = 0; /* scalable algorithm */ > + } > > i.e., if user does NOT pick an algorithm, when ave cols per process > 100k, use scalable implementation; otherwise, non-scalable version. > > 2) We do NOT have scalable implementation for MatPtAP() yet. > We have non-scalable PtAP and interface to Hypre's PtAP. Comparing the two, > Petsc MatPtAP() is approx 3x faster than Hypre's. > > I'm writing a scalable MatPtAP() now. > > Hong > > On Thu, Feb 2, 2017 at 2:54 PM, Stefano Zampini <stefano.zampini(a)gmail.com> wrote: > > > Il 02 Feb 2017 23:43, "Mark Adams" <mfadams(a)lbl.gov> ha scritto: > > > On Thu, Feb 2, 2017 at 12:02 PM, Stefano Zampini <stefano.zampini(a)gmail.com> wrote: > Mark, > > I saw your configuration has hypre. If you could run with master, you may try -matptap_via hypre. > > This is worth trying. Does this even work with GAMG? > > Yes, it should work, except that the block sizes, if any, are not propagated to the resulting matrix. I can add it if you need it. > > > > Treb: try hypre anyway. It has its own RAP code. > > > With that option, you will use hypre's RAP with MATAIJ > > > It uses BoomerAMGBuildCoarseOperator directly with the AIJ matrices. > > Stefano > >> On Feb 2, 2017, at 7:28 PM, Mark Adams <mfadams(a)lbl.gov> wrote: >> >> >> >> On Thu, Feb 2, 2017 at 11:13 AM, Hong <hzhang(a)mcs.anl.gov> wrote: >> Mark: >> Try '-matmatmult_via scalable' first. If this works, should we set it as default? >> >> If it is robust I would say yes unless it is noticeably slower (say >20%) small scale problems. >> > > > > > > > <out_ex56_cetus_short>

1 0

Re: [petsc-dev] no petsc on Edison
by Hong 27 Feb '17

27 Feb '17

Mark, I implemented scalable MatPtAP and did comparisons of three implementations using ex56.c on alcf cetus machine (this machine has small memory, 1GB/core): - nonscalable PtAP: use an array of length PN to do dense axpy - scalable PtAP: do sparse axpy without use of PN array - hypre PtAP. The results are attached. Summary: - nonscalable PtAP is 2x faster than scalable, 8x faster than hypre PtAP - scalable PtAP is 4x faster than hypre PtAP - hypre uses less memory (see job.ne399.n63.np1000.sh) Based on above observation, I set the default PtAP algorithm as 'nonscalable'. When PN > local estimated nonzero of C=PtAP, then switch default to 'scalable'. User can overwrite default. For the case of np=8000, ne=599 (see job.ne599.n500.np8000.sh), I get MatPtAP 3.6224e+01 (nonscalable for small mats, scalable for larger ones) scalable MatPtAP 4.6129e+01 hypre 1.9389e+02 I'm merging this work to next, then to master soon. Hong On Wed, Feb 8, 2017 at 7:25 PM, Mark Adams <mfadams(a)lbl.gov> wrote: > Thanks Hong, that sounds great. > > I am weary of silent optimizations like you suggest but 2x is big! and > failure is very bad. So I would vote for your suggestion. > > ex56 is an elasticity problem. It would be nice, now that you have this > experimental setup in place, to compare with hypre on a 3D scalar > problem. Hypre might not spend much effort optimizing for block matrices. > 3x better than hypre seems large to me. I have to suspect that hypre does > not exploit blocked matrices as much as we do. > > Thanks again, > > > > On Wed, Feb 8, 2017 at 12:07 PM, Hong <hzhang(a)mcs.anl.gov> wrote: > >> I conducted tests on MatMatMult() and MatPtAP() >> using petsc/src/ksp/ksp/examples/tutorials/ex56.c (gamg) on a 8-core >> machine (petsc machine). The output file is attached. >> >> Summary: >> 1) non-scalable MatMatMult() for mpiaij format is 2x faster than scalable >> version. The major difference between the two is dense-axpy vs. sparse-axpy. >> >> Currently, we set non-scalable as default, which leads to problem when >> running large problems. >> How about setting default as >> - non-scalable for small to medium size matrices >> - scalable for larger ones, e.g. >> >> + ierr = PetscOptionsEList("-matmatmult_via","Algorithmic >> approach","MatMatMult",algTypes,nalg,algTypes[1],&alg,&flg); >> >> + if (!flg) { /* set default algorithm based on B->cmap->N */ >> + PetscMPIInt size; >> + ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); >> + if ((PetscReal)(B->cmap->N)/size > 100000.0) alg = 0; /* scalable >> algorithm */ >> + } >> >> i.e., if user does NOT pick an algorithm, when ave cols per process > >> 100k, use scalable implementation; otherwise, non-scalable version. >> >> 2) We do NOT have scalable implementation for MatPtAP() yet. >> We have non-scalable PtAP and interface to Hypre's PtAP. Comparing the >> two, >> Petsc MatPtAP() is approx 3x faster than Hypre's. >> >> I'm writing a scalable MatPtAP() now. >> >> Hong >> >> On Thu, Feb 2, 2017 at 2:54 PM, Stefano Zampini < >> stefano.zampini(a)gmail.com> wrote: >> >>> >>> >>> Il 02 Feb 2017 23:43, "Mark Adams" <mfadams(a)lbl.gov> ha scritto: >>> >>> >>> >>> On Thu, Feb 2, 2017 at 12:02 PM, Stefano Zampini < >>> stefano.zampini(a)gmail.com> wrote: >>> >>>> Mark, >>>> >>>> I saw your configuration has hypre. If you could run with master, you >>>> may try -matptap_via hypre. >>>> >>> >>> This is worth trying. Does this even work with GAMG? >>> >>> >>> Yes, it should work, except that the block sizes, if any, are not >>> propagated to the resulting matrix. I can add it if you need it. >>> >>> >>> >>> Treb: try hypre anyway. It has its own RAP code. >>> >>> >>> >>> With that option, you will use hypre's RAP with MATAIJ >>> >>> >>> It uses BoomerAMGBuildCoarseOperator directly with the AIJ matrices. >>>> >>>> Stefano >>>> >>>> On Feb 2, 2017, at 7:28 PM, Mark Adams <mfadams(a)lbl.gov> wrote: >>>> >>>> >>>> >>>> On Thu, Feb 2, 2017 at 11:13 AM, Hong <hzhang(a)mcs.anl.gov> wrote: >>>> >>>>> Mark: >>>>> Try '-matmatmult_via scalable' first. If this works, should we set it >>>>> as default? >>>>> >>>> >>>> If it is robust I would say yes unless it is noticeably slower (say >>>> >20%) small scale problems. >>>> >>>> >>>> >>> >>> >> >

1 0