[petsc-dev] Problem in VecAXPY

[Matthew Knepley]

Pushed. ML never ever worked right. Ugh.

   Matt

On Wed, Jun 16, 2010 at 7:24 PM, Barry Smith <bsmith at mcs.anl.gov> wrote:

>
> On Jun 16, 2010, at 9:38 AM, Matthew Knepley wrote:
>
> Barry,
>
>   You inserted a check in VecAXPY:563 (which Jed corrected)
>
>  jed at 16195
>
>    563
>
>   if (x == y) SETERRQ(((PetscObject)x)->comm,PETSC_ERR_ARG_IDN,"x and y cannot be the same vector");
>
>
> However, this seems to make no sense in light of mg.c:53
>
>     ierr =
> MatInterpolateAdd(mglevels->interpolate,mgc->x,mglevels->x,mglevels->x);CHKERRQ(ierr);
>
> which will eventually call this with y == w, and has stack
>
> 0]PETSC ERROR: VecAXPY() line 563 in src/vec/vec/interface/rvector.c
> [0]PETSC ERROR: MatMultAdd_ML() line 179 in src/ksp/pc/impls/ml/ml.c
>
>
>    Is this code not completely wrong?
>
>   ML_Operator_Apply(shell->mlmat,x_length,xarray,y_length,yarray);
>   ierr = VecRestoreArray(x,&xarray);CHKERRQ(ierr);
>   ierr = VecRestoreArray(y,&yarray);CHKERRQ(ierr);
>   ierr = VecAXPY(y,1.0,w);CHKERRQ(ierr);
>
>    ML_Operator() overwrites y. Now if w is the same as y that means it
> overwrites w now the VecAXPY(y,1.0,w) does not have the correct w (because y
> values are in w) so it ends up being 2*y which is not waht we want since it
> is suppose to add the original input values of w in.
>
>    I think it needs to have two cases:
>
>    if (w != y) {
>       ML_Operator_Apply(shell->mlmat,x_length,xarray,y_length,yarray);
>       ierr = VecRestoreArray(x,&xarray);CHKERRQ(ierr);
>      ierr = VecRestoreArray(y,&yarray);CHKERRQ(ierr);
>       ierr = VecAXPY(y,1.0,w);CHKERRQ(ierr);
>   } else {
>         copy w into a work vector then do
>       ierr = VecRestoreArray(x,&xarray);CHKERRQ(ierr);
>      ierr = VecRestoreArray(y,&yarray);CHKERRQ(ierr);
>       ierr = VecAXPY(y,1.0,wworkvector);CHKERRQ(ierr);
>   }
>
>   Or do I totally misunderstand the code? Since it appears you are using
> this routine can you fix it?
>
>
>
> Can I remove this check?
>
>
>    No, I really like requiring these various arguments to be different with
> the vector routines. It is one consistent model that keeps the code simple
> and easy to understand. And we don't really need to support having it be the
> same vector. If the only place in ALL of PETSc where it is the same vector
> is actually in a bug I think that demonstrates there is no use case for the
> vectors to be the same.
>
>    As someone else pointed out the BLAS call cannot have the two arguments
> be the same (because it is Fortran) so we've have to handle that case with
> additional ugly checks if we did support the same vector.
>
>    Barry
>
>
>   Thanks,
>
>       Matt
>
> --
> 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/20100617/b6b2c925/attachment.html>