On Mon, Mar 30, 2009 at 6:13 PM, Barry Smith <bsmith at mcs.anl.gov> wrote:
> > I don't think square has anything to do with it. It depends on whether > all the entries in y[i] are hit or if some are skipped, this depends exactly > on the scatter. > > The first time the mult is called, duplicate x and y temporarility, fill > a xdup with all 1s. Do the scatter, locate all entries of ydup that are > zero, keep a list of these. > destroy xdup, ydup. > > Now do the MatMult() by zeroing all the entires in the list and then do > the scatter, now you have an efficient and correct MatMult. The drawback of > course is the setup time, but in fact that cost is only one > scatter (MatMult()) so is likely worth it. > > One could do a similar setup for the transpose? I must misunderstand you. The problem above is that INSERT_VALUES is used, but when two input locations map to the same output location, they should be added. You are right, that square matrices can still have this property. However, I do not see how the above sequence fixes this. I can't see how you would get around an addition. Matt > > Barry > > > On Mar 30, 2009, at 7:04 PM, Matthew Knepley wrote: > > On Mon, Mar 30, 2009 at 5:52 AM, Jed Brown <jed at 59a2.org> wrote: >> The current implementation of MatMult_Scatter and >> MatMultTranspose_Scatter do not actually have matrix semantics unless >> the index sets are permutations. For example, it's generally desirable >> that >> >> MatMult(A,x,y); >> >> agrees with >> >> VecZeroEntries(y); >> MatMultAdd(A,x,y,y); >> >> Similarly for MatMultTranspose. In addition, MatMultTranspose should >> actually be the transpose of MatMult. Consider the scatter defined by >> >> scatter : x -> y >> isx = [0 0] >> isy = [0 1] >> >> where x has length 1 and y has length 2. This forward scatter is >> equivalent to a 2x1 matrix A = [1;1] >> >> Indeed A*[1] = [1;1] >> >> but A'*[1;1] = [2] >> >> where as MatMultTranspose_Scatter(A,y=[1;1],x) gives x=[1]. >> >> This can be corrected by changing the body of MatMultTranspose_Scatter >> from >> >> ierr = >> VecScatterBegin(scatter->scatter,x,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); >> ierr = >> VecScatterEnd(scatter->scatter,x,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); >> >> to >> >> ierr = VecZeroEntries(y);CHKERRQ(ierr); >> ierr = >> VecScatterBegin(scatter->scatter,x,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); >> ierr = >> VecScatterEnd(scatter->scatter,x,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); >> >> and similarly for MatMult_Scatter. >> >> Unfortunately, these modifications roughly double the cost of typical >> sequential scatters and could be much worse (scattering from a small >> vector to a very large one). I think that correctness is more important >> here and users would typically not use MatScatter when this would have >> significant impact or when INSERT_VALUES semantics are desired. Of >> course some performance can be recovered by having MatScatter use >> INSERT_VALUES when the destination index set is a permutation. >> >> I would vote for checking for a square matrix, and otherwise use the >> expensive form. >> >> >> Also, what is the ssh url for the repo? I've tried some variations on >> ssh://petsc at petsc.cs.iit.edu/petsc/petsc-dev without success. >> >> It is ssh://petsc at petsc.cs.iit.edu//hg/petsc/petsc-dev >> >> Matt >> >> >> Jed >> -- >> 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 -------------- next part -------------- An HTML attachment was scrubbed... URL: <http://lists.mcs.anl.gov/pipermail/petsc-dev/attachments/20090330/ea22e9e4/attachment.html>