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>