"Jose E. Roman" <jro...@dsic.upv.es> writes: > I have created a pull request with the change in MatHeaderReplace, > with state+1. Regarding the other changes in MatAXPY, I don't know > where they should go, so I left it unchanged. Hence, the problem still > exists in the case of same and subset nonzero pattern.
Sorry about dropping this. The place where memory is modified should
have the update. Like I would write this for AIJ (if not using an
access interface to manage state).
diff --git i/src/mat/impls/aij/seq/aij.c w/src/mat/impls/aij/seq/aij.c
index 3c8fe3a..dbb12e6 100644
--- i/src/mat/impls/aij/seq/aij.c
+++ w/src/mat/impls/aij/seq/aij.c
@@ -2826,6 +2826,7 @@ PetscErrorCode MatAXPY_SeqAIJ(Mat Y,PetscScalar a,Mat
X,MatStructure str)
PetscScalar alpha = a;
PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&bnz,&alpha,x->a,&one,y->a,&one));
ierr = MatSeqAIJInvalidateDiagonal(Y);CHKERRQ(ierr);
+ ierr = PetscObjectStateIncrease((PetscObject)Y);CHKERRQ(ierr);
} else if (str == SUBSET_NONZERO_PATTERN) { /* nonzeros of X is a subset of
Y's */
if (y->xtoy && y->XtoY != X) {
ierr = PetscFree(y->xtoy);CHKERRQ(ierr);
@@ -2837,6 +2838,7 @@ PetscErrorCode MatAXPY_SeqAIJ(Mat Y,PetscScalar a,Mat
X,MatStructure str)
ierr = PetscObjectReference((PetscObject)X);CHKERRQ(ierr);
}
for (i=0; i<x->nz; i++) y->a[y->xtoy[i]] += a*(x->a[i]);
+ ierr = PetscObjectStateIncrease((PetscObject)Y);CHKERRQ(ierr);
ierr = PetscInfo3(Y,"ratio of nnz(X)/nnz(Y): %D/%D =
%g\n",x->nz,y->nz,(double)((PetscReal)(x->nz)/(y->nz+1)));CHKERRQ(ierr);
} else {
Mat B;
pgpYTaB7IXpd1.pgp
Description: PGP signature
