[petsc-users] MatMult inside a for loop

[Matthew Knepley]

On Sun, Feb 8, 2015 at 7:20 PM, Ronal Celaya <ronalcelayavzla at gmail.com>
wrote:

> I know that. I want to have all the vector x replicated in all processes
> and update it in each iteration, so I don't need to communicate the vector
> x each time MatMult() is called.
> I'm not sure I'm making myself clear, sorry
>

1) This is not a scalable strategy.

2) If you know A and all the updates to x locally, why don't you just
compute y directly?

  Thanks,

     Matt


> On Sun, Feb 8, 2015 at 7:52 PM, Barry Smith <bsmith at mcs.anl.gov> wrote:
>
>>
>> > On Feb 8, 2015, at 6:14 PM, Ronal Celaya <ronalcelayavzla at gmail.com>
>> wrote:
>> >
>> > Thank you Barry.
>> > Is there a way to reuse the vector x? I don't want to gather the vector
>> in each iteration, I'd rather replicate the vector x in each process.
>>
>>   I don't understand. With each new matrix vector product there are new
>> values in x (which are updated from other parts of the CG algorithm), these
>> new values need to be communicated in the MatMult() to where they are
>> needed; you can't just reuse the old values in x.
>>
>>   Barry
>>
>> >
>> > Thanks in advance.
>> >
>> > On Sun, Feb 8, 2015 at 7:17 PM, Barry Smith <bsmith at mcs.anl.gov> wrote:
>> >
>> > > On Feb 8, 2015, at 5:41 PM, Ronal Celaya <ronalcelayavzla at gmail.com>
>> wrote:
>> > >
>> > > Hello
>> > > If I have a MatMult operation inside a for loop (e. g. CG algorithm),
>> and the matrix A is MPIAIJ, vector x is gathered to local process in every
>> loop?
>> >
>> >   Yes, internal to MatMult() it calls MatMult_MPIAIJ() which is in
>> src/mat/impls/aij/mpi/mpiaij,c which has the following code:
>> >
>> > PetscErrorCode MatMult_MPIAIJ(Mat A,Vec xx,Vec yy)
>> > {
>> >   Mat_MPIAIJ     *a = (Mat_MPIAIJ*)A->data;
>> >   PetscErrorCode ierr;
>> >   PetscInt       nt;
>> >
>> >   PetscFunctionBegin;
>> >   ierr = VecGetLocalSize(xx,&nt);CHKERRQ(ierr);
>> >   if (nt != A->cmap->n)
>> SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Incompatible partition of A
>> (%D) and xx (%D)",A->cmap->n,nt);
>> >   ierr =
>> VecScatterBegin(a->Mvctx,xx,a->lvec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
>> >   ierr = (*a->A->ops->mult)(a->A,xx,yy);CHKERRQ(ierr);
>> >   ierr =
>> VecScatterEnd(a->Mvctx,xx,a->lvec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
>> >   ierr = (*a->B->ops->multadd)(a->B,a->lvec,yy,yy);CHKERRQ(ierr);
>> >   PetscFunctionReturn(0);
>> >
>> > The needed values of x are communicated in the VecScatterBegin() to
>> VecScatterEnd(). Note only exactly those values needed by each process are
>> communicated in the scatter so not all values are communicated to all
>> processes. Since the matrix is very sparse (normally) only a small
>> percentage of the values need to be communicated.
>> >
>> >   Barry
>> >
>> > >
>> > > I'm sorry for my English.
>> > >
>> > > Regards,
>> > >
>> > > --
>> > > Ronal Celaya
>> >
>> >
>> >
>> >
>> > --
>> > Ronal Celaya
>>
>>
>
>
> --
> Ronal Celaya
>



-- 
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-users/attachments/20150208/d08a14e8/attachment.html>