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>