// Using PETSc MatIS, how to matmult a global IS matrix and a global vector ?
//
// A 3x3 global matrix is built (diag: 1, 2, 1): it's made of 2 overlapping 2x2 local matrix (diag: 1, 1).
//
// How to multiply this matrix with a global vector (filled with 1.) ?
// Expected result => vector such that {proc0: 1, 2 + proc1: 2, 1} => it's not what I get ?!...
//
// ~> g++ -o matISProdMatVec.exe matISProdMatVec.cpp -lpetsc -lm; mpirun -n 2 matISProdMatVec.exe
// ~> mpirun -n 2 matISProdMatVec.exe

#include <iostream>
#include "petsc.h"

int main(int argc,char **argv) {
  PetscInitialize(&argc, &argv, NULL, NULL);
  int size = 0; MPI_Comm_size(MPI_COMM_WORLD, &size); if (size != 2) return 1;
  int rank = 0; MPI_Comm_rank(MPI_COMM_WORLD, &rank);

  PetscInt localSize = 2, globalSize = localSize*2 /*2 MPI*/;
  PetscInt localIdx[2] = {0, 0};
  if (rank == 0) {localIdx[0] = 0; localIdx[1] = 1;}
  else           {localIdx[0] = 1; localIdx[1] = 2;}
  ISLocalToGlobalMapping map;
  ISLocalToGlobalMappingCreate(PETSC_COMM_WORLD, 1, localSize, localIdx, PETSC_COPY_VALUES, &map);
  Mat A;
  MatCreateIS(PETSC_COMM_WORLD, 1, localSize, localSize, globalSize, globalSize, map, map, &A);
  ISLocalToGlobalMappingDestroy(&map);
  MatISSetPreallocation(A, localSize, NULL, localSize, NULL);
  PetscScalar localVal[4] = {1., 0., 0., 1.};
  MatSetValues(A, 2, localIdx, 2, localIdx, localVal, ADD_VALUES);
  MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY); MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
  if (rank == 0) std::cout << std::endl << "matrix A:" << std::endl << std::endl;
  MatView(A, PETSC_VIEWER_STDOUT_WORLD); PetscViewerFlush(PETSC_VIEWER_STDOUT_WORLD); // Diag: 1, 2, 1

  Vec x, y ;
  MatCreateVecs(A, &x, &y);
  VecSet(x, 1.);
  if (rank == 0) std::cout << std::endl << "vector x:" << std::endl << std::endl;
  VecView(x, PETSC_VIEWER_STDOUT_WORLD);
  MatMult(A, x, y);
  if (rank == 0) std::cout << std::endl << "vector y=Ax:" << std::endl << std::endl;
  VecView(y, PETSC_VIEWER_STDOUT_WORLD); // Expected result => proc0: 1, 2 + proc1: 2, 1 => it's not what I get ?!...

  PetscFinalize();
  return 0;
}