[petsc-users] Create identity IS

[Florian Lindner]

Am Montag, 20. April 2015, 14:57:58 schrieb Lisandro Dalcin:
> On 20 April 2015 at 14:35, Florian Lindner <mailinglists at xgm.de> wrote:
> > Hello,
> >
> > I use an index set for a row mapping but don't want to use one for the column mapping.
> >
> > I try to create an idenitity index set to supply to ISLocalToGlobalMappingCreateIS and MatSetLocalToGlobalMapping:
> >
> >   IS is;
> >   ISCreate(PETSC_COMM_WORLD, &is);
> 
> At this point the index set does not have the type set.
> 
> >   ISSetIdentity(is);
> >
> > but I get an segmentation fault at the last line.
> >
> 
> Bad error checking, you should get an error saying that the type is not yet set.
> 
> > What is the best way to create an identity index set?
> >
> 
> I would use ISCreateStride(), and perhaps next ISToGeneral() if you
> really want a ISGENERAL type.

All I want is an identity index set, that I can use with ISLocalToGlobalMappingCreateIS and MatSetLocalToGlobalMapping.

Problem is when I know only the local size, not the global size of the matrix

ISCreateStride(PETSC_COMM_WORLD, nLocal, 0, 1, &is); followed by ISAllGather(is, &isg); gives an IS like that:

0 0
1 1
2 2
3 3
4 0
5 1
6 2
7 3

so 0 to 3 are mapped twice, while the rest is not mapped. Just using the local IS, without ISAllGather gives an Argument out of range at MatSetValuesLocal.

EDIT: Got it now using:

  ISCreateStride(PETSC_COMM_WORLD, A.ownerRange().second - A.ownerRange().first, A.ownerRange().first, 1, &is);
  ISSetIdentity(is);

  ISView(is, PETSC_VIEWER_STDOUT_WORLD);
  ISAllGather(is, &isg);
  if (rank==0)
    ISView(isg, PETSC_VIEWER_STDOUT_SELF);
  
  ISLocalToGlobalMapping mapping;
  ISLocalToGlobalMappingCreateIS(isg, &mapping);
  MatSetLocalToGlobalMapping(A.matrix, mapping, mapping);

based on the owner range.