[petsc-dev] BlockGetIndices and GetBlockIndices

[Smith, Barry F.]

> 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 {
    ierr = PetscMalloc1(bs*n,&jj);CHKERRQ(ierr);
    *idx = jj;
    k    = 0;
    ii   = sub->idx;
    for (i=0; i<n; i++)
      for (j=0; j<bs; j++)
        jj[k++] = bs*ii[i] + j;
  }
  PetscFunctionReturn(0);
}


> 
>>> Again, my goal is just to avoid having to write something like: https://www.mcs.anl.gov/petsc/petsc-dev/src/vec/is/utils/isltog.c.html#line322
>>> If things were done as in ISLocalToGlobalMapping (which I know is not the same at all as IS/ISBlock), I think this could just read:
>>> ISGetBlockIndices(is,&indices)
>>> ISLocalToGlobalMappingCreate(comm,bs,n/bs,indices,PETSC_COPY_VALUES,mapping)
>> 
>> If you want your code to take advantage of the underlying structure of the IS (be it strided, blocked) for optimizations then you HAVE to go through the subclasses by design and nature (for example VecScatterCreate() takes advantage of the structure of the IS as does ISLocalToGlobalMappingCreate()). There is no universal concept of a block size in IS (despite the crap that other people put in using PetscLayout which has nothing to do with ISBlock and is completely buggy.) nor should there be or is there a need for. 
> 
> OK then.
> 
> Thanks,
> Pierre
> 
>> Barry
>> 
>> 
>> 
>> 
>>> 
>>> Thanks,
>>> Pierre
>>> 
>>>> 
>>>>> On Oct 16, 2019, at 2:50 AM, Pierre Jolivet via petsc-dev <petsc-dev at mcs.anl.gov> wrote:
>>>>> 
>>>>> Hello,
>>>>> I’m trying to understand what is the rationale for naming a function ISBlockGetIndices and another ISLocalToGlobalMappingGetBlockIndices (BlockGet vs. GetBlock).
>>>>> Also, it looks to me like the implementation of ISBlockGetIndices is somehow less versatile than ISLocalToGlobalMappingGetBlockIndices.
>>>>> Indeed, I may call ISLocalToGlobalMappingGetBlockIndices with an underlying non-block IS, while I can’t for ISBlockGetIndices, so a check must be performed by the user, e.g., https://www.mcs.anl.gov/petsc/petsc-current/src/vec/is/is/examples/tutorials/ex3.c.html#line58 or this (IMHO) ugly code duplication https://www.mcs.anl.gov/petsc/petsc-dev/src/vec/is/utils/isltog.c.html#line322.
>>>>> 
>>>>> Thoughts and/or comments? Would it make sense to add an ISGetBlockIndices/ISRestoreBlockIndices or would that be too confusing for the user?
>>>>> 
>>>>> Thanks,
>>>>> Pierre