[petsc-users] IS Invert Not Permutation

[Florian Lindner]

Am 28.11.18 um 18:12 schrieb Matthew Knepley:
> On Wed, Nov 28, 2018 at 12:06 PM Florian Lindner <mailinglists at xgm.de <mailto:mailinglists at xgm.de>> wrote:
> 
>     Hey,
> 
>     thanks for your quick reply!
> 
>     As far as I understand how MatSetLocalToGlobalMapping works, there is no way to create such a mapping that can be used for MatSetValuesLocal?
> 
>     What I need is basically
> 
>     MatSetValuesLocal row/col argument = 3 / 4 / 6 maps to local row/col 0 / 1 / 2
> 
> 
> You are inverting the meaning of MatSetValuesLocal(). It takes in local indices (0, 1, 2) and sets the
> values into a global matrix using global indices (3, 4, 6) which it gets from a compact table.
> 
> It you have global indices (3, 4, 6) use MatSetValues() directly.
> 
> Why are you trying to convert global indices to local indices, which requires a search.

Hey,

we are building a matrix for RBF interpolation. Given two sets of vertices, this requires to evaluate a function on every pair of vertices. The mesh is decomposed. The vertices have a globally unique, but not necessarily continous globalIndex.

A rank 1 has the local portion of the mesh with the globalIndizes = {3, 4, 6}. 

| (1,1) (1,2) | (1,3) (1,4) (1,6)   (1,7) (1,8) (off rank row)
| (2,1) (2,2) | (2,3) (2,4) (2,6)   (2,7) (3,8) (off rank row)
--
  (3,1) (3,2) | (3,3) (3,4) (3,6) | (3,7) (3,8)
  (4,1) (4,2) | (4,3) (4,4) (4,6) | (4,7) (4,8)
  (6,1) (6,2) | (6,3) (6,4) (6,6) | (6,7) (6,8)
--
(more off rank rows)

(marked processor boundaries with | and --)

Going through local vertices, I could just count and use MatSetValuesLocal without a mapping and set the rows as {0, 1, 2}.

But for the cols I need a global enumeration. Just counting is not an option. While the locally owned vertices are disjoint on each rank, each rank has some additional (halo like) vertices vertices. It does not own these (it doesn't build rows for them) but uses them building rows for its own vertices). The size of these halos depend on the basis function used. I can't use globalIndex directly, as it may contain holes and so on.

For that I go through all local vertices and collect their globalIndex in an std::vector. My original implementation was to create an IS from the globalIndizes, ISInvertPermutation, ISAllGather and ISLocalToGlobalMappingCreateIS.

My problem is, that InvertPermutation does not allow the IS to have holes in it.

Following your advice, I use an AO instead and translate the rowIdx (that maybe I can omit) and the colIdx using AOApplicationToPetsc before calling MatSetValues, see below.

Any thoughts on that matter are deeply appreciated!

Best,
Florian



For that, my original implementation was to go ro

> 
>   Thanks,
> 
>     Matt
>  
> 
>     I seems to me, that if I set an index set with MatSetLocalToGlobalMapping like {3, 4, 6} it does the opposite. i.e. maps local {0, 1, 2] to global {3, 4, 6}
> 
>     Therefore I probably need to create my own translation of {3, 4, 6} to {0, 1, 2}
> 
>     A first sketch looks like that
> 
>     mapping[rank] = {3, 4, 6}
> 
>     for (auto i : mapping[rank]) {
>         for (auto j : mapping[rank]) {
>           cout << "Setting " << i << ", " << j << endl;
>           PetscScalar v[] = {i*10.0 + j};
>           PetscInt rows[] = {i};
>           PetscInt cols[] = {j};
>           AOApplicationToPetsc(ao, 1, rows);
>           AOApplicationToPetsc(ao, 1, cols);
>           cout << "Really setting " << rows[0] << ", " << cols[0] << endl << endl;
>           ierr = MatSetValues(matrix, 1, rows, 1, cols, v, INSERT_VALUES); CHKERRQ(ierr);
>         }
>       }
> 
>     which seems to work so far.
> 
>     Any comments from your side are appreciated!
> 
>     Of course, I should set as much values as possible using MatSetValues...
> 
>     Best,
>     Florian
> 
> 
>     Am 27.11.18 um 16:33 schrieb Matthew Knepley:
>     > On Tue, Nov 27, 2018 at 10:17 AM Florian Lindner via petsc-users <petsc-users at mcs.anl.gov <mailto:petsc-users at mcs.anl.gov> <mailto:petsc-users at mcs.anl.gov <mailto:petsc-users at mcs.anl.gov>>> wrote:
>     >
>     >     Hello,
>     >
>     >     I have a range of local input data indices that I want to use for row indexing, say { 3, 4, 6}.
>     >
>     >     For that, I create a matrix with a local number of rows of 3 and map the indices {3, 4, 5} to these rows.
>     >
>     >
>     > It seems like you are trying to create a mapping and its inverse. We do that in 
>     >
>     >   https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/AO/AOCreateMapping.html
>     >
>     > Note that this is not really scalable. Usually there is a way to restructure the code to avoid this.
>     >
>     >   Thanks,
>     >
>     >     Matt
>     >  
>     >
>     >     I create an index set:
>     >
>     >     ISCreateGeneral(comm, myIndizes.size(), myIndizes.data(), PETSC_COPY_VALUES, &ISlocal);
>     >     ISSetPermutation(ISlocal);
>     >     ISInvertPermutation(ISlocal, myIndizes.size(), &ISlocalInv);
>     >     ISAllGather(ISlocalInv, &ISglobal); // Gather the IS from all processors
>     >     ISLocalToGlobalMappingCreateIS(ISglobal, &ISmapping); // Make it a mapping
>     >
>     >     MatSetLocalToGlobalMapping(matrix, ISmapping, ISmapping); // Set mapping for rows and cols
>     >
>     >     The InvertPermutation is required because, well, it seems that the direction the mapping stuff works in PETSc.
>     >
>     >     Now I can use MatSetValues local to set rows {3, 4, 5} and actually set the local rows {1, 2, 3}.
>     >
>     >     The problem is that the range of data vertices is not always contiguous, i.e. could be {3, 4, 6}. This seems to be not a permutation of PETSc, therefore ISSetPermutation fails and subsequently ISInvertPermutation.
>     >
>     >     How can I still create such a mapping, so that:
>     >
>     >     3 -> 1
>     >     4 -> 2
>     >     6 -> 6
>     >
>     >     row 5 is unassigned and will never be addressed.
>     >
>     >     Thanks!
>     >     Florian
>     >
>     >
>     >
>     >
>     > --
>     > What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.
>     > -- Norbert Wiener
>     >
>     > https://www.cse.buffalo.edu/~knepley/ <http://www.cse.buffalo.edu/~knepley/>
> 
> 
> 
> -- 
> What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.
> -- Norbert Wiener
> 
> https://www.cse.buffalo.edu/~knepley/ <http://www.cse.buffalo.edu/~knepley/>