On Mon, Mar 30, 2009 at 6:13 PM, Barry Smith <span dir="ltr"><<a href="mailto:bsmith@mcs.anl.gov">bsmith@mcs.anl.gov</a>></span> wrote:<br><div class="gmail_quote"><blockquote class="gmail_quote" style="border-left: 1px solid rgb(204, 204, 204); margin: 0pt 0pt 0pt 0.8ex; padding-left: 1ex;">
<br>
   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.<br>
<br>
   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.<br>
destroy xdup, ydup.<br>
<br>
   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<br>
scatter (MatMult()) so is likely worth it.<br>
<br>
  One could do a similar setup for the transpose?<font color="#888888"></font></blockquote><div><br>I must misunderstand you. The problem above is that INSERT_VALUES is used, but when two input locations<br>map to the same output location, they should be added. You are right, that square matrices can still have this<br>
property. However, I do not see how the above sequence fixes this. I can't see how you would get around an<br>addition.<br><br>   Matt<br> </div><blockquote class="gmail_quote" style="border-left: 1px solid rgb(204, 204, 204); margin: 0pt 0pt 0pt 0.8ex; padding-left: 1ex;">
<font color="#888888"><br>
   Barry</font><div><div></div><div class="h5"><br>
<br>
On Mar 30, 2009, at 7:04 PM, Matthew Knepley wrote:<br>
<br>
<blockquote class="gmail_quote" style="border-left: 1px solid rgb(204, 204, 204); margin: 0pt 0pt 0pt 0.8ex; padding-left: 1ex;">
On Mon, Mar 30, 2009 at 5:52 AM, Jed Brown <<a href="mailto:jed@59a2.org" target="_blank">jed@59a2.org</a>> wrote:<br>
The current implementation of MatMult_Scatter and<br>
MatMultTranspose_Scatter do not actually have matrix semantics unless<br>
the index sets are permutations.  For example, it's generally desirable<br>
that<br>
<br>
 MatMult(A,x,y);<br>
<br>
agrees with<br>
<br>
 VecZeroEntries(y);<br>
 MatMultAdd(A,x,y,y);<br>
<br>
Similarly for MatMultTranspose.  In addition, MatMultTranspose should<br>
actually be the transpose of MatMult.  Consider the scatter defined by<br>
<br>
scatter : x -> y<br>
isx = [0 0]<br>
isy = [0 1]<br>
<br>
where x has length 1 and y has length 2.  This forward scatter is<br>
equivalent to a 2x1 matrix A = [1;1]<br>
<br>
Indeed A*[1] = [1;1]<br>
<br>
but A'*[1;1] = [2]<br>
<br>
where as MatMultTranspose_Scatter(A,y=[1;1],x) gives x=[1].<br>
<br>
This can be corrected by changing the body of MatMultTranspose_Scatter<br>
from<br>
<br>
 ierr = VecScatterBegin(scatter->scatter,x,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);<br>
 ierr = VecScatterEnd(scatter->scatter,x,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);<br>
<br>
to<br>
<br>
 ierr = VecZeroEntries(y);CHKERRQ(ierr);<br>
 ierr = VecScatterBegin(scatter->scatter,x,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);<br>
 ierr = VecScatterEnd(scatter->scatter,x,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);<br>
<br>
and similarly for MatMult_Scatter.<br>
<br>
Unfortunately, these modifications roughly double the cost of typical<br>
sequential scatters and could be much worse (scattering from a small<br>
vector to a very large one).  I think that correctness is more important<br>
here and users would typically not use MatScatter when this would have<br>
significant impact or when INSERT_VALUES semantics are desired.  Of<br>
course some performance can be recovered by having MatScatter use<br>
INSERT_VALUES when the destination index set is a permutation.<br>
<br>
I would vote for checking for a square matrix, and otherwise use the expensive form.<br>
<br>
<br>
Also, what is the ssh url for the repo?  I've tried some variations on<br>
ssh://<a href="http://petsc@petsc.cs.iit.edu/petsc/petsc-dev" target="_blank">petsc@petsc.cs.iit.edu/petsc/petsc-dev</a> without success.<br>
<br>
It is ssh://<a href="http://petsc@petsc.cs.iit.edu//hg/petsc/petsc-dev" target="_blank">petsc@petsc.cs.iit.edu//hg/petsc/petsc-dev</a><br>
<br>
  Matt<br>
<br>
<br>
Jed<br>
-- <br>
What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.<br>
-- Norbert Wiener<br>
</blockquote>
<br>
</div></div></blockquote></div><br><br clear="all"><br>-- <br>What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.<br>
-- Norbert Wiener<br>