[petsc-dev] BlockGetIndices and GetBlockIndices

[Jed Brown]

Pierre Jolivet via petsc-dev <petsc-dev at mcs.anl.gov> writes:

>> On 21 Oct 2019, at 7:52 AM, Smith, Barry F. <bsmith at mcs.anl.gov> wrote:
>> 
>> 
>> 
>>> On Oct 21, 2019, at 12:23 AM, Pierre Jolivet <pierre.jolivet at enseeiht.fr> wrote:
>>> 
>>> 
>>> 
>>>> On 21 Oct 2019, at 7:11 AM, Smith, Barry F. <bsmith at mcs.anl.gov> wrote:
>>>> 
>>>> 
>>>> 
>>>>> On Oct 20, 2019, at 11:52 PM, Pierre Jolivet <pierre.jolivet at enseeiht.fr> wrote:
>>>>> 
>>>>> 
>>>>> 
>>>>>> On 21 Oct 2019, at 6:42 AM, Smith, Barry F. <bsmith at mcs.anl.gov> wrote:
>>>>>> 
>>>>>> Could you provide a use case where you want to access/have a block size of a IS that is not an ISBlock? 
>>>>> 
>>>>> In the end, all I really want is get access to the underlying is->data->idx without having to worry about the subclass of is.
>>>>> I don’t have such a use case, but I don’t think this is really related to what I want to achieve (or maybe it is…).
>>>> 
>>>> ISGetIndices()
>>> 
>>> Not true for ISBlock with bs > 1.
>> 
>>  Certainly suppose to be, is there a bug?
>> 
>> static PetscErrorCode ISGetIndices_Block(IS in,const PetscInt *idx[])
>> {
>>  IS_Block       *sub = (IS_Block*)in->data;
>>  PetscErrorCode ierr;
>>  PetscInt       i,j,k,bs,n,*ii,*jj;
>> 
>>  PetscFunctionBegin;
>>  ierr = PetscLayoutGetBlockSize(in->map, &bs);CHKERRQ(ierr);
>> 
>> Dang, there is that stupid layout stuff again. Who put this crap in. 
>> 
>>  ierr = PetscLayoutGetLocalSize(in->map, &n);CHKERRQ(ierr);
>>  n   /= bs;
>>  if (bs == 1) *idx = sub->idx;
>>  else {
>
> There it is, I don’t want this if/else. ISGetBlockIndices would have been a function always returning sub->idx.

Your code still can't skip the branch because ISGeneral can have bs>1.
Block size in IS means "these indices go together" while ISBlock imposes
the further constraint: "indices in each block are contiguous".  So you
can't just take the code you quoted where it returns the raw idx:

if (!isblock) {
  ISGetIndices(is,&indices);
  ISLocalToGlobalMappingCreate(comm,1,n,indices,PETSC_COPY_VALUES,mapping);
  ISRestoreIndices(is,&indices);
} else {
  ISGetBlockSize(is,&bs);
  ISBlockGetIndices(is,&indices);
  ISLocalToGlobalMappingCreate(comm,bs,n/bs,indices,PETSC_COPY_VALUES,mapping);
  ISBlockRestoreIndices(is,&indices);
}


PetscErrorCode  ISGetBlockSize(IS is,PetscInt *size)
{
  PetscErrorCode ierr;

  PetscFunctionBegin;
  ierr = PetscLayoutGetBlockSize(is->map, size);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}