[petsc-dev] PetscObjectStateIncrease in MatMult*

[Barry Smith]

> On Oct 19, 2015, at 5:06 AM, Václav Hapla <vaclav.hapla at vsb.cz> wrote:
> 
> Again, I can see that MatShellSetOperation sets mat->ops->mult pointer directly
> (see http://www.mcs.anl.gov/petsc/petsc-dev/src/mat/impls/shell/shell.c.html line 779).
> to user-specified MatMult unless MatShellUseScaledMethods has been called (e.g. from MatScale_Shell)
> (see http://www.mcs.anl.gov/petsc/petsc-dev/src/mat/impls/shell/shell.c.html line 70).
> 
> So your MatMult_Shell and MatMultTranspose_Shell is completely ignored if one has not scaled its Mat.

  Yes, you are right!

  The code is way to wonky at the moment with the MatShellUseScaledMethods() and I think will not work properly with certain combinations of user calls to set operations. All to save a few runtime checks on pointers during matrix operations? If so it is overly complicated micro-optimization that should be removed.

There are predefined routines MatDiagonalScale_Shell, MatScale_Shell, MatShift_Shell, that the user could override individually later but they set variables in the shell that will be used in places like MatShellShiftAndScale() (called by MatMult_Shell()) even if the user overwrote the MatShift_Shell() to just get the shift value and handle it within the shell provided operations (like the MatMult) which means the code would improperly do two shifts and scales.

  So to fix this we could
1) make the current logic even more complicated to handle all the possible combinations (yuck) 

2) discard the complicated logic and somehow? replace it. 
     a) one way to do this is for PETSc to ALWAYS handle the MatDiagonalScale_Shell, MatScale_Shell, MatShift_Shell and simply have MatSetOperation() generate an error if the user tries to set one of these. 
    b) another way might be to have Shell provide MatDiagonalScale_Shell, MatScale_Shell, MatShift_Shell by default (like a) but if the user overwrites
        any of the three then Shell automatically cleans up any of the data set by them and unsets the other ones from the _Shell format. In other words if
        the user wants to overwrite one of them then they need to overwrite all of them if they want the functionality of the other ones. This is still super messy, for example what if the user does
             MatShellCreate()
             MatScale(a1)
             MatOperationSet(MAT_SCALE)
             MatScale(a2) 

      the information set in a1 is lost and they would get a different answer then expected. Or we could have MatOperationSet(MAT_SCALE) generate an error if 
the MatScale() had previously been called. Seems to require a good amount of careful logic (which is never a desirable feature of software).

  What to do?

  Barry



> 
> Sorry Barry for replying this only to you before.
> 
> Vaclav
> 
> 
> Dne 17.10.2015 v 22:08 Barry Smith napsal(a):
>> https://bitbucket.org/petsc/petsc/pull-requests/378/update-the-output-vector-state-in-pc-and/diff
>> 
>> On Oct 17, 2015, at 11:12 AM, Václav Hapla <vaclav.hapla at vsb.cz> wrote:
>>> Yes, exactly, MATSHELL.
>>> But if I am right, MatShellSetOperation sets mat->ops->mult pointer directly to user-specified MatMult unless MatShellUseScaledMethods has been called.
>>> So user must keep in mind he should call PetscObjectStateIncrease in his MatMult.
>>> I feel like this is something one would intuitively await from the interface function...
>>> It would be less error-prone for future and user-implemented MatTypes.
>>> Vaclav
>>> 
>>> Dne 17.10.2015 v 17:58 Barry Smith napsal(a):
>>>>> On Oct 17, 2015, at 5:10 AM, Václav Hapla <vaclav.hapla at vsb.cz> wrote:
>>>>> 
>>>>> Hello,
>>>>> MatMult(), MatMultTranspose(), MatMultAdd(), MatMultTransposeAdd() and PCApply() belong to most common modifiers of Vec values.
>>>>> Maybe it would be nice if these routines called PetscObjectStateIncrease explicitly because of MATSHELL and such stuff.
>>>>> Thanks,
>>>>> Vaclav
>>>>    I assume you are suggesting this because users of MATSHELL or PCSHELL may not realize they need to call PetscObjectStateIncrease() on the output vector if they do not use VecGe/RestoretArray() to access the vector entries? (VecRestoreArray() automatically increases the state)
>>>> 
>>>>    What if we add it to MatMult_Shell(), MatMultTranspose_Shell(), PCApply_Shell() and PCApplyTranspose_Shell()?  It looks like those are the only ones where the state won't get increased automatically
>>>> 
>>>>    Thanks for the suggestion,
>>>> 
>>>>   Barry
>>>> 
>>>>> -- 
>>>>> 
>>>>> Vaclav Hapla
>>>>> Junior researcher
>>>>> IT4Innovations <http://www.it4i.eu/>
>>>>> 
>>>>> 
> 
>