[petsc-dev] Problem in VecAXPY

[Matthew Knepley]

On Thu, Jun 17, 2010 at 4:18 PM, Barry Smith <bsmith at mcs.anl.gov> wrote:

>
> On Jun 17, 2010, at 4:07 PM, Matthew Knepley wrote:
>
> Pushed. ML never ever worked right. Ugh.
>
>
>    Is it now converging better?
>

We are going to test right after this session ends.

   Matt


>
>     Barry
>
>
>
>    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
>
>
>


-- 
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/dab5164f/attachment.html>