[petsc-users] Inverse local to global mapping?

[Florian Lindner]

Hello,

I want to build a global vector (and matrix) from local data. The local data has a global index, which can be non-contiguous i.e. global index (0, 5, 2) is on rank 0, (1, 4, 3) is on rank 1. To keep all local data on the local part of vector I need a mapping:

Application -> PETSc
(0, 5, 2)   -> (0, 1, 2)
(1, 4, 3)   -> (3, 4, 5)

It can happen that rank 1 also need to touch vector[5] (application ordering) i.e. vector[1] (petsc ordering) on rank 0.

I can produce a mapping using index sets:

ierr = ISCreateGeneral(PETSC_COMM_WORLD, indizes.size(), indizes.data(), PETSC_COPY_VALUES, &ISlocal); CHKERRV(ierr);
ierr = ISAllGather(ISlocal, &ISglobal); CHKERRV(ierr); // Gather the IS from all processors
ierr = ISLocalToGlobalMappingCreateIS(ISglobal, &ISmapping); CHKERRV(ierr); // Make it a mapping

I construct std::vector indizes by walking though all local indizes. This way I have a mapping ISmapping that looks like that:

0 0
1 5
2 2
3 1
4 4
5 3

Now I can use that mapping: VecSetLocalToGlobalMapping(vec, mapping);

but this does not work like it should (for me). If I use VecSetValueLocal to set element 5 it gets mapped to 3, element 1 gets mapped to 5.
I want to have it reversed, accessing element 3 gets mapped to 5, element 1 to 3.

This loop:

  for (size_t i = 0; i < indizes.size(); i++) {
    ierr = VecSetValueLocal(vec, indizes[i], i*10, INSERT_VALUES); CHKERRQ(ierr);
  }

should produce a vector = [0, 10, 20, 30, 40, 50] and not [0, 50, 20, 10, 40, 30].

Are AO (application orderings) more suited for that? But afaik I can not produce a ISLocalToGlobalMapping from an AO that can be used with VecSetValueLocal, can I?

Is there a way to inverse a mapping? To produce:

0 0
1 3
2 2
3 5 
4 4
5 1


Or am I down the wrong path?

Thanks,
Florian

Example code:


#include <iostream>
#include <vector>
#include <mpi.h>
#include "petscmat.h"   
#include "petscviewer.h"
#include "petscksp.h"

using namespace std;

int main(int argc, char **argv)
{
  PetscInitialize(&argc, &argv, "", NULL);
  PetscErrorCode ierr;
  Vec vec;
  VecCreate(PETSC_COMM_SELF, &vec);

  std::vector<int> indizes = {0, 5, 2, 1, 4, 3};

  IS is;
  ISLocalToGlobalMapping mapping;
  ierr = ISCreateGeneral(PETSC_COMM_WORLD, indizes.size(),indizes.data(), PETSC_COPY_VALUES, &is); CHKERRQ(ierr);
  ierr = ISLocalToGlobalMappingCreateIS(is, &mapping); CHKERRQ(ierr);
  ierr = VecCreate(PETSC_COMM_WORLD, &vec); CHKERRQ(ierr);
  ierr = VecSetLocalToGlobalMapping(vec, mapping); CHKERRQ(ierr);

  ierr = ISLocalToGlobalMappingView(mapping, PETSC_VIEWER_STDOUT_SELF); CHKERRQ(ierr);

  ierr = VecSetSizes(vec, PETSC_DECIDE, indizes.size()); CHKERRQ(ierr);
  ierr = VecSetFromOptions(vec); CHKERRQ(ierr);

  for (size_t i = 0; i < indizes.size(); i++) {
    ierr = VecSetValueLocal(vec, indizes[i], i*10, INSERT_VALUES); CHKERRQ(ierr);
  }

  ierr = VecView(vec, PETSC_VIEWER_STDOUT_WORLD); CHKERRQ(ierr);

  PetscFinalize();
  return 0;
}