#include <petscksp.h>

int main(int argc, char **argv)
{

  Mat          C, D, nest[4] = {NULL, NULL, NULL, NULL};
  Vec          left, right;
  PetscScalar *a;
  PetscInt     rStart, rEnd;

  PetscCall(PetscInitialize(&argc, &argv, NULL, NULL));
  PetscCall(MatCreateConstantDiagonal(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, 4, 4, 2., &nest[0]));
  PetscCall(MatCreateConstantDiagonal(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, 4, 4, 1., &nest[3]));
  PetscCall(MatCreateNest(PETSC_COMM_WORLD, 2, NULL, 2, NULL, nest, &C));
  PetscCall(MatViewFromOptions(C, NULL, "-nest_view"));
  PetscCall(MatCreateVecs(C, &left, &right));
  PetscCall(MatConvert(C, MATAIJ, MAT_INITIAL_MATRIX, &D));
  PetscCall(MatViewFromOptions(D, NULL, "-full_view"));

  PetscCall(VecGetOwnershipRange(left, &rStart, &rEnd));
	PetscCall(VecGetArray(left, &a));
  for (PetscInt r = 0; r < rEnd - rStart; ++r) {
    a[r] = rStart + r;
  }
	PetscCall(VecRestoreArray(left, &a));
  PetscCall(VecViewFromOptions(left, NULL, "-left_view"));
	PetscCall(MatMult(C, left, right));
	PetscCall(VecViewFromOptions(right, NULL, "-right_view"));

  PetscCall(VecDestroy(&left));
  PetscCall(VecDestroy(&right));
  PetscCall(MatDestroy(&nest[0]));
  PetscCall(MatDestroy(&nest[3]));
  PetscCall(MatDestroy(&C));
  PetscCall(MatDestroy(&D));
  PetscCall(PetscFinalize());
  return 0;
}