Yes, this looks great. Good to get some data on this. Hard problem. On Mon, Feb 27, 2017 at 3:31 PM, Barry Smith <[email protected]> wrote:
> > Hong, > > Very nice, thanks! > > Barry > > > On Feb 27, 2017, at 2:06 PM, Hong <[email protected]> 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 <[email protected]> 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 <[email protected]> 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 < > [email protected]> wrote: > > > > > > Il 02 Feb 2017 23:43, "Mark Adams" <[email protected]> ha scritto: > > > > > > On Thu, Feb 2, 2017 at 12:02 PM, Stefano Zampini < > [email protected]> 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 <[email protected]> wrote: > >> > >> > >> > >> On Thu, Feb 2, 2017 at 11:13 AM, Hong <[email protected]> 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> > >
