[petsc-dev] MatAXPY does not increase state

[Jose E. Roman]

El 02/04/2014, a las 18:21, Jed Brown escribió:

> "Jose E. Roman" <jroman at dsic.upv.es> writes:
> 
>> Hi.
>> 
>> We are having problems since the MatStructure flag was removed from KSPSetOperators.
>> https://bitbucket.org/petsc/petsc/commits/b37f9b8
>> 
>> Tracking our problem leads to blame MatAXPY, which does not increase
>> the state of the first Mat argument. The attached patch fixes the
>> problem for us. Is this the right thing to do?
> 
> Sort of.
> 
>> diff --git a/src/mat/impls/aij/seq/aij.c b/src/mat/impls/aij/seq/aij.c
>> index f05b858..91a51d5 100644
>> --- a/src/mat/impls/aij/seq/aij.c
>> +++ b/src/mat/impls/aij/seq/aij.c
>> @@ -2862,6 +2862,7 @@ PetscErrorCode MatAXPY_SeqAIJ(Mat Y,PetscScalar a,Mat X,MatStructure str)
>>     ierr = MatHeaderReplace(Y,B);CHKERRQ(ierr);
>>     ierr = PetscFree(nnz);CHKERRQ(ierr);
>>   }
>> +  ierr = PetscObjectStateIncrease((PetscObject)Y);CHKERRQ(ierr);
>>   PetscFunctionReturn(0);
>> }
> 
> I would rather move this up to the place where memory is actually
> accessed.  Vec typically uses VecGetArray to manage access (and state),
> while Mat usually uses direct access.  There may be many places with
> this issue.
> 
>> diff --git a/src/mat/utils/gcreate.c b/src/mat/utils/gcreate.c
>> index 55665bf..eb8fd76 100644
>> --- a/src/mat/utils/gcreate.c
>> +++ b/src/mat/utils/gcreate.c
>> @@ -336,8 +336,9 @@ PetscErrorCode MatHeaderMerge(Mat A,Mat C)
>> #define __FUNCT__ "MatHeaderReplace"
>> PETSC_EXTERN PetscErrorCode MatHeaderReplace(Mat A,Mat C)
>> {
>> -  PetscErrorCode ierr;
>> -  PetscInt       refct;
>> +  PetscErrorCode   ierr;
>> +  PetscInt         refct;
>> +  PetscObjectState state;
>> 
>>   PetscFunctionBegin;
>>   PetscValidHeaderSpecific(A,MAT_CLASSID,1);
>> @@ -356,9 +357,11 @@ PETSC_EXTERN PetscErrorCode MatHeaderReplace(Mat A,Mat C)
>> 
>>   /* copy C over to A */
>>   refct = ((PetscObject)A)->refct;
>> +  state = ((PetscObject)A)->state;
>>   ierr  = PetscMemcpy(A,C,sizeof(struct _p_Mat));CHKERRQ(ierr);
>> 
>>   ((PetscObject)A)->refct = refct;
>> +  ((PetscObject)A)->state = state;
> 
> Surely this should increment state so that someone holding a reference
> to A can tell that it has changed.
> 
>>   ierr = PetscFree(C);CHKERRQ(ierr);
>>   PetscFunctionReturn(0);

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.

Jose