MatMult_Scatter woes

[Barry Smith]

On Mar 30, 2009, at 7:48 PM, Matthew Knepley wrote:

> 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.

  No, when two map to the same output the scatter is not well defined  
and hence the multiply is not defined.
When building the scatter we could check for multiple inputs going to  
the same output but we don't, it will end up
just putting one value in and then the other and depending on the  
order of the message passing it could be different.


    Barry

> 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@petsc.cs.iit.edu/petsc/petsc-dev without success.
>
> It is ssh://petsc@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