Re: [petsc-users] Dropping single entries in matrix insertion for block matrices

[Smith, Barry F. via petsc-users]

> On Dec 10, 2018, at 8:30 PM, Jed Brown <j...@jedbrown.org> wrote:
> 
> "Smith, Barry F." <bsm...@mcs.anl.gov> writes:
> 
>>> On Dec 10, 2018, at 7:58 PM, Jed Brown <j...@jedbrown.org> wrote:
>>> 
>>> "Smith, Barry F. via petsc-users" <petsc-users@mcs.anl.gov> writes:
>>> 
>>>>  Well the thing is the indices you pass in are currently always
>>>>  "blocked base", they work automatically for both blocked and non
>>>>  blocked versions of MatSetValuesLocal(). 
>>> 
>>> What do you mean?
>> 
>>   The indices passed to ISLocalToGlobalMappingCreate() are with respect to 
>> the block size of the mapping which has to match the block size of the 
>> matrix.
> 
> You could just create a scalar mapping.
> 
>>> You pass block indices to MatSetValuesBlocked and
>>> scalar indices to MatSetValues.  In MatSetValues_SeqBAIJ, for example
>>> 
>>> for (k=0; k<m; k++) { /* loop over added rows */
>>>   row  = im[k];
>>>   brow = row/bs;
>>> 
>>>>  To support what you want you need to be able to pass in nonblock
>>>>  based version of the indices and (as you noted) keep track of
>>>>  whether you pass in the blocked or non blocked indices so you can
>>>>  apply the scaling or not.
>>>> 
>>>>  Ok, give it a try, it does change the interface a bit but I suppose
>>>>  is worth it since it introduces more functionality.
>>> 
>>> How would the interface change?
>> 
>>    One would now be able to pass a non-blocked ISLocalToGlobalMapping to a 
>> blocked matrix. This is currently not allowed/supported. 
> 
> Is that an interface changing or just not erroring in a circumstance
> where it previously did?

    It is more than not erroring; it could not do what is requested (and it 
printed a useful error message). In the new case where you pass a scalar 
mapping it will not be able to do the MatSetValuesLocalBlocked() while it can 
do MatSetValuesLocal() so it could be slightly confusing. I'm fine with trying 
it as I said.

>  I normally think of interface change meaning
> that previously correct code either behaves differently or fails.