[petsc-dev] Problem in VecAXPY

[Lisandro Dalcin]

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<http://petsc/petsc-dev/annotate/c5b7ef1f9e2d/src/vec/vec/interface/rvector.c#l563>
>
>    563 <#1294132b5614bbcc_l563>
>
>   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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-dev/attachments/20100616/64c4b0c3/attachment.html>