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

Re: [petsc-dev] GASM
by Dmitry Karpeev 29 Mar '12

29 Mar '12

On Thu, Mar 29, 2012 at 8:28 AM, Jed Brown <jedbrown(a)mcs.anl.gov> wrote: > On Thu, Mar 29, 2012 at 07:16, Dmitry Karpeev <karpeev(a)mcs.anl.gov> wrote: > >> Here the problem isn't in needing to call PCReset() twice per se. The >> issue is that if the user uses the same PC object in >> two different KSPs and destroys one ksp, that will cause a PCReset() and >> invalidate the pc for the other ksp -- all despite increfing the pc. >> > > Yeah, this is tricky. In this case, is it acceptable to break the > dependency explicitly using > > PCSetType(ksp,PCNONE); // will decrement the reference count, but not > PCReset() > KSPDestroy(&ksp); > > > Skipping the reset means that KSPReset() is basically always incorrect if > the user created the PC (for example). I think a top-level KSPReset() > should always reset the PC. > I am not sure how to solve the problem of the relationship between KSPReset() and PCReset() when the pc is being shared. How useful are KSPReset() and PCReset() anyway? All but one uses I found were "trivial" -- called from a PCReset_XXX, where XXX uses a KSP internally. The only nontrivial use I found was in VIRS, where the KSP and PC are reset when the active set changes. Would a simple KSPDestroy() do here instead? Perhaps there are uses of KSPReset()/PCReset() in user code, though. > Now maybe ksp->pc should be withheld by KSPDestroy() before calling > KSPReset(); if the PC refct hits 0, it will be destroyed in the subsequent > PCDestroy(). > The present problem seems to come from a misuse of KSPReset()/PCDestroy() as a partial destructor inside KSPDestroy(): the intent here is to destroy the ksp and the underlying pc, and in the case of the pc this should be protected by the refcount. However, PCDestroy() is being entered "from a side" for the intended meaning -- via KSPReset(), not PCDestroy() and the refcount check isn't being triggered. One solution would be not to use KSPReset()/PCReset() as a partial destructor, but that would require duplicating that code directly in KSPDestroy()/PCDestroy(), or we could move this code to KSPReset_Private() that would not cascade to PCReset_Private(). Then we could have: PetscErrorCode KSPReset(KSP ksp) { ... KSPReset_Private(ksp); PCReset_Private(ksp->pc); ... } PetscErrorCode KSPDestroy(KSP *ksp) { ... KSPReset_Private(*ksp); PCDestroy((*ksp)->pc); ... } PetscErrorCode PCDestroy(PC *pc) { if (--((PetscObject)(*pc))->refct > 0) {*pc = 0; PetscFunctionReturn(0);} ... } Dmitry.

1 0

Re: [petsc-dev] GASM
by Jed Brown 29 Mar '12

29 Mar '12

On Thu, Mar 29, 2012 at 07:16, Dmitry Karpeev <karpeev(a)mcs.anl.gov> wrote: > Here the problem isn't in needing to call PCReset() twice per se. The > issue is that if the user uses the same PC object in > two different KSPs and destroys one ksp, that will cause a PCReset() and > invalidate the pc for the other ksp -- all despite increfing the pc. > Yeah, this is tricky. In this case, is it acceptable to break the dependency explicitly using PCSetType(ksp,PCNONE); // will decrement the reference count, but not PCReset() KSPDestroy(&ksp); Skipping the reset means that KSPReset() is basically always incorrect if the user created the PC (for example). I think a top-level KSPReset() should always reset the PC. Now maybe ksp->pc should be withheld by KSPDestroy() before calling KSPReset(); if the PC refct hits 0, it will be destroyed in the subsequent PCDestroy().

1 0

Re: [petsc-dev] GASM
by Dmitry Karpeev 29 Mar '12

29 Mar '12

On Thu, Mar 29, 2012 at 8:07 AM, Jed Brown <jedbrown(a)mcs.anl.gov> wrote: > On Thu, Mar 29, 2012 at 06:44, Dmitry Karpeev <karpeev(a)mcs.anl.gov> wrote: > >> After this fix that should cause a more obvious SEGV, but the real >> problem is that a pc that is supposed remain valid is being reset. > > > It's a bug if you can't call PCReset() more than once without an error. > Here the problem isn't in needing to call PCReset() twice per se. The issue is that if the user uses the same PC object in two different KSPs and destroys one ksp, that will cause a PCReset() and invalidate the pc for the other ksp -- all despite increfing the pc. > Also, any SEGV triggered by PETSc is a bug (unless the user passes in > corrupt pointers). > I was a bit sloppy in describing what will happen here: now the (G)ASM code actually won't throw a SEGV, but will rather rebuild ISs internally using a default partition -- PCReset_ASM/PCReset_GASM does not forget the *number* of subdomains upon a Reset, but destroys the actual subdomains. Dmitry.

1 0

Re: [petsc-dev] GASM
by Jed Brown 29 Mar '12

29 Mar '12

On Thu, Mar 29, 2012 at 06:44, Dmitry Karpeev <karpeev(a)mcs.anl.gov> wrote: > After this fix that should cause a more obvious SEGV, but the real > problem is that a pc that is supposed remain valid is being reset. It's a bug if you can't call PCReset() more than once without an error. Also, any SEGV triggered by PETSc is a bug (unless the user passes in corrupt pointers).

1 0

Re: [petsc-dev] GASM
by Dmitry Karpeev 29 Mar '12

29 Mar '12

Mark, I pushed a fix for that specific bug, but from my reading of the GAMG code this won't solve the real problem: PCSetUp_GAMG() reuses pc as a subpc in eksp, protecting subpc with an incref. That prevents the subpc from being destroyed when eksp is destroyed, but the somewhat unfortunate logic of KSPDestroy(eksp) also calls KSPReset(eksp), which cascades to PCReset(subpc), blowing out its internals -- PCReset() pays no mind to the refcount. The bug was that GASM's guts weren't being destroyed very cleanly, allowing a later read into a freed array. After this fix that should cause a more obvious SEGV, but the real problem is that a pc that is supposed remain valid is being reset. I'm not sure what the best way to address the general problem should be: should KSPDestroy() has to be redesigned? should PCReset() be a noop for a positive refcount? or should this use of a pc in two places be deemed illegal? Dmitry. On Wed, Mar 28, 2012 at 8:24 PM, Mark F. Adams <mark.adams(a)columbia.edu>wrote: > Dmitry, > > I seem to have a problem with GASM. I have code like this to get an eigen > estimate: > > ierr = PetscObjectReference( (PetscObject)subpc ); > CHKERRQ(ierr); > ierr = KSPSetPC( eksp, subpc ); CHKERRQ( ierr ); > < snip > > ierr = KSPDestroy( &eksp ); CHKERRQ(ierr); > > This 'subpc' has GASM and is used in the real solver. And I get this > error when the real solver gets setup, after the above code: > > ==2070== Invalid read of size 4 > ==2070== at 0x9F034C: PCSetUp_GASM (gasm.c:275) > ==2070== by 0xA9620A: PCSetUp (precon.c:832) > ==2070== by 0xB20EE5: KSPSetUp (itfunc.c:266) > ==2070== by 0xA74975: PCSetUp_MG (mg.c:647) > ==2070== by 0xA11A1B: PCSetUp_GAMG (gamg.c:803) > ==2070== by 0xA9620A: PCSetUp (precon.c:832) > ==2070== by 0xB20EE5: KSPSetUp (itfunc.c:266) > ==2070== by 0xB22178: KSPSolve (itfunc.c:390) > ==2070== by 0x95B8C0: kspsolve_ (itfuncf.c:219) > ==2070== by 0x404E: MAIN__ (in ./ex54f) > ==2070== by 0x4A6B: main (in ./ex54f) > ==2070== Address 0x29e6ac0 is 32 bytes inside a block of size 108 free'd > ==2070== at 0x19317D8: free (vg_replace_malloc.c:366) > ==2070== by 0x865AD: PetscFreeAlign (mal.c:75) > ==2070== by 0x9FAB7F: PCGASMDestroySubdomains (gasm.c:1309) > ==2070== by 0x9F3AED: PCReset_GASM (gasm.c:536) > ==2070== by 0xA8B322: PCReset (precon.c:80) > ==2070== by 0xB270D0: KSPReset (itfunc.c:723) > ==2070== by 0xB279BA: KSPDestroy (itfunc.c:760) > ==2070== by 0xA1164D: PCSetUp_GAMG (gamg.c:779) > ==2070== by 0xA9620A: PCSetUp (precon.c:832) > ==2070== by 0xB20EE5: KSPSetUp (itfunc.c:266) > ==2070== by 0xB22178: KSPSolve (itfunc.c:390) > ==2070== by 0x95B8C0: kspsolve_ (itfuncf.c:219) > ==2070== by 0x404E: MAIN__ (in ./ex54f) > ==2070== by 0x4A6B: main (in ./ex54f) > ==2070== > [0]PETSC ERROR: --------------------- Error Message > ------------------------------------ > [0]PETSC ERROR: Null argument, when expecting valid pointer! > [0]PETSC ERROR: Null Object: Parameter # 1! > [0]PETSC ERROR: > ------------------------------------------------------------------------ > [0]PETSC ERROR: Petsc Development HG revision: > e4ca525c39fd35cca849d095040a2de0488f347b HG Date: Fri Mar 23 22:45:44 2012 > -0500 > [0]PETSC ERROR: See docs/changes/index.html for recent updates. > [0]PETSC ERROR: See docs/faq.html for hints about trouble shooting. > [0]PETSC ERROR: See docs/index.html for manual pages. > [0]PETSC ERROR: > ------------------------------------------------------------------------ > [0]PETSC ERROR: ./ex54f on a arch-maco named madams-macbk-3.home by > markadams Wed Mar 28 21:20:16 2012 > [0]PETSC ERROR: Libraries linked from > /Users/markadams/Codes/petsc-dev/arch-macosx-gnu-O/lib > [0]PETSC ERROR: Configure run at Sat Mar 24 09:45:16 2012 > [0]PETSC ERROR: Configure options CXX="/sw/lib/gcc4.6/bin/g++-4 > -malign-double" CC="/sw/lib/gcc4.6/bin/gcc-4 -malign-double" > FC=/sw/lib/gcc4.6/bin/gfortran --download-parmetis=1 --download-metis=1 > --download-hdf5=1 --download-hypre=1 --download-ml=1 --download-triangle=1 > --with-x=0 --with-debugging=1 --download-mpich=1 > PETSC_ARCH=arch-macosx-gnu-O --download-f-blas-lapack=1 > [0]PETSC ERROR: > ------------------------------------------------------------------------ > [0]PETSC ERROR: ISSort() line 761 in > /Users/markadams/Codes/petsc-dev/src/vec/is/interface/index.c > [0]PETSC ERROR: PCSetUp_GASM() line 275 in > /Users/markadams/Codes/petsc-dev/src/ksp/pc/impls/gasm/gasm.c > > THis error goes away if I comment out the above call to KSPDestroy. So it > looks like KSPDestroy is destroying GASM (or part of it) even though I take > a reference to it, to prevent this destruction. > > Any ideas? > > Mark > > > Mark, > > You should pass your subdomain ISs through the first IS list argument, and > leave the second as PETSC_NULL. > In that case GASM assumes that the user passed in a nonoverlapping > decomposition through the first list. > This is the behavior that GASM inherited from ASM and it isn't very well > (or at all) documented in either. > Let me know if that doesn't fix the problem. > > Note also that you don't have to explicitly set overlap to 0: if the user > sets not overlap PETSc won't try to increase it, > and if you set overlap to 0 it will not try to adjust the subdomains to > make sure the overlap is 0. Both ASM and GASM simply assume that the user > does the right thing and passes in a nonoverlapping decomposition. > > Incidentally, it may be time to change this API. I'm trying to debug a > bad convergence problem with GASM -- inexplicably, it converges slower than > ASM with a single process and multiple subdomains -- so now is a good time > to mess with the code, while I still remember what's what. Another thing > I've considered doing is reordering the subdomains internally (on demand > only), to ensure a deadlock-free ordering of the subcomms. > > Dmitry. > > On Tue, Mar 27, 2012 at 6:28 PM, Barry Smith <bsmith(a)mcs.anl.gov> wrote: > >> >> This is why petsc-dev is better than a random list of email addresses. >> No one on your random list uses GASM >> >> >> >> Barry >> >> On Mar 27, 2012, at 4:49 PM, Mark F. Adams wrote: >> >> > I'm trying to make general non-overlapping domains in GASM. I create a >> vector of IS, each with the global indices of the equations in each block. >> > >> > I use code that looks like: >> > >> > ierr = PCSetType( subpc, PCGASM ); CHKERRQ(ierr); >> > ierr = PCGASMSetLocalSubdomains( subpc, sz, PETSC_NULL, is ); >> CHKERRQ(ierr); >> > ierr = PetscFree( is ); CHKERRQ(ierr); >> > ierr = PCGASMSetOverlap( subpc, 0 ); CHKERRQ(ierr); >> > >> > In parallel, I get errors like this: >> > >> > [1]PETSC ERROR: gis_local contains -85 indices outside of gis! >> > <snip> >> > [1]PETSC ERROR: PCSetUp_GASM() line 332 in >> /Users/markadams/Codes/petsc-dev/src/ksp/pc/impls/gasm/gasm.c >> > >> > Anyone have any ideas on what I'm doing wrong here? I've checked the >> input and I think it is fine, >> > >> > Mark >> >> > >

1 0

Re: [petsc-dev] How does this compare...
by Jack Poulson 28 Mar '12

28 Mar '12

I dug in a little deeper, and I recommend reading their previous paper: http://arxiv.org/abs/1011.3534 That paper actually _does_ mention the connection to H-matrix algebra: "SpAMM is similar to the H-matrix algebra of Hackbusch and co-workers, where off-diagonal sub-matrices are treated as reduced rank factorizations (truncated SVD), typically structured and grouped to reflect properties of the underlying operators. For problems with rapid decay, truncated SVD may behave in a similar way to simple dropping schemes. SpAMM is different than the H-matrix algebra as it achieves separation uniquely in the product space and does not rely on a reduced complexity representation of matrices. For very slow decay, the H-matrix algebra may certainly offer a path for intractable problems for which SpAMM is ineffective." So it seems that the claim is that SpAMM might be more computationally efficient than H-matrices for fast decay, but ineffective otherwise. Either way, I don't think it's going to help with classical parallel sparse matrix-matrix multiplication. Jack On Wed, Mar 28, 2012 at 10:28 PM, Jack Poulson <jack.poulson(a)gmail.com>wrote: > Actually, please excuse my misreading. The quoted 2-3x speedup was versus > itself due to hardware prefetch, so it does appear to be asymptotically > better than SGEMM (they quote the crossover at n ~= 1000). > > Jack > > > On Wed, Mar 28, 2012 at 10:23 PM, Jack Poulson <jack.poulson(a)gmail.com>wrote: > >> I just took a look at the paper and, after reading a bit, I am baffled >> that there are no citations of the H-matrix literature (Ctrl+F "Hackbusch" >> turns up nothing...). The paper seems to be in a similar spirit as the >> well-known H-matrix/H-matrix multiplication algorithms, which take O(N >> log^2 N) work, for O(1) ranks, but using truncation instead of low-rank >> approximations in clusters. Since only factors of 2-3x are reported versus >> SGEMM, which is O(N^3), that tells you that this algorithm (or, at least, >> its implementation) is O(N^3). >> >> I do not think a comparison can be drawn to sparse matrix-matrix multiply >> algorithms. I believe Gilbert and Buluc have the most recent publication on >> parallel sparse mat-mat. I have done some work on parallel >> H-matrix/H-matrix multiplication (the code is here: >> http://bitbucket.org/poulson/dmhm) but I haven't pushed the paper out >> the door yet. See http://www.ices.utexas.edu/~poulson/slides/DMHM.pdffor the talk I gave on it. I suspect that the ideas from parallel >> H-matrix/H-matrix multiplication could be used to generate a >> communication-efficient sparse mat-mat if it is okay to generate the result >> in a more general data-sparse form (i.e., low rank blocks) rather than >> specifically as explicitly classically sparse. >> >> Jack >> >> >> On Wed, Mar 28, 2012 at 10:07 PM, Hong Zhang <hzhang(a)mcs.anl.gov> wrote: >> >>> Entire paper discusses SpAMM, and compares with SGEMM. >>> Sorry, I do not see anything about SpMV. >>> Hong >>> >>> On Wed, Mar 28, 2012 at 9:52 PM, Matthew Knepley <knepley(a)gmail.com>wrote: >>> >>>> On Wed, Mar 28, 2012 at 9:39 PM, Hong Zhang <hzhang(a)mcs.anl.gov> wrote: >>>> >>>>> Matt: >>>>> >>>>>> to the new scalable MatMatMult() that is done in PETSc? Is it >>>>>> possible to even compare >>>>>> the flop rates sensibly? >>>>>> >>>>>> http://arxiv.org/abs/1203.1692 >>>>>> >>>>> "An Optimized Sparse Approximate Matrix Multiply" >>>>> computes Approximate matrix product for matrices with decay >>>>> ^^^^^^^^^^^^ >>>>> ^^^^^^^^^^ >>>>> It is for special application in chemistry(?), >>>>> not applicable for general sparse matrix product. >>>>> No comparison can be drawn to petsc's C=A*B. >>>>> >>>> >>>> This is the wrong conclusion. What he does is sparsify the result, so >>>> he is >>>> computing SOME sort of sparse matrix-vector product. The question is: >>>> >>>> Does his SpMV use the same algorithm as PETSc? >>>> >>>> That is why I asked. >>>> >>>> Matt >>>> >>>> >>>>> Hong >>>>> >>>>>> >>>>>> This is what Jeff is planning on using in his SciDAC. >>>>>> >>>>>> 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 >>>>>> >>>>> >>>>> >>>> >>>> >>>> -- >>>> 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] How does this compare...
by Jack Poulson 28 Mar '12

28 Mar '12

Actually, please excuse my misreading. The quoted 2-3x speedup was versus itself due to hardware prefetch, so it does appear to be asymptotically better than SGEMM (they quote the crossover at n ~= 1000). Jack On Wed, Mar 28, 2012 at 10:23 PM, Jack Poulson <jack.poulson(a)gmail.com>wrote: > I just took a look at the paper and, after reading a bit, I am baffled > that there are no citations of the H-matrix literature (Ctrl+F "Hackbusch" > turns up nothing...). The paper seems to be in a similar spirit as the > well-known H-matrix/H-matrix multiplication algorithms, which take O(N > log^2 N) work, for O(1) ranks, but using truncation instead of low-rank > approximations in clusters. Since only factors of 2-3x are reported versus > SGEMM, which is O(N^3), that tells you that this algorithm (or, at least, > its implementation) is O(N^3). > > I do not think a comparison can be drawn to sparse matrix-matrix multiply > algorithms. I believe Gilbert and Buluc have the most recent publication on > parallel sparse mat-mat. I have done some work on parallel > H-matrix/H-matrix multiplication (the code is here: > http://bitbucket.org/poulson/dmhm) but I haven't pushed the paper out > the door yet. See http://www.ices.utexas.edu/~poulson/slides/DMHM.pdf for > the talk I gave on it. I suspect that the ideas from parallel > H-matrix/H-matrix multiplication could be used to generate a > communication-efficient sparse mat-mat if it is okay to generate the result > in a more general data-sparse form (i.e., low rank blocks) rather than > specifically as explicitly classically sparse. > > Jack > > > On Wed, Mar 28, 2012 at 10:07 PM, Hong Zhang <hzhang(a)mcs.anl.gov> wrote: > >> Entire paper discusses SpAMM, and compares with SGEMM. >> Sorry, I do not see anything about SpMV. >> Hong >> >> On Wed, Mar 28, 2012 at 9:52 PM, Matthew Knepley <knepley(a)gmail.com>wrote: >> >>> On Wed, Mar 28, 2012 at 9:39 PM, Hong Zhang <hzhang(a)mcs.anl.gov> wrote: >>> >>>> Matt: >>>> >>>>> to the new scalable MatMatMult() that is done in PETSc? Is it possible >>>>> to even compare >>>>> the flop rates sensibly? >>>>> >>>>> http://arxiv.org/abs/1203.1692 >>>>> >>>> "An Optimized Sparse Approximate Matrix Multiply" >>>> computes Approximate matrix product for matrices with decay >>>> ^^^^^^^^^^^^ >>>> ^^^^^^^^^^ >>>> It is for special application in chemistry(?), >>>> not applicable for general sparse matrix product. >>>> No comparison can be drawn to petsc's C=A*B. >>>> >>> >>> This is the wrong conclusion. What he does is sparsify the result, so he >>> is >>> computing SOME sort of sparse matrix-vector product. The question is: >>> >>> Does his SpMV use the same algorithm as PETSc? >>> >>> That is why I asked. >>> >>> Matt >>> >>> >>>> Hong >>>> >>>>> >>>>> This is what Jeff is planning on using in his SciDAC. >>>>> >>>>> 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 >>>>> >>>> >>>> >>> >>> >>> -- >>> 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] [petsc-users] Bug in pc_right + nonzero initial guess + {tfqmr, tcqmr, cgs}
by Barry Smith 28 Mar '12

28 Mar '12

From KSPSolve_BCGS() /* Compute initial preconditioned residual */ ierr = KSPInitialResidual(ksp,X,V,T,R,B);CHKERRQ(ierr); /* with right preconditioning need to save initial guess to add to final solution */ if (ksp->pc_side == PC_RIGHT && !ksp->guess_zero) { if (!bcgs->guess) { ierr = VecDuplicate(X,&bcgs->guess);CHKERRQ(ierr); } ierr = VecCopy(X,bcgs->guess);CHKERRQ(ierr); ierr = VecSet(X,0.0);CHKERRQ(ierr); } Thus in the right preconditioned, non-zero initial guess case it never computes the original norm of b so the convergence test is "different"; for example if the initial guess is really really good the convergence test doesn't know this so will iterate normally more than one would likely want. We can fix this by computing the norm of b (when using the default convergence test) at the beginning passing it to the convergence test and then pass the || r|| norm to the convergence test but then that is like iteration 1/2.. Do we need to change the handling of the default convergence test to "eliminate" this 1/2 iteration? In answer to how to fix the other methods: I'm thinking maybe this storage of the initial guess, computing first residual, zeroing of solution and addition back on at the end of the solve should be embedded in the KSPInitialResidual() and KSPUnwindPreconditioner() routines. And it could be done perhaps in a more systematic way; for example GMRES does not use the paradigm of bcgs (and the others when we fix them) of simply removing the non-zero guess and changing the right hand side; should we change it? Barry On Mar 13, 2012, at 7:47 AM, Jed Brown wrote: > On Sun, Mar 11, 2012 at 00:38, Tobin Isaac <tisaac(a)ices.utexas.edu> wrote: > > Hi, > > The changes made to bcgs in the commit below should be applied to > tfqmr, tcqmr, and cgs as well. > > Toby, thanks for hunting this down. Barry, do you think we can reasonably share code between these implementations (and any others that might need it)? > > > Cheers, > Toby > > author Barry Smith bsmith(a)mcs.anl.gov > Wed, 19 Jan 2011 20:10:29 -0600 (13 months ago) > changeset 18159 a526acdfdc95 > parent 18153 3c53022524fb > child 18160 4db6569c26af > > fixed KSPBCGS to work with right preconditioning and nonzero initial > guess > > -----BEGIN PGP SIGNATURE----- > Version: GnuPG v1.4.10 (GNU/Linux) > > iEYEARECAAYFAk9cSIAACgkQk/TrNolnueUqMQCcD9gSdldhOY14sIFpwEC4nhpD > pTEAniP556Lh3+5MgYB9cTUzHj7I4jg5 > =FvoR > -----END PGP SIGNATURE----- > >

1 0

Re: [petsc-dev] How does this compare...
by Jack Poulson 28 Mar '12

28 Mar '12

I just took a look at the paper and, after reading a bit, I am baffled that there are no citations of the H-matrix literature (Ctrl+F "Hackbusch" turns up nothing...). The paper seems to be in a similar spirit as the well-known H-matrix/H-matrix multiplication algorithms, which take O(N log^2 N) work, for O(1) ranks, but using truncation instead of low-rank approximations in clusters. Since only factors of 2-3x are reported versus SGEMM, which is O(N^3), that tells you that this algorithm (or, at least, its implementation) is O(N^3). I do not think a comparison can be drawn to sparse matrix-matrix multiply algorithms. I believe Gilbert and Buluc have the most recent publication on parallel sparse mat-mat. I have done some work on parallel H-matrix/H-matrix multiplication (the code is here: http://bitbucket.org/poulson/dmhm) but I haven't pushed the paper out the door yet. See http://www.ices.utexas.edu/~poulson/slides/DMHM.pdf for the talk I gave on it. I suspect that the ideas from parallel H-matrix/H-matrix multiplication could be used to generate a communication-efficient sparse mat-mat if it is okay to generate the result in a more general data-sparse form (i.e., low rank blocks) rather than specifically as explicitly classically sparse. Jack On Wed, Mar 28, 2012 at 10:07 PM, Hong Zhang <hzhang(a)mcs.anl.gov> wrote: > Entire paper discusses SpAMM, and compares with SGEMM. > Sorry, I do not see anything about SpMV. > Hong > > On Wed, Mar 28, 2012 at 9:52 PM, Matthew Knepley <knepley(a)gmail.com>wrote: > >> On Wed, Mar 28, 2012 at 9:39 PM, Hong Zhang <hzhang(a)mcs.anl.gov> wrote: >> >>> Matt: >>> >>>> to the new scalable MatMatMult() that is done in PETSc? Is it possible >>>> to even compare >>>> the flop rates sensibly? >>>> >>>> http://arxiv.org/abs/1203.1692 >>>> >>> "An Optimized Sparse Approximate Matrix Multiply" >>> computes Approximate matrix product for matrices with decay >>> ^^^^^^^^^^^^ >>> ^^^^^^^^^^ >>> It is for special application in chemistry(?), >>> not applicable for general sparse matrix product. >>> No comparison can be drawn to petsc's C=A*B. >>> >> >> This is the wrong conclusion. What he does is sparsify the result, so he >> is >> computing SOME sort of sparse matrix-vector product. The question is: >> >> Does his SpMV use the same algorithm as PETSc? >> >> That is why I asked. >> >> Matt >> >> >>> Hong >>> >>>> >>>> This is what Jeff is planning on using in his SciDAC. >>>> >>>> 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 >>>> >>> >>> >> >> >> -- >> 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] How does this compare...
by Hong Zhang 28 Mar '12

28 Mar '12

Entire paper discusses SpAMM, and compares with SGEMM. Sorry, I do not see anything about SpMV. Hong On Wed, Mar 28, 2012 at 9:52 PM, Matthew Knepley <knepley(a)gmail.com> wrote: > On Wed, Mar 28, 2012 at 9:39 PM, Hong Zhang <hzhang(a)mcs.anl.gov> wrote: > >> Matt: >> >>> to the new scalable MatMatMult() that is done in PETSc? Is it possible >>> to even compare >>> the flop rates sensibly? >>> >>> http://arxiv.org/abs/1203.1692 >>> >> "An Optimized Sparse Approximate Matrix Multiply" >> computes Approximate matrix product for matrices with decay >> ^^^^^^^^^^^^ >> ^^^^^^^^^^ >> It is for special application in chemistry(?), >> not applicable for general sparse matrix product. >> No comparison can be drawn to petsc's C=A*B. >> > > This is the wrong conclusion. What he does is sparsify the result, so he is > computing SOME sort of sparse matrix-vector product. The question is: > > Does his SpMV use the same algorithm as PETSc? > > That is why I asked. > > Matt > > >> Hong >> >>> >>> This is what Jeff is planning on using in his SciDAC. >>> >>> 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 >>> >> >> > > > -- > 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