MatMult_Scatter woes

[Barry Smith]

    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?

    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