[petsc-dev] Problem in VecAXPY

[Lisandro Dalcin]

On 16 June 2010 12:54, Matthew Knepley <knepley at gmail.com> wrote:
>
> That cannot be right because this check was added two days ago.
>

I'm with you. VecAXPY should accept x==y. My point is that in such
case, perhaps we should not actually use an "axpy". Using a 1.0+alpha
"scale"seems both correct and faster. Am I missing something?


>
>    Matt
>
> On Wed, Jun 16, 2010 at 10:40 AM, Lisandro Dalcin <dalcinl at gmail.com> wrote:
>>
>>
>> On 16 June 2010 11:38, Matthew Knepley <knepley at gmail.com> 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
>>> [0]PETSC ERROR: MatMultAdd() line 2041 in src/mat/interface/matrix.c
>>> [0]PETSC ERROR: MatInterpolateAdd() line 6535 in src/mat/interface/matrix.c
>>> [0]PETSC ERROR: PCMGMCycle_Private() line 53 in src/ksp/pc/impls/mg/mg.c
>>> [0]PETSC ERROR: PCMGMCycle_Private() line 50 in src/ksp/pc/impls/mg/mg.c
>>> [0]PETSC ERROR: PCMGMCycle_Private() line 50 in src/ksp/pc/impls/mg/mg.c
>>> [0]PETSC ERROR: PCApply_MG() line 288 in src/ksp/pc/impls/mg/mg.c
>>> [0]PETSC ERROR: PCApply() line 358 in src/ksp/pc/interface/precon.c
>>> [0]PETSC ERROR: KSPSolve_PREONLY() line 27 in src/ksp/ksp/impls/preonly/preonly.c
>>> [0]PETSC ERROR: KSPSolve() line 427 in src/ksp/ksp/interface/itfunc.c
>>> [0]PETSC ERROR: PCApply_FieldSplit() line 530 in src/ksp/pc/impls/fieldsplit/fieldsplit.c
>>> [0]PETSC ERROR: PCApply() line 358 in src/ksp/pc/interface/precon.c
>>> [0]PETSC ERROR: KSPInitialResidual() line 65 in src/ksp/ksp/interface/itres.c
>>> [0]PETSC ERROR: KSPSolve_GMRES() line 236 in src/ksp/ksp/impls/gmres/gmres.c
>>> [0]PETSC ERROR: KSPSolve() line 427 in src/ksp/ksp/interface/itfunc.c
>>> Can I remove this check?
>>>   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
>>
>> I think x==y have to be special-cased (potential issues with arg aliasing in BLAS?). What about:  if (x==y) VecScale(x, 1+alpha)
>>
>> --
>> Lisandro Dalcin
>> ---------------
>> CIMEC (INTEC/CONICET-UNL)
>> Predio CONICET-Santa Fe
>> Colectora RN 168 Km 472, Paraje El Pozo
>> Tel: +54-342-4511594 (ext 1011)
>> Tel/Fax: +54-342-4511169
>>
>
>
>
> --
> 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



--
Lisandro Dalcin
---------------
CIMEC (INTEC/CONICET-UNL)
Predio CONICET-Santa Fe
Colectora RN 168 Km 472, Paraje El Pozo
Tel: +54-342-4511594 (ext 1011)
Tel/Fax: +54-342-4511169