MatMult_Scatter woes

[Matthew Knepley]

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@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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-dev/attachments/20090330/ea22e9e4/attachment.html>