[petsc-dev] BlockGetIndices and GetBlockIndices

[Pierre Jolivet]

> On 21 Oct 2019, at 3:01 PM, Jed Brown <jed at jedbrown.org> wrote:
> 
> 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.

Ah, OK, got confused by the lack of bs in the ISCreateGeneral constructor.

Thanks,
Pierre

> 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);
> }
>