[petsc-dev] Problem in VecAXPY

[Barry Smith]

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

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-dev/attachments/20100616/84aa5f10/attachment.html>